c

real dx(1),dy(1),da

integer i,incx,incy,m,mp1,n c

if(n.le.0)return

if (da .eq. 0.0e0) return

if(incx.eq.1.and.incy.eq.1)go to 20 c

c code for unequal increments or equal increments c not equal to 1 c

ix = 1 iy = 1

if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n

dx(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c

c code for both increments equal to 1 c c

c clean-up loop c

20 m = mod(n,4)

if( m .eq. 0 ) go to 40 do 30 i = 1,m

dx(i) = dy(i) + da*dx(i) 30 continue

if( n .lt. 4 ) return 40 mp1 = m + 1

do 50 i = mp1,n,4

dx(i) = dy(i) + da*dx(i)

dx(i + 1) = dy(i + 1) + da*dx(i + 1) dx(i + 2) = dy(i + 2) + da*dx(i + 2) dx(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end

real function dot(a,n,ia,ja,i,js,je,b) c

integer n,ia(1),ja(1),i,js,je real a(1),b(1) real t c

c this routine computes an inner product for a row of a sparse matri c with a vector. c

t = 0.0e0 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 t = t + a(ir)*b(j) 10 continue 20 continue dot = t return end

subroutine inits(a,n,ia,ja) c

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

c this routine initializes storage for the sparse c matrix arrays.

c note: the matrix is initialized to have zeroes c on the diagonal. c

do 10 i = 1,n a(i) = 0.0e0 ja(i) = i ia(i) = i 10 continue

ia(n+1) = n+1 return end

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

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

c this routine rearranges the elements in arrays a and ja c and updates array ia for the new element. c

cc write(6,1000)i,j,ia(i),ia(i+1),t

1000 format(' *** from insert i,j,ia(i),ia(i+1),t ',4i5,1x,e15.8) l1 = k

l2 = ia(n+1) - 1 do 10 lb = l1,l2 l = l2 - lb + l1 a(l+1) = a(l) ja(l+1) = ja(l) 10 continue a(k) = t ja(k) = j ip1 = i + 1 np1 = n + 1

do 20 l = ip1,np1 ia(l) = ia(l) + 1 20 continue return end

logical function locate(a,n,ia,ja,i,j,k) c

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

c this routine will locate the i,j-th element in the c sparse matrix structure. c

is = ia(i)

ie = ia(i+1) - 1 c

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

do 10 l = is,ie

if( j .gt. ja(l) ) go to 10 if( j .ne. ja(l) ) go to 5 locate = .true.

k = l ！l为对应的i,j元素在向量a中的位置 go to 20 5 continue

if( j .ge. ja(l) ) go to 10 locate = .false. k = 0

go to 20 10 continue c

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

locate = .false. k = 0 20 continue return end

subroutine lu(a,n,ia,ja) c

logical locate

integer n,ia(1),ja(1) integer ikp1,kp1 real a(1) c

c this subroutine does incomplete gaussian elimenation c on a sparse matrix. the matrix is stored in a sparse c data structure.

c note: no pivoting is done

c and the factorization is incomplete,

c i.e., only where there exists a storage location c will the operation take place. c

do 40 k = 1,n

if( .not. locate(a,n,ia,ja,k,k,l2) ) go to 50 if( a(l2) .eq. 0.0e0 ) go to 50 kp1 = k + 1 do 30 i = kp1,n

if( .not. locate(a,n,ia,ja,i,k,ik) ) go to 25 is = ik

ie = ia(i+1) - 1 kj = l2

ke = ia(k+1) - 1 a(ik) = a(ik)/a(kj)

if( kj .eq. ke ) go to 30 kj = kj + 1 ikp1 = ik + 1

do 20 j = ikp1,ie 10 continue

if( kj .gt. ke ) go to 30

if( ja(kj) .ge. ja(j) ) go to 15 kj = kj + 1