fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

psblas3:

Browse Source 
base/serial/f77/Makefile
 base/serial/f77/cnumbmm.f
 base/serial/f77/dnumbmm.f
 base/serial/f77/smmp.f
 base/serial/f77/snumbmm.f
 base/serial/f77/symbmm.f
 base/serial/f77/znumbmm.f

Split smmp.f
psblas-3.2.0
Salvatore Filippone 11 years ago
parent af8ae96628
commit eb63bf8c4b
7 changed files with 480 additions and 431 deletions
Show Stats Download Patch File Download Diff File

3
base/serial/f77/Makefile
Unescape Escape View File

@ -4,7 +4,8 @@ include ../../../Make.inc
# The object files # The object files
# #
FOBJS = iaxpby.o daxpby.o saxpby.o \ FOBJS = iaxpby.o daxpby.o saxpby.o \
caxpby.o zaxpby.o smmp.o caxpby.o zaxpby.o symbmm.o \
cnumbmm.o dnumbmm.o snumbmm.o znumbmm.o
OBJS=$(FOBJS) OBJS=$(FOBJS)

85
base/serial/f77/cnumbmm.f
Unescape Escape View File

@ -0,0 +1,85 @@
c == =====================================================================
c Sparse Matrix Multiplication Package
c
c Randolph E. Bank and Craig C. Douglas
c
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
c
c Compile this with the following command (or a similar one):
c
c f77 -c -O smmp.f
c
c == =====================================================================
subroutine cnumbmm(n, m, l,
* ia, ja, diaga, a,
* ib, jb, diagb, b,
* ic, jc, diagc, c,
* temp)
c
use psb_const_mod
integer(psb_ipk_) :: ia(*), ja(*), diaga,
* ib(*), jb(*), diagb,
* ic(*), jc(*), diagc
c
complex(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
c
c numeric matrix multiply c=a*b
c
maxlmn = max(l,m,n)
do 10 i = 1,maxlmn
temp(i) = 0.
10 continue
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
c
c c = a*b
c
do 50 i = 1,n
do 30 jj = ia(i),ia(i+1)
c a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) goto 30
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
c 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 20 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
20 continue
30 continue
c c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
c$$$ if (mod(i,100) == 1)
c$$$ + write(psb_err_unit,*)
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do 40 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
40 continue
50 continue
return
end

85
base/serial/f77/dnumbmm.f
Unescape Escape View File

@ -0,0 +1,85 @@
c == =====================================================================
c Sparse Matrix Multiplication Package
c
c Randolph E. Bank and Craig C. Douglas
c
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
c
c Compile this with the following command (or a similar one):
c
c f77 -c -O smmp.f
c
c == =====================================================================
subroutine dnumbmm(n, m, l,
* ia, ja, diaga, a,
* ib, jb, diagb, b,
* ic, jc, diagc, c,
* temp)
use psb_const_mod
c
integer(psb_ipk_) :: ia(*), ja(*), diaga,
* ib(*), jb(*), diagb,
* ic(*), jc(*), diagc
c
real(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
c
c numeric matrix multiply c=a*b
c
maxlmn = max(l,m,n)
do 10 i = 1,maxlmn
temp(i) = 0.
10 continue
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
c
c c = a*b
c
do 50 i = 1,n
do 30 jj = ia(i),ia(i+1)
c a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) goto 30
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
c 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 20 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
20 continue
30 continue
c c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
c$$$ if (mod(i,100) == 1)
c$$$ + write(psb_err_unit,*)
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do 40 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
40 continue
50 continue
return
end

430
base/serial/f77/smmp.f
Unescape Escape View File

@ -1,430 +0,0 @@
c == =====================================================================
c Sparse Matrix Multiplication Package
c
c Randolph E. Bank and Craig C. Douglas
c
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
c
c Compile this with the following command (or a similar one):
c
c f77 -c -O smmp.f
c
c == =====================================================================
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
c
integer(psb_ipk_) :: ia(*), ja(*), diaga,
* ib(*), jb(*), diagb,
* diagc,
* index(*)
integer(psb_ipk_), allocatable :: ic(:),jc(:)
integer(psb_ipk_) :: nze, info
c
c symbolic matrix multiply c=a*b
c
c$$$ open(iunit)
c$$$ write(iunit,*) 'SYMBMM: ',n,m,l
c$$$ write(iunit,*) 'SYMBMM: A:'
c$$$ do i=1,n
c$$$ write(iunit,*) 'Row:',i,' : ',ja(ia(i):ia(i+1)-1)
c$$$ enddo
c$$$ write(iunit,*) 'SYMBMM: B:'
c$$$ do i=1,m
c$$$ write(iunit,*) 'Row:',i,' : ',jb(ib(i):ib(i+1)-1)
c$$$ enddo
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 10 i=1,maxlmn
index(i)=0
10 continue
if (diagc.eq.0) then
ic(1)=1
else
ic(1)=n+2
endif
minlm = min(l,m)
minmn = min(m,n)
c
c main loop
c
do 50 i=1,n
c$$$ write(psb_err_unit,*)
c$$$ 'SYMBMM: 1 loop ',i,n,ia(i),ia(i+1)
istart=-1
length=0
c
c merge row lists
c
do 30 jj=ia(i),ia(i+1)
c a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) goto 30
j = i
else
j=ja(jj)
endif
c 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 20 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
20 continue
30 continue
c
c row i of jc
c
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 40 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
40 continue
call psb_msort(jc(ic(i):ic(i)+length -1))
c$$$ write(iunit,*) length,' : ',jc(ic(i):ic(i)+length-1)
index(i) = 0
50 continue
c$$$ close(iunit)
c$$$ iunit = iunit + 1
c$$$ write(psb_err_unit,*)
c$$$ 'SYMBMM: on exit',ic(n+1)-1,jc(ic(n+1)-1)
return
end
subroutine snumbmm(n, m, l,
* ia, ja, diaga, a,
* ib, jb, diagb, b,
* ic, jc, diagc, c,
* temp)
use psb_const_mod
c
integer(psb_ipk_) :: ia(*), ja(*), diaga,
* ib(*), jb(*), diagb,
* ic(*), jc(*), diagc
c
real(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
c
c numeric matrix multiply c=a*b
c
maxlmn = max(l,m,n)
do 10 i = 1,maxlmn
temp(i) = 0.
10 continue
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
c
c c = a*b
c
do 50 i = 1,n
do 30 jj = ia(i),ia(i+1)
c a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) goto 30
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
c 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 20 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
20 continue
30 continue
c c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
c$$$ if (mod(i,100) == 1)
c$$$ + write(psb_err_unit,*)
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do 40 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
40 continue
50 continue
return
end
subroutine dnumbmm(n, m, l,
* ia, ja, diaga, a,
* ib, jb, diagb, b,
* ic, jc, diagc, c,
* temp)
use psb_const_mod
c
integer(psb_ipk_) :: ia(*), ja(*), diaga,
* ib(*), jb(*), diagb,
* ic(*), jc(*), diagc
c
real(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
c
c numeric matrix multiply c=a*b
c
maxlmn = max(l,m,n)
do 10 i = 1,maxlmn
temp(i) = 0.
10 continue
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
c
c c = a*b
c
do 50 i = 1,n
do 30 jj = ia(i),ia(i+1)
c a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) goto 30
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
c 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 20 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
20 continue
30 continue
c c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
c$$$ if (mod(i,100) == 1)
c$$$ + write(psb_err_unit,*)
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do 40 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
40 continue
50 continue
return
end
subroutine cnumbmm(n, m, l,
* ia, ja, diaga, a,
* ib, jb, diagb, b,
* ic, jc, diagc, c,
* temp)
c
use psb_const_mod
integer(psb_ipk_) :: ia(*), ja(*), diaga,
* ib(*), jb(*), diagb,
* ic(*), jc(*), diagc
c
complex(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
c
c numeric matrix multiply c=a*b
c
maxlmn = max(l,m,n)
do 10 i = 1,maxlmn
temp(i) = 0.
10 continue
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
c
c c = a*b
c
do 50 i = 1,n
do 30 jj = ia(i),ia(i+1)
c a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) goto 30
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
c 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 20 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
20 continue
30 continue
c c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
c$$$ if (mod(i,100) == 1)
c$$$ + write(psb_err_unit,*)
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do 40 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
40 continue
50 continue
return
end
subroutine znumbmm(n, m, l,
* ia, ja, diaga, a,
* ib, jb, diagb, b,
* ic, jc, diagc, c,
* temp)
c
use psb_const_mod
integer(psb_ipk_) :: ia(*), ja(*), diaga,
* ib(*), jb(*), diagb,
* ic(*), jc(*), diagc
c
complex(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
c
c numeric matrix multiply c=a*b
c
maxlmn = max(l,m,n)
do 10 i = 1,maxlmn
temp(i) = 0.
10 continue
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
c
c c = a*b
c
do 50 i = 1,n
do 30 jj = ia(i),ia(i+1)
c a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) goto 30
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
c 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 20 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
20 continue
30 continue
c c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
c$$$ if (mod(i,100) == 1)
c$$$ + write(psb_err_unit,*)
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do 40 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
40 continue
50 continue
return
end

85
base/serial/f77/snumbmm.f
Unescape Escape View File

@ -0,0 +1,85 @@
c == =====================================================================
c Sparse Matrix Multiplication Package
c
c Randolph E. Bank and Craig C. Douglas
c
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
c
c Compile this with the following command (or a similar one):
c
c f77 -c -O smmp.f
c
c == =====================================================================
subroutine snumbmm(n, m, l,
* ia, ja, diaga, a,
* ib, jb, diagb, b,
* ic, jc, diagc, c,
* temp)
use psb_const_mod
c
integer(psb_ipk_) :: ia(*), ja(*), diaga,
* ib(*), jb(*), diagb,
* ic(*), jc(*), diagc
c
real(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
c
c numeric matrix multiply c=a*b
c
maxlmn = max(l,m,n)
do 10 i = 1,maxlmn
temp(i) = 0.
10 continue
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
c
c c = a*b
c
do 50 i = 1,n
do 30 jj = ia(i),ia(i+1)
c a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) goto 30
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
c 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 20 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
20 continue
30 continue
c c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
c$$$ if (mod(i,100) == 1)
c$$$ + write(psb_err_unit,*)
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do 40 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
40 continue
50 continue
return
end

138
base/serial/f77/symbmm.f
Unescape Escape View File

@ -0,0 +1,138 @@
c == =====================================================================
c Sparse Matrix Multiplication Package
c
c Randolph E. Bank and Craig C. Douglas
c
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
c
c Compile this with the following command (or a similar one):
c
c f77 -c -O smmp.f
c
c == =====================================================================
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
c
integer(psb_ipk_) :: ia(*), ja(*), diaga,
* ib(*), jb(*), diagb,
* diagc,
* index(*)
integer(psb_ipk_), allocatable :: ic(:),jc(:)
integer(psb_ipk_) :: nze, info
c
c symbolic matrix multiply c=a*b
c
c$$$ open(iunit)
c$$$ write(iunit,*) 'SYMBMM: ',n,m,l
c$$$ write(iunit,*) 'SYMBMM: A:'
c$$$ do i=1,n
c$$$ write(iunit,*) 'Row:',i,' : ',ja(ia(i):ia(i+1)-1)
c$$$ enddo
c$$$ write(iunit,*) 'SYMBMM: B:'
c$$$ do i=1,m
c$$$ write(iunit,*) 'Row:',i,' : ',jb(ib(i):ib(i+1)-1)
c$$$ enddo
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 10 i=1,maxlmn
index(i)=0
10 continue
if (diagc.eq.0) then
ic(1)=1
else
ic(1)=n+2
endif
minlm = min(l,m)
minmn = min(m,n)
c
c main loop
c
do 50 i=1,n
c$$$ write(psb_err_unit,*)
c$$$ 'SYMBMM: 1 loop ',i,n,ia(i),ia(i+1)
istart=-1
length=0
c
c merge row lists
c
do 30 jj=ia(i),ia(i+1)
c a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) goto 30
j = i
else
j=ja(jj)
endif
c 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 20 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
20 continue
30 continue
c
c row i of jc
c
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 40 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
40 continue
call psb_msort(jc(ic(i):ic(i)+length -1))
c$$$ write(iunit,*) length,' : ',jc(ic(i):ic(i)+length-1)
index(i) = 0
50 continue
c$$$ close(iunit)
c$$$ iunit = iunit + 1
c$$$ write(psb_err_unit,*)
c$$$ 'SYMBMM: on exit',ic(n+1)-1,jc(ic(n+1)-1)
return
end

85
base/serial/f77/znumbmm.f
Unescape Escape View File

@ -0,0 +1,85 @@
c == =====================================================================
c Sparse Matrix Multiplication Package
c
c Randolph E. Bank and Craig C. Douglas
c
c na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
c
c Compile this with the following command (or a similar one):
c
c f77 -c -O smmp.f
c
c == =====================================================================
subroutine znumbmm(n, m, l,
* ia, ja, diaga, a,
* ib, jb, diagb, b,
* ic, jc, diagc, c,
* temp)
c
use psb_const_mod
integer(psb_ipk_) :: ia(*), ja(*), diaga,
* ib(*), jb(*), diagb,
* ic(*), jc(*), diagc
c
complex(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
c
c numeric matrix multiply c=a*b
c
maxlmn = max(l,m,n)
do 10 i = 1,maxlmn
temp(i) = 0.
10 continue
minlm = min(l,m)
minln = min(l,n)
minmn = min(m,n)
c
c c = a*b
c
do 50 i = 1,n
do 30 jj = ia(i),ia(i+1)
c a = d + ...
if (jj.eq.ia(i+1)) then
if (diaga.eq.0 .or. i.gt.minmn) goto 30
j = i
ajj = a(i)
else
j=ja(jj)
ajj = a(jj)
endif
c 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 20 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
20 continue
30 continue
c c = d + ...
if (diagc.eq.1 .and. i.le.minln) then
c(i) = temp(i)
temp(i) = 0.
endif
c$$$ if (mod(i,100) == 1)
c$$$ + write(psb_err_unit,*)
c$$$ ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
do 40 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
40 continue
50 continue
return
end
Write Preview
Loading…
Cancel
Save