发布时间 : 星期三 文章不完全cholesky分解实现更新完毕开始阅读
c is a flag to signal if incomplete factorization c aready exits in array af.
c job = 0 perform incomplete factorization c job = 1 skip incomplete factorization c
c on output c
c af real ()
c contains the incomplete factorization of the matrix c contained in the array a. c
c x real (n)
c contains the solution. c
c r,p,s real (n)
c these are scratch work arrays. c
c its integer
c contains the number of iterations need to converge. c
c info integer
c signals if normal termination.
c info = 0 method converged in its iterations
c info = 1 method converged, but exit occurred because c residual norm was less than sqrt(rtxmin).
c info = -999 method didnot converge in maxits iterations. c c
c the algorithm has the following form. c
c form incomplete factors l and u c x(0) <- initial estiate c r(0) <- b - a*x(0)
c r(0) <- trans(a)*invtrans(l)*invtrans(u)*inv(u)*inv(l)*r(0) c p(0) <- r(0) c i <- 0
c while r(i) > tol do
c s <- inv(u)*inv(l)*a*p(i)
c a(i) <- trans(r(i))*r(i)/(trans(s)*s) c x(i+1) <- x(i) + a(i)*p(i)
c r(i+1) <- r(i) - a(i)*trans(a)*invtrans(l)*invtrans(u)*s c b(i) <- trans(r(i+1))*r(i+1)/(trans(r(i))*r(i)) c p(i+1) <- r(i+1) + b(i)*p(i) c i <- i + 1
c end c
real ai,bi,rowold,rownew,xnrm,anrm real sdot
real rtxmax,rtxmin
common /gear14/ rtxmax,rtxmin data rtxmin/1.0e-16/ info = 0
anrm = abs(a(isamax(ia(n+1)-1,a,1))) c
c form incomplete factors l and u c
if( job .ne. 0 ) go to 5
call scopy(ia(n+1)-1,a,1,af,1) call lu(af,n,ia,ja) 5 continue c
c r(0) <- b - a*x(0) c
call cgres(a,n,ia,ja,x,b,r) c
c r(0) <- trans(a)*invtrans(l)*invtrans(u)*inv(u)*inv(l)*r(0) c
c inv(u)*inv(l)*r(0) c
call scopy(n,r,1,s,1)
call ssol(af,n,ia,ja,s,1) call ssol(af,n,ia,ja,s,2) c
c invtrans(l)*invtrans(u)*above c
call ssol(af,n,ia,ja,s,-2) call ssol(af,n,ia,ja,s,-1) c
c trans(a)*above c
call mmult(a,n,ia,ja,r,s,-1) c
c p(0) <- r(0) c
call scopy(n,r,1,p,1) rowold = sdot(n,r,1,r,1) i = 0 c
c while r(i) > tol do c
ai = 1.0d0 10 continue
xnrm = abs(x(isamax(n,x,1)))
cc write(6,*) ' iter residual xnrm ai',i,sqrt(rowold)/(anrm*xnrm), cc $ xnrm,ai
if( sqrt(rowold) .le. tol*(anrm*xnrm) ) go to 12 cc if (rowold .le. rtxmin) go to 25 13 continue c
c s <- inv(u)*inv(l)*a*p(i) c
call mmult(a,n,ia,ja,s,p,1) call ssol(af,n,ia,ja,s,1) call ssol(af,n,ia,ja,s,2) c
c a(i) <- trans(r(i))*r(i)/(trans(s)*s) c
ai = rowold/sdot(n,s,1,s,1) c
c x(i+1) <- x(i) + a(i)*p(i) c
call saxpy(n,ai,p,1,x,1) c c
c r(i+1) <- r(i) - a(i)*trans(a)*invtrans(l)*invtrans(u)*s c
call ssol(af,n,ia,ja,s,-2) call ssol(af,n,ia,ja,s,-1) call mmult(a,n,ia,ja,b,s,-1) call saxpy(n,-ai,b,1,r,1) c
c b(i) <- trans(r(i+1))*r(i+1)/(trans(r(i))*r(i)) c
rownew = sdot(n,r,1,r,1) bi = rownew/rowold c
c p(i+1) <- r(i+1) + b(i)*p(i) c
call saxpy2(n,bi,p,1,r,1) rowold = rownew c
c i <- i + 1
c
i = i + 1
if( i .gt. maxits ) go to 20 go to 10 12 continue 15 continue its = i return 20 continue info = -999 its = maxits return c 25 continue c info = 1 c its = i c return end
subroutine axpy(a,n,ia,ja,i,js,je,t,y) c
integer n,ia(1),ja(1),i,js,je real a(1),y(1) real t c
c this routine computes an axpy for a row of a sparse matrix c with a vector.
c an axpy is a multiple of a row of the matrix is added to the c vector y. c
is = ia(i)
ie = ia(i+1) - 1 do 10 ir = is,ie j = ja(ir)
if( j .gt. je ) go to 20 if( j .lt. js ) go to 10 y(j) = y(j) + t*a(ir) 10 continue 20 continue return end
subroutine saxpy2(n,da,dx,incx,dy,incy) c
c constant times a vector plus a vector.
c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78.