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