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