不完全cholesky分解实现

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