c
c takes the sum of the absolute values.
c jack dongarra, linpack, 3/11/78.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision function dasum(n,dx,incx)
double precision dx(*),dtemp
integer i,incx,m,mp1,n,nincx
c
dasum = 0.0d0
dtemp = 0.0d0
if( n.le.0 .or. incx.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
dtemp = dtemp + dabs(dx(i))
10 continue
dasum = dtemp
return
c
c code for increment equal to 1
c
c
c clean-up loop
c
20 m = mod(n,6)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dabs(dx(i))
30 continue
if( n .lt. 6 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,6
dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2))
* + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5))
50 continue
60 dasum = dtemp
return
end
c
c constant times a vector plus a vector.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
subroutine daxpy(n,da,dx,incx,dy,incy)
double precision dx(*),dy(*),da
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if (da .eq. 0.0d0) 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
dy(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
dy(i) = dy(i) + da*dx(i)
30 continue
if( n .lt. 4 ) return
40 mp1 = m + 1
do 50 i = mp1,n,4
dy(i) = dy(i) + da*dx(i)
dy(i + 1) = dy(i + 1) + da*dx(i + 1)
dy(i + 2) = dy(i + 2) + da*dx(i + 2)
dy(i + 3) = dy(i + 3) + da*dx(i + 3)
50 continue
return
end
c
c copies a vector, x, to a vector, y.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
subroutine dcopy(n,dx,incx,dy,incy)
double precision dx(*),dy(*)
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
dy(iy) = 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,7)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dy(i) = dx(i)
30 continue
if( n .lt. 7 ) return
40 mp1 = m + 1
do 50 i = mp1,n,7
dy(i) = dx(i)
dy(i + 1) = dx(i + 1)
dy(i + 2) = dx(i + 2)
dy(i + 3) = dx(i + 3)
dy(i + 4) = dx(i + 4)
dy(i + 5) = dx(i + 5)
dy(i + 6) = dx(i + 6)
50 continue
return
end
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 modified 12/3/93, array(1) declarations changed to array(*)
c
double precision function ddot(n,dx,incx,dy,incy)
double precision dx(*),dy(*),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
c
ddot = 0.0d0
dtemp = 0.0d0
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
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
10 continue
ddot = dtemp
return
c
c code for both increments equal to 1
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dx(i)*dy(i)
30 continue
if( n .lt. 5 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,5
dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
* dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
50 continue
60 ddot = dtemp
return
end
c
c smach computes machine parameters of floating point
c arithmetic for use in testing only. not required by
c linpack proper.
c
c if trouble with automatic computation of these quantities,
c they can be set by direct assignment statements.
c assume the computer has
c
c b = base of arithmetic
c t = number of base b digits
c l = smallest possible exponent
c u = largest possible exponent
c
c then
c
c eps = b**(1-t)
c tiny = 100.0*b**(-l+t)
c huge = 0.01*b**(u-t)
c
c dmach same as smach except t, l, u apply to
c double precision.
c
c cmach same as smach except if complex division
c is done by
c
c 1/(x+i*y) = (x-i*y)/(x**2+y**2)
c
c then
c
c tiny = sqrt(tiny)
c huge = sqrt(huge)
c
c
c job is 1, 2 or 3 for epsilon, tiny and huge, respectively.
c
double precision function dmach(job)
integer job
double precision eps,tiny,huge,s
c
eps = 1.0d0
10 eps = eps/2.0d0
s = 1.0d0 + eps
if (s .gt. 1.0d0) go to 10
eps = 2.0d0*eps
c
s = 1.0d0
20 tiny = s
s = s/16.0d0
if (s*1.0 .ne. 0.0d0) go to 20
tiny = (tiny/eps)*100.0
huge = 1.0d0/tiny
c
if (job .eq. 1) dmach = eps
if (job .eq. 2) dmach = tiny
if (job .eq. 3) dmach = huge
return
end
*
* DNRM2 returns the euclidean norm of a vector via the function
* name, so that
*
* DNRM2 := sqrt( x'*x )
*
*
*
* -- This version written on 25-October-1982.
* Modified on 14-October-1993 to inline the call to DLASSQ.
* Sven Hammarling, Nag Ltd.
*
*
DOUBLE PRECISION FUNCTION DNRM2 ( N, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, N
* .. Array Arguments ..
DOUBLE PRECISION X( * )
* .. Parameters ..
DOUBLE PRECISION ONE , ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* .. Local Scalars ..
INTEGER IX
DOUBLE PRECISION ABSXI, NORM, SCALE, SSQ
* .. Intrinsic Functions ..
INTRINSIC ABS, SQRT
* ..
* .. Executable Statements ..
IF( N.LT.1 .OR. INCX.LT.1 )THEN
NORM = ZERO
ELSE IF( N.EQ.1 )THEN
NORM = ABS( X( 1 ) )
ELSE
SCALE = ZERO
SSQ = ONE
* The following loop is equivalent to this call to the LAPACK
* auxiliary routine:
* CALL DLASSQ( N, X, INCX, SCALE, SSQ )
*
DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX
IF( X( IX ).NE.ZERO )THEN
ABSXI = ABS( X( IX ) )
IF( SCALE.LT.ABSXI )THEN
SSQ = ONE + SSQ*( SCALE/ABSXI )**2
SCALE = ABSXI
ELSE
SSQ = SSQ + ( ABSXI/SCALE )**2
END IF
END IF
10 CONTINUE
NORM = SCALE * SQRT( SSQ )
END IF
*
DNRM2 = NORM
RETURN
*
* End of DNRM2.
*
END
c
c applies a plane rotation.
c jack dongarra, linpack, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
subroutine drot (n,dx,incx,dy,incy,c,s)
double precision dx(*),dy(*),dtemp,c,s
integer i,incx,incy,ix,iy,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
dtemp = c*dx(ix) + s*dy(iy)
dy(iy) = c*dy(iy) - s*dx(ix)
dx(ix) = dtemp
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
20 do 30 i = 1,n
dtemp = c*dx(i) + s*dy(i)
dy(i) = c*dy(i) - s*dx(i)
dx(i) = dtemp
30 continue
return
end
c
c construct givens plane rotation.
c jack dongarra, linpack, 3/11/78.
c
subroutine drotg(da,db,c,s)
double precision da,db,c,s,roe,scale,r,z
c
roe = db
if( dabs(da) .gt. dabs(db) ) roe = da
scale = dabs(da) + dabs(db)
if( scale .ne. 0.0d0 ) go to 10
c = 1.0d0
s = 0.0d0
r = 0.0d0
z = 0.0d0
go to 20
10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2)
r = dsign(1.0d0,roe)*r
c = da/r
s = db/r
z = 1.0d0
if( dabs(da) .gt. dabs(db) ) z = s
if( dabs(db) .ge. dabs(da) .and. c .ne. 0.0d0 ) z = 1.0d0/c
20 da = r
db = z
return
end
c
c scales a vector by a constant.
c uses unrolled loops for increment equal to one.
c jack dongarra, linpack, 3/11/78.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
subroutine dscal(n,da,dx,incx)
double precision da,dx(*)
integer i,incx,m,mp1,n,nincx
c
if( n.le.0 .or. incx.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
dx(i) = da*dx(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
dx(i) = da*dx(i)
30 continue
if( n .lt. 5 ) return
40 mp1 = m + 1
do 50 i = mp1,n,5
dx(i) = da*dx(i)
dx(i + 1) = da*dx(i + 1)
dx(i + 2) = da*dx(i + 2)
dx(i + 3) = da*dx(i + 3)
dx(i + 4) = da*dx(i + 4)
50 continue
return
end
c
c interchanges two vectors.
c uses unrolled loops for increments equal one.
c jack dongarra, linpack, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
subroutine dswap (n,dx,incx,dy,incy)
double precision dx(*),dy(*),dtemp
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
dtemp = dx(ix)
dx(ix) = dy(iy)
dy(iy) = dtemp
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
c clean-up loop
c
20 m = mod(n,3)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dx(i)
dx(i) = dy(i)
dy(i) = dtemp
30 continue
if( n .lt. 3 ) return
40 mp1 = m + 1
do 50 i = mp1,n,3
dtemp = dx(i)
dx(i) = dy(i)
dy(i) = dtemp
dtemp = dx(i + 1)
dx(i + 1) = dy(i + 1)
dy(i + 1) = dtemp
dtemp = dx(i + 2)
dx(i + 2) = dy(i + 2)
dy(i + 2) = dtemp
50 continue
return
end
c
c finds the index of element having max. absolute value.
c jack dongarra, linpack, 3/11/78.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
integer function idamax(n,dx,incx)
double precision dx(*),dmax
integer i,incx,ix,n
c
idamax = 0
if( n.lt.1 .or. incx.le.0 ) return
idamax = 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
dmax = dabs(dx(1))
ix = ix + incx
do 10 i = 2,n
if(dabs(dx(ix)).le.dmax) go to 5
idamax = i
dmax = dabs(dx(ix))
5 ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 dmax = dabs(dx(1))
do 30 i = 2,n
if(dabs(dx(i)).le.dmax) go to 30
idamax = i
dmax = dabs(dx(i))
30 continue
return
end