不完全cholesky分解实现

go to 10 15 continue

if( ja(kj) .gt. ja(j) ) go to 20 a(j) = a(j) - a(ik)*a(kj) 20 continue 25 continue 30 continue 40 continue return 50 continue

cc write(6,*)' value of k and a(k,k)',k,a(l2) cc write(6,*)' error zero on diagonal'

cc write(6,*)' matrix probably specified incorrectly or'

cc write(6,*)' incomplete factorization produces a singular matrix' return end

subroutine mmult(a,n,ia,ja,b,x,job) c

integer n,ia(1),ja(1),job real a(1),b(1),x(1) c

c this routine performs matrix vector multiple c for a sparse matrix structure c

c job has the following input meanings: c job = 1 matrix vector multiple

c = -1 matrix transpose vector multiple c = 2 unit lower matrix vector multiple

c = -2 unit lower matrix transpose vector multiple c = 3 upper matrix vector multiple

c = -3 upper matrix transpose vector multiple real dot c

c a*x c

if( job .ne. 1 ) go to 15 do 10 i = 1,n

b(i) = dot(a,n,ia,ja,i,1,n,x) 10 continue c

c trans(a)*x c

15 continue

if( job .ne. -1 ) go to 35

do 20 i = 1,n b(i) = 0.0e0 20 continue

do 30 i = 1,n

call axpy(a,n,ia,ja,i,1,n,x(i),b) 30 continue c

c l*x when l is unit lower c

35 continue

if( job .ne. 2 ) go to 45 do 40 i = 1,n

b(i) = x(i) + dot(a,n,ia,ja,i,1,i-1,x) 40 continue c

c trans(l)*x when l is unit lower c

45 continue

if( job .ne. -2 ) go to 55 do 49 i = 1,n b(i) = x(i) 49 continue

do 50 i = 1,n

call axpy(a,n,ia,ja,i,i,n,x(i),b) 50 continue c

c u*x c

55 continue

if( job .ne. 3 ) go to 65 do 60 i = 1,n

b(i) = dot(a,n,ia,ja,i,i,n,x) 60 continue c

c trans(u)*x c

65 continue

if( job .ne. -3 ) go to 85 do 70 i = 1,n b(i) = 0.0e0 70 continue

do 80 i = 1,n

call axpy(a,n,ia,ja,i,1,i,x(i),b) 80 continue

85 continue return end

subroutine put(t,a,n,ia,ja,i,j) c

integer n,ia(1),ja(1),i,j real t,a(1) c

c this routine will insert an element into the sparse matrix c structure. c

is = ia(i)

ie = ia(i+1) - 1

cc write(6,100)i,j,ia(i),ia(i+1)

100 format(' *** from put i,j,ia(i),ia(i+1) ',4i7) c

c row i is from is to ie in array a. c

do 10 k = is,ie

if( j .gt. ja(k) ) go to 10 if( j .ne. ja(k) ) go to 5 a(k) = t go to 20 5 continue

if( j .ge. ja(k) ) go to 12

call insert(t,a,n,ia,ja,i,j,k) go to 20 12 continue 10 continue c

c get here if should be at the end of a row. c

k = ie + 1

call insert(t,a,n,ia,ja,i,j,k) go to 30 20 continue 30 continue return end

subroutine cgres(a,n,ia,ja,x,b,r) c

integer n,ia(1),ja(1) real a(1),x(1),b(1),r(1) c

c this routine computes a residual for a*x=b where c a is in a sparse structure c

real dot

do 10 i = 1,n

r(i) = b(i) - dot(a,n,ia,ja,i,1,n,x) 10 continue return end

subroutine ssol(a,n,ia,ja,b,job) c

integer n,ia(1),ja(1),job real a(1),b(1) c

c this routine solves a system of equations based on a sparse c matrix date structure.

c the array b contains the right hand side on input c and on output has the solution c

c job has the value 1 if l*x = b is to be solved.

c job has the value -1 if trans(l)*x = b is to be solved. c job has the value 2 if u*x = b is to be solved.

c job has the value -2 if trans(u)*x = b is to be solved. c

logical locate real t real dot c

c job = 1 solve l*x = b c

if( job .ne. 1 ) go to 15 c

c solve l*y=b c

do 10 i = 2,n

b(i) = b(i) - dot(a,n,ia,ja,i,1,i-1,b) 10 continue 15 continue

if( job .ne. 2 ) go to 25 c

c solve u*x=y c

do 20 ib = 1,n i = n - ib + 1