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
integer function isamax(n,sx,incx)
c
c finds the index of element having max. absolute value.
c jack dongarra, linpack, 3/11/78.
c modified to correct problem with negative increments, 9/29/88.
c
real sx(1),smax
integer i,incx,ix,n
c
isamax = 0
if( n .lt. 1 ) return
isamax = 1
if(n.eq.1)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
smax = abs(sx(ix))
ix = ix + incx
do 10 i = 2,n
if(abs(sx(ix)).le.smax) go to 5
isamax = i
smax = abs(sx(ix))
5 ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 smax = abs(sx(1))
do 30 i = 2,n
if(abs(sx(i)).le.smax) go to 30
isamax = i
smax = abs(sx(i))
30 continue
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 snrm2 ( n, sx, incx)
integer next
real sx(1), cutlo, cuthi, hitest, sum, xmax, zero, one
data zero, one /0.0e0, 1.0e0/
c
c euclidean norm of the n-vector stored in sx() with storage
c increment incx .
c if n .le. 0 return with result = 0.
c if n .ge. 1 then incx must be .ge. 1
c
c c.l.lawson, 1978 jan 08
c
c four phase method using two built-in constants that are
c hopefully applicable to all machines.
c cutlo = maximum of sqrt(u/eps) over all known machines.
c cuthi = minimum of sqrt(v) over all known machines.
c where
c eps = smallest no. such that eps + 1. .gt. 1.
c u = smallest positive no. (underflow limit)
c v = largest no. (overflow limit)
c
c brief outline of algorithm..
c
c phase 1 scans zero components.
c move to phase 2 when a component is nonzero and .le. cutlo
c move to phase 3 when a component is .gt. cutlo
c move to phase 4 when a component is .ge. cuthi/m
c where m = n for x() real and m = 2*n for complex.
c
c values for cutlo and cuthi..
c from the environmental parameters listed in the imsl converter
c document the limiting values are as follows..
c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are
c univac and dec at 2**(-103)
c thus cutlo = 2**(-51) = 4.44089e-16
c cuthi, s.p. v = 2**127 for univac, honeywell, and dec.
c thus cuthi = 2**(63.5) = 1.30438e19
c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec.
c thus cutlo = 2**(-33.5) = 8.23181d-11
c cuthi, d.p. same as s.p. cuthi = 1.30438d19
c data cutlo, cuthi / 8.232d-11, 1.304d19 /
c data cutlo, cuthi / 4.441e-16, 1.304e19 /
data cutlo, cuthi / 4.441e-16, 1.304e19 /
c
if(n .gt. 0) go to 10
snrm2 = zero
go to 300
c
10 assign 30 to next
sum = zero
nn = n * incx
c begin main loop
i = 1
20 go to next,(30, 50, 70, 110)
30 if( abs(sx(i)) .gt. cutlo) go to 85
assign 50 to next
xmax = zero
c
c phase 1. sum is zero
c
50 if( sx(i) .eq. zero) go to 200
if( abs(sx(i)) .gt. cutlo) go to 85
c
c prepare for phase 2.
assign 70 to next
go to 105
c
c prepare for phase 4.
c
100 i = j
assign 110 to next
sum = (sum / sx(i)) / sx(i)
105 xmax = abs(sx(i))
go to 115
c
c phase 2. sum is small.
c scale to avoid destructive underflow.
c
70 if( abs(sx(i)) .gt. cutlo ) go to 75
c
c common code for phases 2 and 4.
c in phase 4 sum is large. scale to avoid overflow.
c
110 if( abs(sx(i)) .le. xmax ) go to 115
sum = one + sum * (xmax / sx(i))**2
xmax = abs(sx(i))
go to 200
c
115 sum = sum + (sx(i)/xmax)**2
go to 200
c
c
c prepare for phase 3.
c
75 sum = (sum * xmax) * xmax
c
c
c for real or d.p. set hitest = cuthi/n
c for complex set hitest = cuthi/(2*n)
c
85 hitest = cuthi/float( n )
c
c phase 3. sum is mid-range. no scaling.
c
do 95 j =i,nn,incx
if(abs(sx(j)) .ge. hitest) go to 100
95 sum = sum + sx(j)**2
snrm2 = sqrt( sum )
go to 300
c
200 continue
i = i + incx
if ( i .le. nn ) go to 20
c
c end of main loop.
c
c compute square root and adjust for scaling.
c
snrm2 = xmax * sqrt(sum)
300 continue
return
end
subroutine sscal(n,sa,sx,incx)
c
c scales a vector by a constant.
c uses unrolled loops for increment equal to 1.
c jack dongarra, linpack, 3/11/78.
c
real sa,sx(1)
integer i,incx,m,mp1,n,nincx
c
if(n.le.0)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
nincx = n*incx
do 10 i = 1,nincx,incx
sx(i) = sa*sx(i)
10 continue
return
c
c code for increment 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
sx(i) = sa*sx(i)
30 continue
if( n .lt. 5 ) return
40 mp1 = m + 1
do 50 i = mp1,n,5
sx(i) = sa*sx(i)
sx(i + 1) = sa*sx(i + 1)
sx(i + 2) = sa*sx(i + 2)
sx(i + 3) = sa*sx(i + 3)
sx(i + 4) = sa*sx(i + 4)
50 continue
return
end
subroutine sswap (n,sx,incx,sy,incy)
c
c interchanges two vectors.
c uses unrolled loops for increments equal to 1.
c jack dongarra, linpack, 3/11/78.
c
real sx(1),sy(1),stemp
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 not equal
c 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 = sx(ix)
sx(ix) = sy(iy)
sy(iy) = stemp
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,3)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
stemp = sx(i)
sx(i) = sy(i)
sy(i) = stemp
30 continue
if( n .lt. 3 ) return
40 mp1 = m + 1
do 50 i = mp1,n,3
stemp = sx(i)
sx(i) = sy(i)
sy(i) = stemp
stemp = sx(i + 1)
sx(i + 1) = sy(i + 1)
sy(i + 1) = stemp
stemp = sx(i + 2)
sx(i + 2) = sy(i + 2)
sy(i + 2) = stemp
50 continue
return
end