不完全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.