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