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 irot if 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 zero info = ',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 10 continue 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 15 continue 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 20 continue 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 30 continue c write(6,*)' continue 30' 45 continue 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 50 continue c write(6,*)' continue 50' 55 continue 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 60 continue c write(6,*)' continue 60' 65 continue 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' 70 continue c write(6,*)' continue 70' 75 continue 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 80 continue c write(6,*)' continue 80' 85 continue 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' 90 continue c write(6,*)' continue 90' 95 continue 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' 100 continue c write(6,*)' continue 100' 105 continue 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' c load 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 220 write(6,*) 'ifli or iflo too large' 215 continue 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 c ineum5=0,dir on bot; =1,neuman on bot c ineum6=0,dir on top; =1,neuman on top c ising =0,regular; =1 singular if neuman on top & bot c imx = # mesh cells in x dir c imy = # mesh cells in y dir c imz = # mesh cells in z dir c iflg =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 a real() c contains the elements of the matrix. c c n integer 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 = 0 perform incomplete factorization c job = 1 skip 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 c 25 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 1000 format(' *** 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 = 1 matrix vector multiple c = -1 matrix transpose vector multiple c = 2 unit lower matrix vector multiple c = -2 unit lower matrix transpose vector multiple c = 3 upper 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 = 1 solve l*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 = -2 solve trans(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