fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
repackage
oacc_loloum
development
psblas-bgmres
repack-csga
dev-openmp
cuda-multivect
repack-ovrlp
repack-llvm
repack-nvid
repack-track
nond-rep
psblas-anderson
non-diag
rename-cuda
tcud
maint-3.8.1
anderson
repack-precuda
master
omp-walther
sp3mm-interface
omp-threadsafe
lambdaI
maint-3.8.0
tspmm
maint-3.7.1
remotebuild
split_sum
maint-3.7.0
remap-coarse
newG2L
new-context
implement-ainv
mat-allocate
unify_aggr_bld
pizdaint-runs
psblas-3.6-maint
merge-paraggr-newops
fnd_owner
merge-paraggr
psblas-paraggr
psblas-3.6-openmp
scr-persistent-collective
new-parstruct
3.6-recursive
sphalo-a2av
MixedI8
psblas-3.5-maint
ILmat
scsr
waitall
test-fcg
CGS2
scrasmussen/test-restructuring
psblas3-mcbind
psblas-caf-xp
psblas-caf
psblas-caf-ext
psblas-caf-csr
trunk
psblas3-caf
psblas-3.4-maint
psblas3-dense
psblas3-newdesc
psblas3-pattern
psblas-3.3-maint
psblas3-submodules
psblas3-accel
psblas-3.3.1-1
tamestoy
psblas-3.2.0
psblas-3.2.0-maint
psblas-3.1-maint
psblas3-openmp
psblas-testmv
psblas3-final
psblas3-type-indexed
psblas-3.0-maint
v3.5.0
3.7.0-rc1
V3.5.2-2
V3.6.0
V3.6.0-RC2
V3.6.1-1
V3.6.1-2
V3.6.1-3
V3.6.1-4
V3.7.0-1
V3.7.0.1
V3.7.1
v3.5.0-1
v3.5.0-2
v3.5.0-3
v3.5.1
v3.5.1-1
v3.5.2
v3.5.2-1
v3.6.0-1
v3.6.0-rc1
v3.6.1
v3.6.1-rc1
v3.6.1-rc2
v3.7.0
v3.7.0.2
v3.8.0
v3.8.0-1
v3.8.0-2
v3.8.1
v3.8.1-2
v3.8.1-rc1-1
v3.8.1-rc1-2
v3.8.1-rc1-3
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '025350a361'
${ noResults }
psblas3/base/serial/smmp.f90

493 lines
13 KiB
Fortran
Raw Blame History

!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
! Original code adapted from:
! == =====================================================================
! Sparse Matrix Multiplication Package
!
! Randolph E. Bank and Craig C. Douglas
!
! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
!
! Compile this with the following command (or a similar one):
!
! f77 -c -O smmp.f
!
! == =====================================================================
subroutine symbmm(n, m, l, ia, ja, diaga, ib, jb, diagb,&
& ic, jc, diagc, index)
use psb_const_mod
use psb_realloc_mod
use psb_sort_mod, only: psb_msort
implicit none
!
integer(psb_ipk_) :: n,m,l
integer(psb_ipk_) :: ia(*), ja(*), diaga, &
& ib(*), jb(*), diagb, diagc, index(*)
integer(psb_ipk_), allocatable :: ic(:),jc(:)
integer(psb_ipk_) :: nze, info,i,j,k,jj
integer(psb_ipk_) :: minlm,minln,minmn,maxlmn, istart,length
!
! symbolic matrix multiply c=a*b
!
if (size(ic) < n+1) then
write(psb_err_unit,*)&
& 'Called realloc in SYMBMM '
call psb_realloc(n+1,ic,info)
if (info /= psb_success_) then
write(psb_err_unit,*)&
& 'realloc failed in SYMBMM ',info
end if
endif
maxlmn = max(l,m,n)
do i=1,maxlmn
index(i)=0
end do
if (diagc.eq.0) then
ic(1)=1
else
ic(1)=n+2
endif
minlm = min(l,m)
minmn = min(m,n)
!
! main loop
!
do i=1,n
istart=-1
length=0
!
! merge row lists
!
rowi: do jj=ia(i),ia(i+1)
! a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
j = i
else
j=ja(jj)
endif
! b = d + ...
if (index(j).eq.0 .and. diagb.eq.1 .and. j.le.minlm)then
index(j)=istart
istart=j
length=length+1
endif
if ((j<1).or.(j>m)) then
write(psb_err_unit,*)&
& ' SymbMM: Problem with A ',i,jj,j,m
endif
do k=ib(j),ib(j+1)-1
if ((jb(k)<1).or.(jb(k)>maxlmn)) then
write(psb_err_unit,*)&
& 'Problem in SYMBMM 1:',j,k,jb(k),maxlmn
else
if(index(jb(k)).eq.0) then
index(jb(k))=istart
istart=jb(k)
length=length+1
endif
endif
end do
end do rowi
!
! row i of jc
!
if (diagc.eq.1 .and. index(i).ne.0) length = length - 1
ic(i+1)=ic(i)+length
if (ic(i+1) > size(jc)) then
if (n > (2*i)) then
nze = max(ic(i+1), ic(i)*((n+i-1)/i))
else
nze = max(ic(i+1), nint((dble(ic(i))*(dble(n)/i))) )
endif
call psb_realloc(nze,jc,info)
end if
do j= ic(i),ic(i+1)-1
if (diagc.eq.1 .and. istart.eq.i) then
istart = index(istart)
index(i) = 0
endif
jc(j)=istart
istart=index(istart)
index(jc(j))=0
end do
call psb_msort(jc(ic(i):ic(i)+length -1))
index(i) = 0
end do
return
end subroutine symbmm
! == =====================================================================
! Sparse Matrix Multiplication Package
!
! Randolph E. Bank and Craig C. Douglas
!
! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
!
! Compile this with the following command (or a similar one):
!
! f77 -c -O smmp.f
!
! == =====================================================================
subroutine cnumbmm(n, m, l, ia, ja, diaga, a, ib, jb, diagb, b,&
& ic, jc, diagc, c, temp)
!
use psb_const_mod
implicit none
integer(psb_ipk_) :: n,m,l
integer(psb_ipk_) :: ia(*), ja(*), diaga,&
& ib(*), jb(*), diagb, ic(*), jc(*), diagc
!
complex(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
integer(psb_ipk_) :: i,j,k,jj,minlm,minln,minmn,maxlmn
!
! numeric matrix multiply c=a*b
!
maxlmn = max(l,m,n)
do i = 1,maxlmn
temp(i) = 0.
end do
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
!
! c = a*b
!
do i = 1,n
rowi: do jj = ia(i),ia(i+1)
! a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
! b = d + ...
if (diagb.eq.1 .and. j.le.minlm) &
& temp(j) = temp(j) + ajj * b(j)
if ((j<1).or.(j>m)) then
write(psb_err_unit,*)&
& ' NUMBMM: Problem with A ',i,jj,j,m
endif
do k = ib(j),ib(j+1)-1
if((jb(k)<1).or. (jb(k) > maxlmn)) then
write(psb_err_unit,*)&
& ' NUMBMM: jb problem',j,k,jb(k),maxlmn
else
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
endif
end do
end do rowi
! c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
!$$$ if (mod(i,100) == 1)
!$$$ + write(psb_err_unit,*)
!$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do j = ic(i),ic(i+1)-1
if((jc(j)<1).or. (jc(j) > maxlmn)) then
write(psb_err_unit,*)&
& ' NUMBMM: output problem',i,j,jc(j),maxlmn
else
c(j) = temp(jc(j))
temp(jc(j)) = 0.
endif
end do
end do
return
end subroutine cnumbmm
! == =====================================================================
! Sparse Matrix Multiplication Package
!
! Randolph E. Bank and Craig C. Douglas
!
! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
!
! Compile this with the following command (or a similar one):
!
! f77 -c -O smmp.f
!
! == =====================================================================
subroutine dnumbmm(n, m, l, ia, ja, diaga, a, ib, jb, diagb, b,&
& ic, jc, diagc, c, temp)
use psb_const_mod
implicit none
!
integer(psb_ipk_) :: n,m,l
integer(psb_ipk_) :: ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& ic(*), jc(*), diagc
!
real(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
integer(psb_ipk_) :: i,j,k,jj,minlm,minln,minmn,maxlmn
!
! numeric matrix multiply c=a*b
!
maxlmn = max(l,m,n)
do i = 1,maxlmn
temp(i) = 0.
end do
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
!
! c = a*b
!
do i = 1,n
rowi: do jj = ia(i),ia(i+1)
! a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
! b = d + ...
if (diagb.eq.1 .and. j.le.minlm) &
& temp(j) = temp(j) + ajj * b(j)
if ((j<1).or.(j>m)) then
write(psb_err_unit,*)&
& ' NUMBMM: Problem with A ',i,jj,j,m
endif
do k = ib(j),ib(j+1)-1
if((jb(k)<1).or. (jb(k) > maxlmn)) then
write(psb_err_unit,*)&
& ' NUMBMM: jb problem',j,k,jb(k),maxlmn
else
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
endif
end do
end do rowi
! c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
!$$$ if (mod(i,100) == 1)
!$$$ + write(psb_err_unit,*)
!$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do j = ic(i),ic(i+1)-1
if((jc(j)<1).or. (jc(j) > maxlmn)) then
write(psb_err_unit,*)&
& ' NUMBMM: output problem',i,j,jc(j),maxlmn
else
c(j) = temp(jc(j))
temp(jc(j)) = 0.
endif
end do
end do
return
end subroutine dnumbmm
! == =====================================================================
! Sparse Matrix Multiplication Package
!
! Randolph E. Bank and Craig C. Douglas
!
! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
!
! Compile this with the following command (or a similar one):
!
! f77 -c -O smmp.f
!
! == =====================================================================
subroutine snumbmm(n, m, l, ia, ja, diaga, a, ib, jb, diagb, b,&
& ic, jc, diagc, c, temp)
use psb_const_mod
implicit none
!
integer(psb_ipk_) :: n,m,l
integer(psb_ipk_) :: ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& ic(*), jc(*), diagc
!
real(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
integer(psb_ipk_) :: i,j,k,jj,minlm,minln,minmn,maxlmn
!
! numeric matrix multiply c=a*b
!
maxlmn = max(l,m,n)
do i = 1,maxlmn
temp(i) = 0.
end do
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
!
! c = a*b
!
do i = 1,n
rowi: do jj = ia(i),ia(i+1)
! a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
! b = d + ...
if (diagb.eq.1 .and. j.le.minlm) &
& temp(j) = temp(j) + ajj * b(j)
if ((j<1).or.(j>m)) then
write(psb_err_unit,*)&
& ' NUMBMM: Problem with A ',i,jj,j,m
endif
do k = ib(j),ib(j+1)-1
if((jb(k)<1).or. (jb(k) > maxlmn)) then
write(psb_err_unit,*)&
& ' NUMBMM: jb problem',j,k,jb(k),maxlmn
else
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
endif
end do
end do rowi
! c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
!$$$ if (mod(i,100) == 1)
!$$$ + write(psb_err_unit,*)
!$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do j = ic(i),ic(i+1)-1
if((jc(j)<1).or. (jc(j) > maxlmn)) then
write(psb_err_unit,*)&
& ' NUMBMM: output problem',i,j,jc(j),maxlmn
else
c(j) = temp(jc(j))
temp(jc(j)) = 0.
endif
end do
end do
return
end subroutine snumbmm
! == =====================================================================
! Sparse Matrix Multiplication Package
!
! Randolph E. Bank and Craig C. Douglas
!
! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
!
! Compile this with the following command (or a similar one):
!
! f77 -c -O smmp.f
!
! == =====================================================================
subroutine znumbmm(n, m, l, ia, ja, diaga, a, ib, jb, diagb, b,&
& ic, jc, diagc, c, temp)
!
use psb_const_mod
implicit none
integer(psb_ipk_) :: n,m,l
integer(psb_ipk_) :: ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
& ic(*), jc(*), diagc
!
complex(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
integer(psb_ipk_) :: i,j,k,jj,minlm,minln,minmn,maxlmn
!
! numeric matrix multiply c=a*b
!
maxlmn = max(l,m,n)
do i = 1,maxlmn
temp(i) = 0.
end do
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
!
! c = a*b
!
do i = 1,n
rowi: do jj = ia(i),ia(i+1)
! a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
! b = d + ...
if (diagb.eq.1 .and. j.le.minlm) &
& temp(j) = temp(j) + ajj * b(j)
if ((j<1).or.(j>m)) then
write(psb_err_unit,*)&
& ' NUMBMM: Problem with A ',i,jj,j,m
endif
do k = ib(j),ib(j+1)-1
if((jb(k)<1).or. (jb(k) > maxlmn)) then
write(psb_err_unit,*)&
& ' NUMBMM: jb problem',j,k,jb(k),maxlmn
else
temp(jb(k)) = temp(jb(k)) + ajj * b(k)
endif
end do
end do rowi
! c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
do j = ic(i),ic(i+1)-1
if((jc(j)<1).or. (jc(j) > maxlmn)) then
write(psb_err_unit,*)&
& ' NUMBMM: output problem',i,j,jc(j),maxlmn
else
c(j) = temp(jc(j))
temp(jc(j)) = 0.
endif
end do
end do
return
end subroutine znumbmm
Reference in New Issue View Git Blame Copy Permalink