不完全cholesky分解实现

t = dot(a,n,ia,ja,i,i+1,n,b)

if( .not. locate(a,n,ia,ja,i,i,k) ) go to 30 b(i) = (b(i) - t)/a(k) 20 continue c

c job = -2 solve trans(u)*x = b c

25 continue

if( job .ne. -2 ) go to 35 c

c solve trans(u)*y=b c

do 21 i = 1,n

if( .not. locate(a,n,ia,ja,i,i,k) ) go to 30 b(i) = b(i)/a(k)

call axpy(a,n,ia,ja,i,i+1,n,-b(i),b) 21 continue c

c solve trans(l)*x=y c

35 continue

if( job .ne. -1 ) go to 45 do 22 ib = 2,n i = n - ib + 2

call axpy(a,n,ia,ja,i,1,i-1,-b(i),b) 22 continue 45 continue return 30 continue

cc write(6,*)' error no diagonal element: from solve' return end

subroutine saxpy(n,sa,sx,incx,sy,incy) c

c constant times a vector plus a vector.

c uses unrolled loop for increments equal to one. c jack dongarra, linpack, 3/11/78. c

real sx(1),sy(1),sa

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

if(n.le.0)return

if (sa .eq. 0.0) 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

sy(iy) = sy(iy) + sa*sx(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

sy(i) = sy(i) + sa*sx(i) 30 continue

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

do 50 i = mp1,n,4

sy(i) = sy(i) + sa*sx(i)

sy(i + 1) = sy(i + 1) + sa*sx(i + 1) sy(i + 2) = sy(i + 2) + sa*sx(i + 2) sy(i + 3) = sy(i + 3) + sa*sx(i + 3) 50 continue return end

subroutine scopy(n,sx,incx,sy,incy) c

c copies a vector, x, to a vector, y.

c uses unrolled loops for increments equal to 1. c jack dongarra, linpack, 3/11/78. c

real sx(1),sy(1)

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

if(n.le.0)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 sy(iy) = sx(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,7)

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

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

do 50 i = mp1,n,7 sy(i) = sx(i)

sy(i + 1) = sx(i + 1) sy(i + 2) = sx(i + 2) sy(i + 3) = sx(i + 3) sy(i + 4) = sx(i + 4) sy(i + 5) = sx(i + 5) sy(i + 6) = sx(i + 6) 50 continue return end

real function sdot(n,sx,incx,sy,incy) c

c forms the dot product of two vectors.

c uses unrolled loops for increments equal to one.

c jack dongarra, linpack, 3/11/78. c

real sx(1),sy(1),stemp

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

stemp = 0.0e0 sdot = 0.0e0 if(n.le.0)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

stemp = stemp + sx(ix)*sy(iy) ix = ix + incx iy = iy + incy 10 continue sdot = stemp return c

c code for both increments equal to 1 c c

c clean-up loop c

20 m = mod(n,5)

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

stemp = stemp + sx(i)*sy(i) 30 continue

if( n .lt. 5 ) go to 60 40 mp1 = m + 1

do 50 i = mp1,n,5

stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) +

* sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4) 50 continue 60 sdot = stemp return end