发布时间 : 星期三 文章不完全cholesky分解实现更新完毕开始阅读
求一个用于稀疏矩阵求解的ICCG的fortran程序
implicit real(a-h,o-z)
real a(200000),af(200000),b(28000),x(28000),
$ anrm,xnrm,resd,g(28000),p(28000),t1(28000),t2(28000) integer ia(28001),ja(200000),nit(3) real d,c2,a2,cf
common /tomm/ a,af,ia,ja common /area2/ d,c2,a2,cf common /order/ iorflg data nit/1,5,10/ c
c this is a main program for running the incomplete LU c conjugate gradient for solving linear systems. c the matrix need not be symmetric. c
c ilucg this flag if 1 will cause incomplete factoization c if 0 no factorization just straight cg c
c write(6,*)' start in main' c
c call trapov(100) loca = 200000 1 continue
c write(6,*) 'input ineum5,ineum6,irot,iorflg' c write(6,*)'set irot<0 to stop'
c read(5,*) ineum5,ineum6,irot,iorflg c if( irot .lt. 0 ) then c write(6,*)' end test' c stop c end if c write(6,*)
'ineum5,ineum6,irot,iorflg',ineum5,ineum6,irot,iorflg c write(6,*)'input mesh cells for x, y, and z' c read(5,*)imx,imy,imz
c write(6,*)'inner, middle, and outer fill-in integers' c read(5,*)ifli,iflm,iflo
cc write(6,*)'ising (normally zero)' cc read(5,*)ising ineum5 = 1 ineum6 = 1 irot = 1 iorflg = 2 imx = 20
imy = 20 imz = 20
cc if (ineum5 .ne. 1 .or. ineum6 .ne. 1) ising=0 ising = 0 cc imx = 5 cc imy = 5 cc imz = 5
n = imx*imy*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,' *****' c
c ifli is the number of fillin stripes allowed for each of the c inner stripes; if zero no fillin.
c iflm is the number of fillin stripes allowed for each of the c middle stripes; if zero no fillin.
c iflo is the number of fillin stripes allowed for each of the c outer stripes; if zero no fillin here. c inum5 if 1 then neumann on bottom c if 0 then dirchlet on bottom c inum6 if 1 then neumann on top c if 0 then dirchlet on top c irot if 0 then no rot flow field c if 1 then rot flow field c ising if zero regular
c if 1 then signular if neumann top and bottom c imx # of mesh cells in the x direction c imy # of mesh cells in the y direction
c imz # of mesh cells in the z direction cc ifli = 0 cc iflm = 0 cc iflo = 0
cc write(6,*)' call prbtyp'
call prbtyp(ineum5,ineum6,ising,imx,imy,imz,0) cc write(6,*)' call matgen' ti1 = second(ii)
call matgen(ifli,iflm,iflo,irot,nonz,ia,ja,a,b) c write(6,*)'n,nonz,ia(n+1)' c write(6,*)n,nonz,ia(n+1) cc write(6,*)' after matgen' ti2 = second(ii) - ti1
if( ia(n+1) .gt. loca ) then
write(6,*)' not enough storage',loca,ia(n+1) stop endif c stop job = 0
anrm = abs(a(isamax(ia(n+1)-1,a,1))) call scopy(n,b,1,t2,1) maxits = n tol = 1.0e-6 write(6,101) ti2
101 format(/20x,'time to generate matrix = ',5x,1pe10.3) do 888 ics = 1,3 ti3 = second(ii)
do 777 mm = 1, nit(ics) call scopy(n,t2,1,x,1) call scopy(n,t2,1,b,1)
call iccglu(a,n,ia,ja,af,b,x,g,p,t1,tol,maxits,its,job,info) c write(6,*)'info should be zero info = ',info job = 1 777 continue
ti4 = second(ii) - ti3 call scopy(n,t2,1,b,1)
call cgres(a,n,ia,ja,x,b,p) resd = abs(p(isamax(n,p,1))) xnrm = abs(x(isamax(n,x,1)))
c write(6,*)'anrm,xnrm,resd',anrm,xnrm,resd sres = resd/(anrm*xnrm) nonzr = ia(n+1) - 1 c ti5 = ti2 + ti4
c write(6,100) sres,ti5,ti2,ti4,nonzr,its
c 100 format(20x,'scaled residual = ',13x,1pe10.3, c 1 /20x,'total time = ',18x,1pe10.3,
c 2 /20x,'time to generate matrix = ',5x,1pe10.3, c 3 /20x,'time for cg iterations = ',6x,1pe10.3, c 4 /20x,'number of nonzero elements = ',2x,i10, c 5 /20x,'number of cg iterations = ',5x,i10) write(6,102) nit(ics),its,sres,ti4
102 format(20x,'results for ',i3,' sets of ',i4,' cg iterations ', 1 'for each set',
2 //20x,'scaled residual = ',13x,1pe10.3,
3 /20x,'time for cg iterations = ',6x,1pe10.3,//) 888 continue c go to 1 stop end
subroutine matgen(ifli,iflm,iflo,irot,nonz,ia,ja,a,f) implicit real(a-h,o-z)
integer ia(1),ja(1),ifli,iflo real a(1)
integer bwi,bwo
real f(28000),acof0(28000),acof1(28000),acof2(28000), # acof3(28000),
# acof4(28000),acof5(28000),acof6(28000),gl(1000),gu(1000) real cvx(30,30,30),cvy(30,30,30),cvz(30,30,30) # ,omg(30,30,30)
common/coef/acof0,acof1,acof2,acof3,acof4,acof5,acof6 common /order/ iorflg
c ineum5=1 neuman on bot;=0 dir on bot. c ineum6=1 neuman on top;=0 dir on top. c irot = 1, rotational flow field c irot = 0, non-rotational flow field flx=1.e0 fly=1.e0 flz=1.e0
c write(6,*)' call prbtyp'
call prbtyp(ineum5,ineum6,ising,imx,jmx,kmx,1) c write(6,*)' out of prbtyp in matgen' id=1
idp1=id+1 iptflg=0
go to (3,5,7),iorflg 3 continue
c write(6,*)' continue 3' bwo=imx*jmx