不完全cholesky分解实现

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