求一个用于稀疏矩阵求解的ICCG的fortran程序
rt，谢谢！ http://www.netlib.org/linalg/iccg回复 #2 风花雪月 的帖子
打不开啊 我这里能够正常打开 implicit real(a-h,o-z)real a(200000),af(200000),b(28000),x(28000),
$ anrm,xnrm,resd,g(28000),p(28000),t1(28000),t2(28000)
integer ia(28001),ja(200000),nit(3)
real d,c2,a2,cf
common /tomm/ a,af,ia,ja
common /area2/ d,c2,a2,cf
common /order/ iorflg
data nit/1,5,10/
c
c this is a main program for running the incomplete LU
c conjugate gradient for solving linear systems.
c the matrix need not be symmetric.
c
c ilucg this flag if 1 will cause incomplete factoization
c if 0 no factorization just straight cg
c
c write(6,*)' start in main'
c
c call trapov(100)
loca = 200000
1 continue
c write(6,*) 'input ineum5,ineum6,irot,iorflg'
c write(6,*)'set irot<0 to stop'
c read(5,*) ineum5,ineum6,irot,iorflg
c if( irot .lt. 0 ) then
c write(6,*)' end test'
c stop
c end if
c write(6,*) 'ineum5,ineum6,irot,iorflg',ineum5,ineum6,irot,iorflg
c write(6,*)'input mesh cells for x, y, and z'
c read(5,*)imx,imy,imz
c write(6,*)'inner, middle, and outer fill-in integers'
c read(5,*)ifli,iflm,iflo
cc write(6,*)'ising (normally zero)'
cc read(5,*)ising
ineum5 = 1
ineum6 = 1
irot = 1
iorflg = 2
imx = 20
imy = 20
imz = 20
cc if (ineum5 .ne. 1 .or. ineum6 .ne. 1) ising=0
ising = 0
cc imx = 5
cc imy = 5
cc imz = 5
n = imx*imy*imz
c write(6,*) ' *********** problem type *************'
c if( ineum5 .eq. 0 ) then
c write(6,*)' ***** dir on bottom *****'
c else
c write(6,*)' ***** neuman on bottom *****'
c endif
c if( ineum6 .eq. 0 ) then
c write(6,*)' ***** dir on top *****'
c else
c write(6,*)' ***** neuman on top *****'
c endif
c if( ising .eq. 0 ) then
c write(6,*)' ***** regular *****'
c else
c write(6,*)' ***** singular if neuman on top and bottom *****'
c endif
c write(6,*) ' ***** mesh cells in x dir',imx,' *****'
c write(6,*) ' ***** mesh cells in y dir',imy,' *****'
c write(6,*) ' ***** mesh cells in z dir',imz,' *****'
c
c ifli is the number of fillin stripes allowed for each of the
c inner stripes; if zero no fillin.
c iflm is the number of fillin stripes allowed for each of the
c middle stripes; if zero no fillin.
c iflo is the number of fillin stripes allowed for each of the
c outer stripes; if zero no fillin here.
c inum5 if 1 then neumann on bottom
c if 0 then dirchlet on bottom
c inum6 if 1 then neumann on top
c if 0 then dirchlet on top
c irotif 0 then no rot flow field
c if 1 then rot flow field
c ising if zero regular
c if 1 then signular if neumann top and bottom
c imx # of mesh cells in the x direction
c imy # of mesh cells in the y direction
c imz # of mesh cells in the z direction
cc ifli = 0
cc iflm = 0
cc iflo = 0
cc write(6,*)' call prbtyp'
call prbtyp(ineum5,ineum6,ising,imx,imy,imz,0)
cc write(6,*)' call matgen'
ti1 = second(ii)
call matgen(ifli,iflm,iflo,irot,nonz,ia,ja,a,b)
c write(6,*)'n,nonz,ia(n+1)'
c write(6,*)n,nonz,ia(n+1)
cc write(6,*)' after matgen'
ti2 = second(ii) - ti1
if( ia(n+1) .gt. loca ) then
write(6,*)' not enough storage',loca,ia(n+1)
stop
endif
c stop
job = 0
anrm = abs(a(isamax(ia(n+1)-1,a,1)))
call scopy(n,b,1,t2,1)
maxits = n
tol = 1.0e-6
write(6,101) ti2
101 format(/20x,'time to generate matrix = ',5x,1pe10.3)
do 888 ics = 1,3
ti3 = second(ii)
do 777 mm = 1, nit(ics)
call scopy(n,t2,1,x,1)
call scopy(n,t2,1,b,1)
call iccglu(a,n,ia,ja,af,b,x,g,p,t1,tol,maxits,its,job,info)
c write(6,*)'info should be zeroinfo = ',info
job = 1
777 continue
ti4 = second(ii) - ti3
call scopy(n,t2,1,b,1)
call cgres(a,n,ia,ja,x,b,p)
resd = abs(p(isamax(n,p,1)))
xnrm = abs(x(isamax(n,x,1)))
c write(6,*)'anrm,xnrm,resd',anrm,xnrm,resd
sres = resd/(anrm*xnrm)
nonzr = ia(n+1) - 1
c ti5 = ti2 + ti4
c write(6,100) sres,ti5,ti2,ti4,nonzr,its
c 100 format(20x,'scaled residual = ',13x,1pe10.3,
c 1 /20x,'total time = ',18x,1pe10.3,
c 2 /20x,'time to generate matrix = ',5x,1pe10.3,
c 3 /20x,'time for cg iterations = ',6x,1pe10.3,
c 4 /20x,'number of nonzero elements = ',2x,i10,
c 5 /20x,'number of cg iterations = ',5x,i10)
write(6,102) nit(ics),its,sres,ti4
102 format(20x,'results for ',i3,' sets of ',i4,' cg iterations ',
1 'for each set',
2 //20x,'scaled residual = ',13x,1pe10.3,
3 /20x,'time for cg iterations = ',6x,1pe10.3,//)
888 continue
c go to 1
stop
end
subroutine matgen(ifli,iflm,iflo,irot,nonz,ia,ja,a,f)
implicit real(a-h,o-z)
integer ia(1),ja(1),ifli,iflo
real a(1)
integer bwi,bwo
real f(28000),acof0(28000),acof1(28000),acof2(28000),
# acof3(28000),
# acof4(28000),acof5(28000),acof6(28000),gl(1000),gu(1000)
real cvx(30,30,30),cvy(30,30,30),cvz(30,30,30)
#,omg(30,30,30)
common/coef/acof0,acof1,acof2,acof3,acof4,acof5,acof6
common /order/ iorflg
c ineum5=1 neuman on bot;=0 dir on bot.
c ineum6=1 neuman on top;=0 dir on top.
c irot = 1, rotational flow field
c irot = 0, non-rotational flow field
flx=1.e0
fly=1.e0
flz=1.e0
c write(6,*)' call prbtyp'
call prbtyp(ineum5,ineum6,ising,imx,jmx,kmx,1)
c write(6,*)' out of prbtyp in matgen'
id=1
idp1=id+1
iptflg=0
go to (3,5,7),iorflg
3 continue
c write(6,*)' continue 3'
bwo=imx*jmx
bwi=imx
ni=1
nj=bwi
nk=bwo
nc=-bwo-bwi
nib=1
njb=imx
ncb=-imx
c (ijk) (5310246)
go to 9
5 continue
c write(6,*)' continue 5'
c (kij) (3150624)
bwo=kmx*imx
bwi=kmx
ni=bwi
nj=bwo
nk=1
nc=-bwo-bwi
nib=1
njb=imx
ncb=-imx
go to 9
7 continue
c write(6,*)' continue 7'
c (jki) (1530462)
bwo=jmx*kmx
bwi=jmx
ni=bwo
nj=1
nk=bwi
nc=-bwo-bwi
nib=jmx
njb=1
ncb=-jmx
9 continue
c write(6,*)' continue 9'
c write(6,*)'imx,jmx,kmx',imx,jmx,kmx
dx=flx/float(imx)
dy=fly/float(jmx)
dz=flz/float(kmx)
rdx2=1.e0/dx**2
rdy2=1.e0/dy**2
rdz2=1.e0/dz**2
mx=imx*jmx*kmx
imxm1=imx-1
jmxm1=jmx-1
kmxm1=kmx-1
nonz=3*mx-2+2*(mx-bwi+ifli)+2*(mx-bwo+iflo)
nonzt=nonz
c write(6,*)' call inits'
call inits(a,mx,ia,ja)
if(ifli .ge. bwi-1) go to 220
if(iflo .ge. bwo-bwi) go to 220
do 10 m=1,mx
f(m)=0.e0
10continue
c write(6,*)' continue 10'
do 11 i=1,imx
do 11 j=1,jmx
do 11 k=1,kmx
xi=(float(i)-0.5e0)*dx
yj=(float(j)-0.5e0)*dy
zk=(float(k)-0.5e0)*dz
facx = 1.0e0
facy = 1.0e0
if( irot .eq. 0 ) go to 12
facx = xi - .05e0
facy = yj - .05e0
12 continue
dum=32.0e0*xi*(1.e0-xi)*yj*(1.e0-yj)*zk
dum1=2.0e0*xi*yj*zk**2
dum=25.e0*dum
dum1=2.0e0*dum1
cvx(i,j,k)=dum*facx
cvy(i,j,k)=dum*facy
cvz(i,j,k)=dum1
omg(i,j,k)=0.e0
dum2 = (xi*xi)*yj*zk
m = i*ni + j*nj + k*nk + nc
f(m) = dum2
11 continue
c write(6,*)' continue 11'
do 15 j=1,jmx
do 15 i=1,imx
mb=i*nib+j*njb+ncb
gl(mb)=1.e0
gu(mb)=2.0e0
15continue
c write(6,*)' continue 15'
do 20 m=1,mx
acof0(m)=0.e0
acof1(m)=0.e0
acof2(m)=0.e0
acof3(m)=0.e0
acof4(m)=0.e0
acof5(m)=0.e0
acof6(m)=0.e0
20continue
c write(6,*)' continue 20'
if(imx .lt. 3) go to 45
if(jmx .lt. 3) go to 45
if(kmx .lt. 3) go to 45
do 30 j=2,jmxm1
do 30 i=2,imxm1
do 30 k=2,kmxm1
m=i*ni+j*nj+k*nk+nc
acof0(m)=2.0e0*(rdx2+rdy2+rdz2)+omg(i,j,k)
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
30continue
c write(6,*)' continue 30'
45continue
c write(6,*)' continue 45'
if((jmx.lt.3).or.(kmx.lt.3)) go to 55
do 50 j=2,jmxm1
do 50 k=2,kmxm1
i=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+2.0e0*(rdy2+rdz2)+omg(i,j,k)-0.5e0*cvx(i,j,k)/dx
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
i=imx
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+2.0e0*(rdy2+rdz2)+omg(i,j,k) +0.5e0*cvx(i,j,k)/dx
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
50continue
c write(6,*)' continue 50'
55continue
c write(6,*)' continue 55'
if((imx.lt.3).or.(kmx.lt.3)) go to 65
do 60 i=2,imxm1
do 60 k=2,kmxm1
j=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdy2+ 2.0e0*(rdx2+rdz2)+omg(i,j,k) -0.5e0*cvy(i,j,k)/dy
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
j=jmx
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdy2+2.0e0*(rdx2+rdz2)+omg(i,j,k)+0.5e0*cvy(i,j,k)/dy
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
60continue
c write(6,*)' continue 60'
65continue
c write(6,*)' continue 65'
if((imx.lt.3).or.(jmx.lt.3)) go to 75
do 70 i=2,imxm1
do 70 j=2,jmxm1
k=1
m=(j-1)*bwo+(i-1)*bwi +k
acof0(m)=2.0e0*(rdx2+rdy2)+3.0e0*rdz2 +omg(i,j,k)
#+cvz(i,j,k)/(2.0e0*dz)
if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) go to 71
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
71 continue
c write(6,*)' continue 71'
k=kmx
m=i*ni+j*nj+k*nk+nc
acof0(m)=2.0e0*(rdx2+rdy2)+3.0e0*rdz2+omg(i,j,k)
#-cvz(i,j,k)/(2.0e0*dz)
if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) go to 72
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
72 continue
c write(6,*)' continue 72'
70continue
c write(6,*)' continue 70'
75continue
c write(6,*)' continue 75'
if(kmx.lt.3) go to 85
do 80 k=2,kmxm1
i=1
j=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+2.0e0*rdz2+omg(i,j,k)
#-0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
j=jmx
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+2.0e0*rdz2+omg(i,j,k)
#-0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
i=imx
m=(j-1)*bwo+(i-1)*bwi +k
acof0(m)=rdx2+rdy2+2.0e0*rdz2+omg(i,j,k)
#+0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
acof6(m)=-rdz2+ 0.5e0*cvz(i,j,k)/dz
j=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+2.0e0*rdz2+omg(i,j,k)
#+0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
80continue
c write(6,*)' continue 80'
85continue
c write(6,*)' continue 85'
if(imx.lt.3)go to 95
do 90 i=2,imxm1
k=1
j=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=2.0e0*rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#-0.5e0*cvy(i,j,k)/dy+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) go to 91
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
91 continue
c write(6,*)' continue 91'
j=jmx
m=i*ni+j*nj+k*nk+nc
acof0(m)=2.0e0*rdx2 +rdy2 +3.0e0*rdz2+omg(i,j,k)
#+0.5e0*cvy(i,j,k)/dy+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) go to 92
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2+ cvz(i,j,k)/dz)*gl(mb)
92 continue
c write(6,*)' continue 92'
k=kmx
m=i*ni+j*nj+k*nk+nc
acof0(m)=2.0e0*rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#+0.5e0*cvy(i,j,k)/dy-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) go to 93
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
93 continue
c write(6,*)' continue 93'
j=1
m=(j-1)*bwo+(i-1)*bwi +k
acof0(m)=2.0e0*rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#-0.5e0*cvy(i,j,k)/dy-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) go to 94
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
94 continue
c write(6,*)' continue 94'
90continue
c write(6,*)' continue 90'
95continue
c write(6,*)' continue 95'
if(jmx.lt.3)go to 105
do 100 j=2,jmxm1
k=1
i=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+ 2.0e0*rdy2+3.0e0*rdz2+omg(i,j,k)
# -0.5e0*cvx(i,j,k)/dx+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) go to 101
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
101 continue
c write(6,*)' continue 101'
i=imx
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+2.0e0*rdy2+3.0e0*rdz2+omg(i,j,k)
# +0.5e0*cvx(i,j,k)/dx+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) go to 102
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
102 continue
c write(6,*)' continue 102'
k=kmx
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+2.0e0*rdy2+3.0e0*rdz2+omg(i,j,k)
# +0.5e0*cvx(i,j,k)/dx-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) go to 103
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
103 continue
c write(6,*)' continue 103'
i=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+2.0e0*rdy2+3.0e0*rdz2+omg(i,j,k)
# -0.5e0*cvx(i,j,k)/dx-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) go to 104
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
104 continue
c write(6,*)' continue 104'
100continue
c write(6,*)' continue 100'
105continue
c write(6,*)' continue 105'
i=1
j=1
k=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#-0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
# +0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dx
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) go to 106
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
106 continue
c write(6,*)' continue 106'
k=kmx
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#-0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
# -0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
if (ineum6. eq. 1) go to 107
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
107 continue
c write(6,*)' continue 107'
i=1
j=jmx
k=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#-0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
# +0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) go to 108
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
108 continue
c write(6,*)' continue 108'
k=kmx
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#-0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
# -0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
acof2(m)=-rdx2+0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) go to 109
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
109 continue
c write(6,*)' continue 109'
i=imx
j=1
k=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#+0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
# +0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) go to 111
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
111 continue
c write(6,*)' continue 111'
k=kmx
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#+0.5e0*cvx(i,j,k)/dx-0.5e0*cvy(i,j,k)/dy
# -0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof4(m)=-rdy2+0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) go to 112
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
112 continue
c write(6,*)' continue 112'
i=imx
j=jmx
k=1
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#+0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
# +0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 -cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof6(m)=-rdz2+0.5e0*cvz(i,j,k)/dz
if (ineum5 .eq. 1) go to 113
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2+cvz(i,j,k)/dz)*gl(mb)
113 continue
c write(6,*)' continue 113'
k=kmx
m=i*ni+j*nj+k*nk+nc
acof0(m)=rdx2+rdy2+3.0e0*rdz2+omg(i,j,k)
#+0.5e0*cvx(i,j,k)/dx+0.5e0*cvy(i,j,k)/dy
# -0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) acof0(m) = acof0(m) - 2.0e0*rdz2 +cvz(i,j,k)/dz
acof1(m)=-rdx2-0.5e0*cvx(i,j,k)/dx
acof3(m)=-rdy2-0.5e0*cvy(i,j,k)/dy
acof5(m)=-rdz2-0.5e0*cvz(i,j,k)/dz
if (ineum6 .eq. 1) go to 114
mb=i*nib+j*njb+ncb
f(m)=f(m)+(2.0e0*rdz2-cvz(i,j,k)/dz)*gu(mb)
114 continue
c write(6,*)' continue 114'
if(ineum5*ineum6 .eq. 0) go to 117
if( ising .eq. 1 ) go to 117
i=1
j=1
k=1
m=i*ni+j*nj+k*nk+nc
acof2(m)=0.e0
acof4(m)=0.e0
acof6(m)=0.e0
f(m)=0.e0
i=2
m=i*ni+j*nj+k*nk+nc
acof1(m)=0.e0
i=1
j=2
m=i*ni+j*nj+k*nk+nc
acof3(m)=0.e0
j=1
k=2
m=i*ni+j*nj+k*nk+nc
acof5(m)=0.e0
117 continue
c write(6,*)' continue 117'
cload coef
c write(6,*)'load matrix'
do 150 i=1+bwo,mx
j=i-bwo
t=acof3(i)
call put(t,a,mx,ia,ja,i,j)
150 continue
c write(6,*)' continue 150'
do 155 i=1+bwi,mx
j=i-bwi
t=acof1(i)
call put(t,a,mx,ia,ja,i,j)
155 continue
c write(6,*)' continue 155'
do 160 i=2,mx
j=i-1
t=acof5(i)
call put(t,a,mx,ia,ja,i,j)
160 continue
c write(6,*)' continue 160'
do 165 i=1,mx
j=i
t=acof0(i)
call put(t,a,mx,ia,ja,i,j)
165 continue
c write(6,*)' continue 165'
do 170 i=1,mx-1
j=i+1
t=acof6(i)
call put(t,a,mx,ia,ja,i,j)
170 continue
c write(6,*)' continue 170'
do 175 i=1,mx-bwi
j=i+bwi
t=acof2(i)
call put(t,a,mx,ia,ja,i,j)
175 continue
c write(6,*)' continue 175'
do 180 i=1,mx-bwo
j=i+bwo
t=acof4(i)
call put(t,a,mx,ia,ja,i,j)
180 continue
c write(6,*)' continue 180'
c load zeros for fill-in stripes
if(iflo .eq. 0) go to 195
do 192 k = 1,iflo
do 185 i=1+bwo-k,mx
j=i-bwo+k
t=0.0e0
call put(t,a,mx,ia,ja,i,j)
185 continue
c write(6,*)' continue 185'
do 190 i=1,mx-bwo+k
j=i+bwo-k
t=0.0e0
call put(t,a,mx,ia,ja,i,j)
190 continue
192 continue
c write(6,*)' continue 190'
195 continue
c write(6,*)' continue 195'
if(ifli .eq. 0) go to 210
do 208 k = 1,ifli
do 200 i=1+bwi-k,mx
j=i-bwi+k
t=0.0e0
call put(t,a,mx,ia,ja,i,j)
200 continue
c write(6,*)' continue 200'
do 205 i=1,mx-bwi+k
j=i+bwi-k
t=0.0e0
call put(t,a,mx,ia,ja,i,j)
205 continue
208 continue
c write(6,*)' continue 205'
210 continue
c write(6,*)' continue 210'
if( iflm .eq. 0 ) go to 219
do 217 k = 1,iflm
do 212 i = 1+bwi+k,mx
j = i - bwi - k
call put(0.0e0,a,mx,ia,ja,i,j)
212 continue
do 213 i = 1,mx-bwi-k
j = i + bwi + k
call put(0.0e0,a,mx,ia,ja,i,j)
213 continue
217 continue
219 continue
c have now loaded the a,ia,ja arrays
c write(6,*) 'formula,actual',nonzt,ia(mx+1)
go to 215
220write(6,*) 'ifli or iflo too large'
215continue
c write(6,*)' continue 215'
return
end
subroutine prbtyp(ineum5,ineum6,ising,imx,imy,imz,iflg)
integer ineum5,ineum6,ising,imx,imy,imz,iflg
integer n5,n6,isg,ix,iy,iz
cineum5=0,dir on bot; =1,neuman on bot
cineum6=0,dir on top; =1,neuman on top
cising =0,regular; =1 singular if neuman on top & bot
cimx = # mesh cells in x dir
cimy = # mesh cells in y dir
cimz = # mesh cells in z dir
ciflg =0 sets local values;=1 returns local values
c
c
if(iflg .eq. 0) go to 10
c write(6,*)' in prbtyp ix,iy,iz recall',ix,iy,iz
ineum5=n5
ineum6=n6
ising=isg
imx=ix
imy=iy
imz=iz
c write(6,*)' in prbtyp recall imx,imy,imz',imx,imy,imz
go to 20
10 continue
c write(6,*)' continue #'
c write(6,*)' in prbtyp init imx,imy,imz',imx,imy,imz
n5=ineum5
n6=ineum6
isg=ising
ix=imx
iy=imy
iz=imz
c write(6,*) ' *********** problem type *************'
c if( ineum5 .eq. 0 ) then
c write(6,*)' ***** dir on bottom *****'
c else
c write(6,*)' ***** neuman on bottom *****'
c endif
c if( ineum6 .eq. 0 ) then
c write(6,*)' ***** dir on top *****'
c else
c write(6,*)' ***** neuman on top *****'
c endif
c if( ising .eq. 0 ) then
c write(6,*)' ***** regular *****'
c else
c write(6,*)' ***** singular if neuman on top and bottom *****'
c endif
c write(6,*) ' ***** mesh cells in x dir',imx,' *****'
c write(6,*) ' ***** mesh cells in y dir',imy,' *****'
c write(6,*) ' ***** mesh cells in z dir',imz,' *****'
20 continue
c write(6,*)' continue #'
return
end
subroutine iccglu(a,n,ia,ja,af,b,x,r,p,s,tol,maxits,its,job,info)
c
integer n,ia(1),ja(1),maxits,its,job,info
real a(1),af(1),x(1),r(1),p(1),s(1),b(1),tol
c
c this routine performs preconditioned conjugate gradient on a
c sparse matrix. the preconditioner is an incomplete lu of
c the matrix.
c
c on entry:
c
c areal()
c contains the elements of the matrix.
c
c ninteger
c is the order of the matrix.
c
c ia integer(n+1)
c contains pointers to the start of the rows in the arrays
c a and ja.
c
c ja integer()
c contains the column location of the corresponding elements
c in the array a.
c
c b real (n)
c contains the right hand side.
c
c x real (n)
c contains an estimate of the solution, the closer it
c is to the true solution the faster tthe method will
c converge.
c
c tol real
c is the accuracy desired in the solution.
c
c maxits integer
c is the maximun number of iterations to be taken
c by the routine.
c
c job integer
c is a flag to signal if incomplete factorization
c aready exits in array af.
c job = 0perform incomplete factorization
c job = 1skip incomplete factorization
c
c on output
c
c af real ()
c contains the incomplete factorization of the matrix
c contained in the array a.
c
c x real (n)
c contains the solution.
c
c r,p,s real (n)
c these are scratch work arrays.
c
c its integer
c contains the number of iterations need to converge.
c
c info integer
c signals if normal termination.
c info = 0 method converged in its iterations
c info = 1 method converged, but exit occurred because
c residual norm was less than sqrt(rtxmin).
c info = -999 method didnot converge in maxits iterations.
c
c
c the algorithm has the following form.
c
c form incomplete factors l and u
c x(0) <- initial estiate
c r(0) <- b - a*x(0)
c r(0) <- trans(a)*invtrans(l)*invtrans(u)*inv(u)*inv(l)*r(0)
c p(0) <- r(0)
c i <- 0
c while r(i) > tol do
c s <- inv(u)*inv(l)*a*p(i)
c a(i) <- trans(r(i))*r(i)/(trans(s)*s)
c x(i+1) <- x(i) + a(i)*p(i)
c r(i+1) <- r(i) - a(i)*trans(a)*invtrans(l)*invtrans(u)*s
c b(i) <- trans(r(i+1))*r(i+1)/(trans(r(i))*r(i))
c p(i+1) <- r(i+1) + b(i)*p(i)
c i <- i + 1
c end
c
real ai,bi,rowold,rownew,xnrm,anrm
real sdot
real rtxmax,rtxmin
common /gear14/ rtxmax,rtxmin
data rtxmin/1.0e-16/
info = 0
anrm = abs(a(isamax(ia(n+1)-1,a,1)))
c
c form incomplete factors l and u
c
if( job .ne. 0 ) go to 5
call scopy(ia(n+1)-1,a,1,af,1)
call lu(af,n,ia,ja)
5 continue
c
c r(0) <- b - a*x(0)
c
call cgres(a,n,ia,ja,x,b,r)
c
c r(0) <- trans(a)*invtrans(l)*invtrans(u)*inv(u)*inv(l)*r(0)
c
c inv(u)*inv(l)*r(0)
c
call scopy(n,r,1,s,1)
call ssol(af,n,ia,ja,s,1)
call ssol(af,n,ia,ja,s,2)
c
c invtrans(l)*invtrans(u)*above
c
call ssol(af,n,ia,ja,s,-2)
call ssol(af,n,ia,ja,s,-1)
c
c trans(a)*above
c
call mmult(a,n,ia,ja,r,s,-1)
c
c p(0) <- r(0)
c
call scopy(n,r,1,p,1)
rowold = sdot(n,r,1,r,1)
i = 0
c
c while r(i) > tol do
c
ai = 1.0d0
10 continue
xnrm = abs(x(isamax(n,x,1)))
cc write(6,*) ' iter residual xnrm ai',i,sqrt(rowold)/(anrm*xnrm),
cc $ xnrm,ai
if( sqrt(rowold) .le. tol*(anrm*xnrm) ) go to 12
cc if (rowold .le. rtxmin) go to 25
13 continue
c
c s <- inv(u)*inv(l)*a*p(i)
c
call mmult(a,n,ia,ja,s,p,1)
call ssol(af,n,ia,ja,s,1)
call ssol(af,n,ia,ja,s,2)
c
c a(i) <- trans(r(i))*r(i)/(trans(s)*s)
c
ai = rowold/sdot(n,s,1,s,1)
c
c x(i+1) <- x(i) + a(i)*p(i)
c
call saxpy(n,ai,p,1,x,1)
c
c
c r(i+1) <- r(i) - a(i)*trans(a)*invtrans(l)*invtrans(u)*s
c
call ssol(af,n,ia,ja,s,-2)
call ssol(af,n,ia,ja,s,-1)
call mmult(a,n,ia,ja,b,s,-1)
call saxpy(n,-ai,b,1,r,1)
c
c b(i) <- trans(r(i+1))*r(i+1)/(trans(r(i))*r(i))
c
rownew = sdot(n,r,1,r,1)
bi = rownew/rowold
c
c p(i+1) <- r(i+1) + b(i)*p(i)
c
call saxpy2(n,bi,p,1,r,1)
rowold = rownew
c
c i <- i + 1
c
i = i + 1
if( i .gt. maxits ) go to 20
go to 10
12 continue
15 continue
its = i
return
20 continue
info = -999
its = maxits
return
c25 continue
c info = 1
c its = i
c return
end
subroutine axpy(a,n,ia,ja,i,js,je,t,y)
c
integer n,ia(1),ja(1),i,js,je
real a(1),y(1)
real t
c
c this routine computes an axpy for a row of a sparse matrix
c with a vector.
c an axpy is a multiple of a row of the matrix is added to the
c vector y.
c
is = ia(i)
ie = ia(i+1) - 1
do 10 ir = is,ie
j = ja(ir)
if( j .gt. je ) go to 20
if( j .lt. js ) go to 10
y(j) = y(j) + t*a(ir)
10 continue
20 continue
return
end
subroutine saxpy2(n,da,dx,incx,dy,incy)
c
c constant times a vector plus a vector.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c
real dx(1),dy(1),da
integer i,incx,incy,m,mp1,n
c
if(n.le.0)return
if (da .eq. 0.0e0) return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dx(iy) = dy(iy) + da*dx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,4)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dx(i) = dy(i) + da*dx(i)
30 continue
if( n .lt. 4 ) return
40 mp1 = m + 1
do 50 i = mp1,n,4
dx(i) = dy(i) + da*dx(i)
dx(i + 1) = dy(i + 1) + da*dx(i + 1)
dx(i + 2) = dy(i + 2) + da*dx(i + 2)
dx(i + 3) = dy(i + 3) + da*dx(i + 3)
50 continue
return
end
real function dot(a,n,ia,ja,i,js,je,b)
c
integer n,ia(1),ja(1),i,js,je
real a(1),b(1)
real t
c
c this routine computes an inner product for a row of a sparse matri
c with a vector.
c
t = 0.0e0
is = ia(i)
ie = ia(i+1) - 1
do 10 ir = is,ie
j = ja(ir)
if( j .gt. je ) go to 20
if( j .lt. js ) go to 10
t = t + a(ir)*b(j)
10 continue
20 continue
dot = t
return
end
subroutine inits(a,n,ia,ja)
c
integer n,ia(1),ja(1)
real a(1)
c
c this routine initializes storage for the sparse
c matrix arrays.
c note: the matrix is initialized to have zeroes
c on the diagonal.
c
do 10 i = 1,n
a(i) = 0.0e0
ja(i) = i
ia(i) = i
10 continue
ia(n+1) = n+1
return
end
subroutine insert(t,a,n,ia,ja,i,j,k)
c
integer n,ia(1),ja(1),i,j,k
integer ip1,np1
real t,a(1)
c
c this routine rearranges the elements in arrays a and ja
c and updates array ia for the new element.
c
cc write(6,1000)i,j,ia(i),ia(i+1),t
1000format('*** from insert i,j,ia(i),ia(i+1),t ',4i5,1x,e15.8)
l1 = k
l2 = ia(n+1) - 1
do 10 lb = l1,l2
l = l2 - lb + l1
a(l+1) = a(l)
ja(l+1) = ja(l)
10 continue
a(k) = t
ja(k) = j
ip1 = i + 1
np1 = n + 1
do 20 l = ip1,np1
ia(l) = ia(l) + 1
20 continue
return
end
logical function locate(a,n,ia,ja,i,j,k)
c
integer n,ia(1),ja(1),i,j
real a(1)
c
c this routine will locate the i,j-th element in the
c sparse matrix structure.
c
is = ia(i)
ie = ia(i+1) - 1
c
c row i is from is to ie in array a.
c
do 10 l = is,ie
if( j .gt. ja(l) ) go to 10
if( j .ne. ja(l) ) go to 5
locate = .true.
k = l
go to 20
5 continue
if( j .ge. ja(l) ) go to 10
locate = .false.
k = 0
go to 20
10 continue
c
c get here if should be at the end of a row.
c
locate = .false.
k = 0
20 continue
return
end
subroutine lu(a,n,ia,ja)
c
logical locate
integer n,ia(1),ja(1)
integer ikp1,kp1
real a(1)
c
c this subroutine does incomplete gaussian elimenation
c on a sparse matrix. the matrix is stored in a sparse
c data structure.
c note: no pivoting is done
c and the factorization is incomplete,
c i.e., only where there exists a storage location
c will the operation take place.
c
do 40 k = 1,n
if( .not. locate(a,n,ia,ja,k,k,l2) ) go to 50
if( a(l2) .eq. 0.0e0 ) go to 50
kp1 = k + 1
do 30 i = kp1,n
if( .not. locate(a,n,ia,ja,i,k,ik) ) go to 25
is = ik
ie = ia(i+1) - 1
kj = l2
ke = ia(k+1) - 1
a(ik) = a(ik)/a(kj)
if( kj .eq. ke ) go to 30
kj = kj + 1
ikp1 = ik + 1
do 20 j = ikp1,ie
10 continue
if( kj .gt. ke ) go to 30
if( ja(kj) .ge. ja(j) ) go to 15
kj = kj + 1
go to 10
15 continue
if( ja(kj) .gt. ja(j) ) go to 20
a(j) = a(j) - a(ik)*a(kj)
20 continue
25 continue
30 continue
40 continue
return
50 continue
cc write(6,*)' value of k and a(k,k)',k,a(l2)
cc write(6,*)' error zero on diagonal'
cc write(6,*)' matrix probably specified incorrectly or'
cc write(6,*)' incomplete factorization produces a singular matrix'
return
end
subroutine mmult(a,n,ia,ja,b,x,job)
c
integer n,ia(1),ja(1),job
real a(1),b(1),x(1)
c
c this routine performs matrix vector multiple
c for a sparse matrix structure
c
c job has the following input meanings:
c job = 1matrix vector multiple
c = -1 matrix transpose vector multiple
c = 2unit lower matrix vector multiple
c = -2 unit lower matrix transpose vector multiple
c = 3upper matrix vector multiple
c = -3 upper matrix transpose vector multiple
real dot
c
c a*x
c
if( job .ne. 1 ) go to 15
do 10 i = 1,n
b(i) = dot(a,n,ia,ja,i,1,n,x)
10 continue
c
c trans(a)*x
c
15 continue
if( job .ne. -1 ) go to 35
do 20 i = 1,n
b(i) = 0.0e0
20 continue
do 30 i = 1,n
call axpy(a,n,ia,ja,i,1,n,x(i),b)
30 continue
c
c l*x when l is unit lower
c
35 continue
if( job .ne. 2 ) go to 45
do 40 i = 1,n
b(i) = x(i) + dot(a,n,ia,ja,i,1,i-1,x)
40 continue
c
c trans(l)*x when l is unit lower
c
45 continue
if( job .ne. -2 ) go to 55
do 49 i = 1,n
b(i) = x(i)
49 continue
do 50 i = 1,n
call axpy(a,n,ia,ja,i,i,n,x(i),b)
50 continue
c
c u*x
c
55 continue
if( job .ne. 3 ) go to 65
do 60 i = 1,n
b(i) = dot(a,n,ia,ja,i,i,n,x)
60 continue
c
c trans(u)*x
c
65 continue
if( job .ne. -3 ) go to 85
do 70 i = 1,n
b(i) = 0.0e0
70 continue
do 80 i = 1,n
call axpy(a,n,ia,ja,i,1,i,x(i),b)
80 continue
85 continue
return
end
subroutine put(t,a,n,ia,ja,i,j)
c
integer n,ia(1),ja(1),i,j
real t,a(1)
c
c this routine will insert an element into the sparse matrix
c structure.
c
is = ia(i)
ie = ia(i+1) - 1
cc write(6,100)i,j,ia(i),ia(i+1)
100 format('*** from put i,j,ia(i),ia(i+1) ',4i7)
c
c row i is from is to ie in array a.
c
do 10 k = is,ie
if( j .gt. ja(k) ) go to 10
if( j .ne. ja(k) ) go to 5
a(k) = t
go to 20
5 continue
if( j .ge. ja(k) ) go to 12
call insert(t,a,n,ia,ja,i,j,k)
go to 20
12 continue
10 continue
c
c get here if should be at the end of a row.
c
k = ie + 1
call insert(t,a,n,ia,ja,i,j,k)
go to 30
20 continue
30 continue
return
end
subroutine cgres(a,n,ia,ja,x,b,r)
c
integer n,ia(1),ja(1)
real a(1),x(1),b(1),r(1)
c
c this routine computes a residual for a*x=b where
c a is in a sparse structure
c
real dot
do 10 i = 1,n
r(i) = b(i) - dot(a,n,ia,ja,i,1,n,x)
10 continue
return
end
subroutine ssol(a,n,ia,ja,b,job)
c
integer n,ia(1),ja(1),job
real a(1),b(1)
c
c this routine solves a system of equations based on a sparse
c matrix date structure.
c the array b contains the right hand side on input
c and on output has the solution
c
c job has the value 1 if l*x = b is to be solved.
c job has the value -1 if trans(l)*x = b is to be solved.
c job has the value 2 if u*x = b is to be solved.
c job has the value -2 if trans(u)*x = b is to be solved.
c
logical locate
real t
real dot
c
c job = 1solvel*x = b
c
if( job .ne. 1 ) go to 15
c
c solve l*y=b
c
do 10 i = 2,n
b(i) = b(i) - dot(a,n,ia,ja,i,1,i-1,b)
10 continue
15 continue
if( job .ne. 2 ) go to 25
c
c solve u*x=y
c
do 20 ib = 1,n
i = n - ib + 1
t = dot(a,n,ia,ja,i,i+1,n,b)
if( .not. locate(a,n,ia,ja,i,i,k) ) go to 30
b(i) = (b(i) - t)/a(k)
20 continue
c
c job = -2solvetrans(u)*x = b
c
25 continue
if( job .ne. -2 ) go to 35
c
c solve trans(u)*y=b
c
do 21 i = 1,n
if( .not. locate(a,n,ia,ja,i,i,k) ) go to 30
b(i) = b(i)/a(k)
call axpy(a,n,ia,ja,i,i+1,n,-b(i),b)
21 continue
c
c solve trans(l)*x=y
c
35 continue
if( job .ne. -1 ) go to 45
do 22 ib = 2,n
i = n - ib + 2
call axpy(a,n,ia,ja,i,1,i-1,-b(i),b)
22 continue
45 continue
return
30 continue
cc write(6,*)' error no diagonal element: from solve'
return
end
subroutine saxpy(n,sa,sx,incx,sy,incy)
c
c constant times a vector plus a vector.
c uses unrolled loop for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c
real sx(1),sy(1),sa
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if (sa .eq. 0.0) return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
sy(iy) = sy(iy) + sa*sx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,4)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
sy(i) = sy(i) + sa*sx(i)
30 continue
if( n .lt. 4 ) return
40 mp1 = m + 1
do 50 i = mp1,n,4
sy(i) = sy(i) + sa*sx(i)
sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
50 continue
return
end
subroutine scopy(n,sx,incx,sy,incy)
c
c copies a vector, x, to a vector, y.
c uses unrolled loops for increments equal to 1.
c jack dongarra, linpack, 3/11/78.
c
real sx(1),sy(1)
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
sy(iy) = sx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,7)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
sy(i) = sx(i)
30 continue
if( n .lt. 7 ) return
40 mp1 = m + 1
do 50 i = mp1,n,7
sy(i) = sx(i)
sy(i + 1) = sx(i + 1)
sy(i + 2) = sx(i + 2)
sy(i + 3) = sx(i + 3)
sy(i + 4) = sx(i + 4)
sy(i + 5) = sx(i + 5)
sy(i + 6) = sx(i + 6)
50 continue
return
end
real function sdot(n,sx,incx,sy,incy)
c
c forms the dot product of two vectors.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c
real sx(1),sy(1),stemp
integer i,incx,incy,ix,iy,m,mp1,n
c
stemp = 0.0e0
sdot = 0.0e0
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
stemp = stemp + sx(ix)*sy(iy)
ix = ix + incx
iy = iy + incy
10 continue
sdot = stemp
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
stemp = stemp + sx(i)*sy(i)
30 continue
if( n .lt. 5 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,5
stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) +
* sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4)
50 continue
60 sdot = stemp
return
end
integer function isamax(n,sx,incx)
c
c finds the index of element having max. absolute value.
c jack dongarra, linpack, 3/11/78.
c
real sx(1),smax
integer i,incx,ix,n
c
isamax = 0
if( n .lt. 1 ) return
isamax = 1
if(n.eq.1)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix = 1
smax = abs(sx(1))
ix = ix + incx
do 10 i = 2,n
if(abs(sx(ix)).le.smax) go to 5
isamax = i
smax = abs(sx(ix))
5 ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 smax = abs(sx(1))
do 30 i = 2,n
if(abs(sx(i)).le.smax) go to 30
isamax = i
smax = abs(sx(i))
30 continue
return
end
页:
[1]