不完全cholesky分解实现

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