c
c takes the sum of the absolute values.
c jack dongarra, 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 dzasum(n,zx,incx)
double complex zx(*)
double precision stemp,abs
integer i,incx,ix,n
c
dzasum = 0.0d0
stemp = 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
ix = 1
do 10 i = 1,n
stemp = stemp + abs(zx(ix))
ix = ix + incx
10 continue
dzasum = stemp
return
c
c code for increment equal to 1
c
20 do 30 i = 1,n
stemp = stemp + abs(zx(i))
30 continue
dzasum = stemp
return
end
*
* DZNRM2 returns the euclidean norm of a vector via the function
* name, so that
*
* DZNRM2 := sqrt( conjg( x' )*x )
*
*
*
* -- This version written on 25-October-1982.
* Modified on 14-October-1993 to inline the call to ZLASSQ.
* Sven Hammarling, Nag Ltd.
*
*
DOUBLE PRECISION FUNCTION DZNRM2( N, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, N
* .. Array Arguments ..
COMPLEX*16 X( * )
* ..
* .. Parameters ..
DOUBLE PRECISION ONE , ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* .. Local Scalars ..
INTEGER IX
DOUBLE PRECISION NORM, SCALE, SSQ, TEMP
* .. Intrinsic Functions ..
INTRINSIC ABS, DIMAG, DBLE, SQRT
* ..
* .. Executable Statements ..
IF( N.LT.1 .OR. INCX.LT.1 )THEN
NORM = ZERO
ELSE
SCALE = ZERO
SSQ = ONE
* The following loop is equivalent to this call to the LAPACK
* auxiliary routine:
* CALL ZLASSQ( N, X, INCX, SCALE, SSQ )
*
DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX
IF( DBLE( X( IX ) ).NE.ZERO )THEN
TEMP = ABS( DBLE( X( IX ) ) )
IF( SCALE.LT.TEMP )THEN
SSQ = ONE + SSQ*( SCALE/TEMP )**2
SCALE = TEMP
ELSE
SSQ = SSQ + ( TEMP/SCALE )**2
END IF
END IF
IF( DIMAG( X( IX ) ).NE.ZERO )THEN
TEMP = ABS( DIMAG( X( IX ) ) )
IF( SCALE.LT.TEMP )THEN
SSQ = ONE + SSQ*( SCALE/TEMP )**2
SCALE = TEMP
ELSE
SSQ = SSQ + ( TEMP/SCALE )**2
END IF
END IF
10 CONTINUE
NORM = SCALE * SQRT( SSQ )
END IF
*
DZNRM2 = NORM
RETURN
END
c
c finds the index of element having max. absolute value.
c jack dongarra, 1/15/85.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
integer function izamax(n,zx,incx)
double complex zx(*)
double precision smax
integer i,incx,ix,n
c
izamax = 0
if( n.lt.1 .or. incx.le.0 )return
izamax = 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
smax = abs(zx(1))
ix = ix + incx
do 10 i = 2,n
if(abs(zx(ix)).le.smax) go to 5
izamax = i
smax = abs(zx(ix))
5 ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 smax = abs(zx(1))
do 30 i = 2,n
if(abs(zx(i)).le.smax) go to 30
izamax = i
smax = abs(zx(i))
30 continue
return
end
c
c constant times a vector plus a vector.
c jack dongarra, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
subroutine zaxpy(n,za,zx,incx,zy,incy)
double complex zx(*),zy(*),za
integer i,incx,incy,ix,iy,n
double precision abs
if(n.le.0)return
if (abs(za) .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
zy(iy) = zy(iy) + za*zx(ix)
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
zy(i) = zy(i) + za*zx(i)
30 continue
return
end
c
c copies a vector, x, to a vector, y.
c jack dongarra, linpack, 4/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
subroutine zcopy(n,zx,incx,zy,incy)
double complex zx(*),zy(*)
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
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
zy(iy) = zx(ix)
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
zy(i) = zx(i)
30 continue
return
end
c
c forms the dot product of a vector.
c jack dongarra, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double complex function zdotc(n,zx,incx,zy,incy)
double complex zx(*),zy(*),ztemp
integer i,incx,incy,ix,iy,n
ztemp = (0.0d0,0.0d0)
zdotc = (0.0d0,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
ztemp = ztemp + conjg(zx(ix))*zy(iy)
ix = ix + incx
iy = iy + incy
10 continue
zdotc = ztemp
return
c
c code for both increments equal to 1
c
20 do 30 i = 1,n
ztemp = ztemp + conjg(zx(i))*zy(i)
30 continue
zdotc = ztemp
return
end
c
c forms the dot product of two vectors.
c jack dongarra, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double complex function zdotu(n,zx,incx,zy,incy)
double complex zx(*),zy(*),ztemp
integer i,incx,incy,ix,iy,n
ztemp = (0.0d0,0.0d0)
zdotu = (0.0d0,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
ztemp = ztemp + zx(ix)*zy(iy)
ix = ix + incx
iy = iy + incy
10 continue
zdotu = ztemp
return
c
c code for both increments equal to 1
c
20 do 30 i = 1,n
ztemp = ztemp + zx(i)*zy(i)
30 continue
zdotu = ztemp
return
end
c
c scales a vector by a constant.
c jack dongarra, 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 zdscal(n,da,zx,incx)
double complex zx(*)
double precision da
integer i,incx,ix,n
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
ix = 1
do 10 i = 1,n
zx(ix) = dcmplx(da,0.0d0)*zx(ix)
ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 do 30 i = 1,n
zx(i) = dcmplx(da,0.0d0)*zx(i)
30 continue
return
end
subroutine zrotg(ca,cb,c,s)
double complex ca,cb,s
double precision c
double precision norm,scale
double complex alpha
if (abs(ca) .ne. 0.0d0) go to 10
c = 0.0d0
s = (1.0d0,0.0d0)
ca = cb
go to 20
10 continue
scale = abs(ca) + abs(cb)
norm = scale*dsqrt((abs(ca/cmplx(scale,0.0d0)))**2 +
* (abs(cb/cmplx(scale,0.0d0)))**2)
alpha = ca /abs(ca)
c = abs(ca) / norm
s = alpha * conjg(cb) / norm
ca = alpha * norm
20 continue
return
end
c
c scales a vector by a constant.
c jack dongarra, 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 zscal(n,za,zx,incx)
double complex za,zx(*)
integer i,incx,ix,n
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
ix = 1
do 10 i = 1,n
zx(ix) = za*zx(ix)
ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 do 30 i = 1,n
zx(i) = za*zx(i)
30 continue
return
end
c
c interchanges two vectors.
c jack dongarra, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
subroutine zswap (n,zx,incx,zy,incy)
double complex zx(*),zy(*),ztemp
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
ztemp = zx(ix)
zx(ix) = zy(iy)
zy(iy) = ztemp
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
20 do 30 i = 1,n
ztemp = zx(i)
zx(i) = zy(i)
zy(i) = ztemp
30 continue
return
end