不完全cholesky分解实现

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