double precision function dcabs1(z)
double complex z,zz
double precision t(2)
equivalence (zz,t(1))
zz = z
dcabs1 = dabs(t(1)) + dabs(t(2))
return
end
double precision function dzasum(n,zx,incx)
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 complex zx(*)
double precision stemp,dcabs1
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 + dcabs1(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 + dcabs1(zx(i))
30 continue
dzasum = stemp
return
end
DOUBLE PRECISION FUNCTION DZNRM2( N, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, N
* .. Array Arguments ..
COMPLEX*16 X( * )
* ..
*
* 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.
*
*
* .. 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 of DZNRM2.
*
END
integer function izamax(n,zx,incx)
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
double complex zx(*)
double precision smax
integer i,incx,ix,n
double precision dcabs1
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 = dcabs1(zx(1))
ix = ix + incx
do 10 i = 2,n
if(dcabs1(zx(ix)).le.smax) go to 5
izamax = i
smax = dcabs1(zx(ix))
5 ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 smax = dcabs1(zx(1))
do 30 i = 2,n
if(dcabs1(zx(i)).le.smax) go to 30
izamax = i
smax = dcabs1(zx(i))
30 continue
return
end
subroutine zaxpy(n,za,zx,incx,zy,incy)
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
double complex zx(*),zy(*),za
integer i,incx,incy,ix,iy,n
double precision dcabs1
if(n.le.0)return
if (dcabs1(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
subroutine zcopy(n,zx,incx,zy,incy)
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
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
double complex function zdotc(n,zx,incx,zy,incy)
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 zx(*),zy(*),ztemp
integer i,incx,incy,ix,iy,n
intrinsic dconjg
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 + dconjg(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 + dconjg(zx(i))*zy(i)
30 continue
zdotc = ztemp
return
end
double complex function zdotu(n,zx,incx,zy,incy)
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 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
subroutine zdscal(n,da,zx,incx)
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
intrinsic dcmplx
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 ZGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
$ BETA, Y, INCY )
* .. Scalar Arguments ..
COMPLEX*16 ALPHA, BETA
INTEGER INCX, INCY, LDA, M, N
CHARACTER*1 TRANS
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * ), Y( * )
* ..
*
* Purpose
* =======
*
* ZGEMV performs one of the matrix-vector operations
*
* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, or
*
* y := alpha*conjg( A' )*x + beta*y,
*
* where alpha and beta are scalars, x and y are vectors and A is an
* m by n matrix.
*
* Parameters
* ==========
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
*
* TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
*
* TRANS = 'C' or 'c' y := alpha*conjg( A' )*x + beta*y.
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of the matrix A.
* M must be at least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry, the leading m by n part of the array A must
* contain the matrix of coefficients.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, m ).
* Unchanged on exit.
*
* X - COMPLEX*16 array of DIMENSION at least
* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
* and at least
* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
* Before entry, the incremented array X must contain the
* vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* BETA - COMPLEX*16 .
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then Y need not be set on input.
* Unchanged on exit.
*
* Y - COMPLEX*16 array of DIMENSION at least
* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
* and at least
* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
* Before entry with BETA non-zero, the incremented array Y
* must contain the vector y. On exit, Y is overwritten by the
* updated vector y.
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
LOGICAL NOCONJ
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( TRANS, 'N' ).AND.
$ .NOT.LSAME( TRANS, 'T' ).AND.
$ .NOT.LSAME( TRANS, 'C' ) )THEN
INFO = 1
ELSE IF( M.LT.0 )THEN
INFO = 2
ELSE IF( N.LT.0 )THEN
INFO = 3
ELSE IF( LDA.LT.MAX( 1, M ) )THEN
INFO = 6
ELSE IF( INCX.EQ.0 )THEN
INFO = 8
ELSE IF( INCY.EQ.0 )THEN
INFO = 11
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZGEMV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
$ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
*
* Set LENX and LENY, the lengths of the vectors x and y, and set
* up the start points in X and Y.
*
IF( LSAME( TRANS, 'N' ) )THEN
LENX = N
LENY = M
ELSE
LENX = M
LENY = N
END IF
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( LENX - 1 )*INCX
END IF
IF( INCY.GT.0 )THEN
KY = 1
ELSE
KY = 1 - ( LENY - 1 )*INCY
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
* First form y := beta*y.
*
IF( BETA.NE.ONE )THEN
IF( INCY.EQ.1 )THEN
IF( BETA.EQ.ZERO )THEN
DO 10, I = 1, LENY
Y( I ) = ZERO
10 CONTINUE
ELSE
DO 20, I = 1, LENY
Y( I ) = BETA*Y( I )
20 CONTINUE
END IF
ELSE
IY = KY
IF( BETA.EQ.ZERO )THEN
DO 30, I = 1, LENY
Y( IY ) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40, I = 1, LENY
Y( IY ) = BETA*Y( IY )
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF( ALPHA.EQ.ZERO )
$ RETURN
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form y := alpha*A*x + y.
*
JX = KX
IF( INCY.EQ.1 )THEN
DO 60, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*X( JX )
DO 50, I = 1, M
Y( I ) = Y( I ) + TEMP*A( I, J )
50 CONTINUE
END IF
JX = JX + INCX
60 CONTINUE
ELSE
DO 80, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*X( JX )
IY = KY
DO 70, I = 1, M
Y( IY ) = Y( IY ) + TEMP*A( I, J )
IY = IY + INCY
70 CONTINUE
END IF
JX = JX + INCX
80 CONTINUE
END IF
ELSE
*
* Form y := alpha*A'*x + y or y := alpha*conjg( A' )*x + y.
*
JY = KY
IF( INCX.EQ.1 )THEN
DO 110, J = 1, N
TEMP = ZERO
IF( NOCONJ )THEN
DO 90, I = 1, M
TEMP = TEMP + A( I, J )*X( I )
90 CONTINUE
ELSE
DO 100, I = 1, M
TEMP = TEMP + DCONJG( A( I, J ) )*X( I )
100 CONTINUE
END IF
Y( JY ) = Y( JY ) + ALPHA*TEMP
JY = JY + INCY
110 CONTINUE
ELSE
DO 140, J = 1, N
TEMP = ZERO
IX = KX
IF( NOCONJ )THEN
DO 120, I = 1, M
TEMP = TEMP + A( I, J )*X( IX )
IX = IX + INCX
120 CONTINUE
ELSE
DO 130, I = 1, M
TEMP = TEMP + DCONJG( A( I, J ) )*X( IX )
IX = IX + INCX
130 CONTINUE
END IF
Y( JY ) = Y( JY ) + ALPHA*TEMP
JY = JY + INCY
140 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZGEMV .
*
END
SUBROUTINE ZGERC ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
* .. Scalar Arguments ..
COMPLEX*16 ALPHA
INTEGER INCX, INCY, LDA, M, N
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * ), Y( * )
* ..
*
* Purpose
* =======
*
* ZGERC performs the rank 1 operation
*
* A := alpha*x*conjg( y' ) + A,
*
* where alpha is a scalar, x is an m element vector, y is an n element
* vector and A is an m by n matrix.
*
* Parameters
* ==========
*
* M - INTEGER.
* On entry, M specifies the number of rows of the matrix A.
* M must be at least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( m - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the m
* element vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* Y - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCY ) ).
* Before entry, the incremented array Y must contain the n
* element vector y.
* Unchanged on exit.
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry, the leading m by n part of the array A must
* contain the matrix of coefficients. On exit, A is
* overwritten by the updated matrix.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, m ).
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, J, JY, KX
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( M.LT.0 )THEN
INFO = 1
ELSE IF( N.LT.0 )THEN
INFO = 2
ELSE IF( INCX.EQ.0 )THEN
INFO = 5
ELSE IF( INCY.EQ.0 )THEN
INFO = 7
ELSE IF( LDA.LT.MAX( 1, M ) )THEN
INFO = 9
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZGERC ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
$ RETURN
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF( INCY.GT.0 )THEN
JY = 1
ELSE
JY = 1 - ( N - 1 )*INCY
END IF
IF( INCX.EQ.1 )THEN
DO 20, J = 1, N
IF( Y( JY ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( Y( JY ) )
DO 10, I = 1, M
A( I, J ) = A( I, J ) + X( I )*TEMP
10 CONTINUE
END IF
JY = JY + INCY
20 CONTINUE
ELSE
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( M - 1 )*INCX
END IF
DO 40, J = 1, N
IF( Y( JY ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( Y( JY ) )
IX = KX
DO 30, I = 1, M
A( I, J ) = A( I, J ) + X( IX )*TEMP
IX = IX + INCX
30 CONTINUE
END IF
JY = JY + INCY
40 CONTINUE
END IF
*
RETURN
*
* End of ZGERC .
*
END
SUBROUTINE ZHEMV ( UPLO, N, ALPHA, A, LDA, X, INCX,
$ BETA, Y, INCY )
* .. Scalar Arguments ..
COMPLEX*16 ALPHA, BETA
INTEGER INCX, INCY, LDA, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * ), Y( * )
* ..
*
* Purpose
* =======
*
* ZHEMV performs the matrix-vector operation
*
* y := alpha*A*x + beta*y,
*
* where alpha and beta are scalars, x and y are n element vectors and
* A is an n by n hermitian matrix.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the array A is to be referenced as
* follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of A
* is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of A
* is to be referenced.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array A must contain the upper
* triangular part of the hermitian matrix and the strictly
* lower triangular part of A is not referenced.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array A must contain the lower
* triangular part of the hermitian matrix and the strictly
* upper triangular part of A is not referenced.
* Note that the imaginary parts of the diagonal elements need
* not be set and are assumed to be zero.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, n ).
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* BETA - COMPLEX*16 .
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then Y need not be set on input.
* Unchanged on exit.
*
* Y - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCY ) ).
* Before entry, the incremented array Y must contain the n
* element vector y. On exit, Y is overwritten by the updated
* vector y.
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP1, TEMP2
INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX, DBLE
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
$ .NOT.LSAME( UPLO, 'L' ) )THEN
INFO = 1
ELSE IF( N.LT.0 )THEN
INFO = 2
ELSE IF( LDA.LT.MAX( 1, N ) )THEN
INFO = 5
ELSE IF( INCX.EQ.0 )THEN
INFO = 7
ELSE IF( INCY.EQ.0 )THEN
INFO = 10
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZHEMV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* Set up the start points in X and Y.
*
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( N - 1 )*INCX
END IF
IF( INCY.GT.0 )THEN
KY = 1
ELSE
KY = 1 - ( N - 1 )*INCY
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through the triangular part
* of A.
*
* First form y := beta*y.
*
IF( BETA.NE.ONE )THEN
IF( INCY.EQ.1 )THEN
IF( BETA.EQ.ZERO )THEN
DO 10, I = 1, N
Y( I ) = ZERO
10 CONTINUE
ELSE
DO 20, I = 1, N
Y( I ) = BETA*Y( I )
20 CONTINUE
END IF
ELSE
IY = KY
IF( BETA.EQ.ZERO )THEN
DO 30, I = 1, N
Y( IY ) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40, I = 1, N
Y( IY ) = BETA*Y( IY )
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF( ALPHA.EQ.ZERO )
$ RETURN
IF( LSAME( UPLO, 'U' ) )THEN
*
* Form y when A is stored in upper triangle.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 60, J = 1, N
TEMP1 = ALPHA*X( J )
TEMP2 = ZERO
DO 50, I = 1, J - 1
Y( I ) = Y( I ) + TEMP1*A( I, J )
TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( I )
50 CONTINUE
Y( J ) = Y( J ) + TEMP1*DBLE( A( J, J ) ) + ALPHA*TEMP2
60 CONTINUE
ELSE
JX = KX
JY = KY
DO 80, J = 1, N
TEMP1 = ALPHA*X( JX )
TEMP2 = ZERO
IX = KX
IY = KY
DO 70, I = 1, J - 1
Y( IY ) = Y( IY ) + TEMP1*A( I, J )
TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( IX )
IX = IX + INCX
IY = IY + INCY
70 CONTINUE
Y( JY ) = Y( JY ) + TEMP1*DBLE( A( J, J ) ) + ALPHA*TEMP2
JX = JX + INCX
JY = JY + INCY
80 CONTINUE
END IF
ELSE
*
* Form y when A is stored in lower triangle.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 100, J = 1, N
TEMP1 = ALPHA*X( J )
TEMP2 = ZERO
Y( J ) = Y( J ) + TEMP1*DBLE( A( J, J ) )
DO 90, I = J + 1, N
Y( I ) = Y( I ) + TEMP1*A( I, J )
TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( I )
90 CONTINUE
Y( J ) = Y( J ) + ALPHA*TEMP2
100 CONTINUE
ELSE
JX = KX
JY = KY
DO 120, J = 1, N
TEMP1 = ALPHA*X( JX )
TEMP2 = ZERO
Y( JY ) = Y( JY ) + TEMP1*DBLE( A( J, J ) )
IX = JX
IY = JY
DO 110, I = J + 1, N
IX = IX + INCX
IY = IY + INCY
Y( IY ) = Y( IY ) + TEMP1*A( I, J )
TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( IX )
110 CONTINUE
Y( JY ) = Y( JY ) + ALPHA*TEMP2
JX = JX + INCX
JY = JY + INCY
120 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZHEMV .
*
END
SUBROUTINE ZHER2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
* .. Scalar Arguments ..
COMPLEX*16 ALPHA
INTEGER INCX, INCY, LDA, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * ), Y( * )
* ..
*
* Purpose
* =======
*
* ZHER2 performs the hermitian rank 2 operation
*
* A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
*
* where alpha is a scalar, x and y are n element vectors and A is an n
* by n hermitian matrix.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the array A is to be referenced as
* follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of A
* is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of A
* is to be referenced.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* Y - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCY ) ).
* Before entry, the incremented array Y must contain the n
* element vector y.
* Unchanged on exit.
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array A must contain the upper
* triangular part of the hermitian matrix and the strictly
* lower triangular part of A is not referenced. On exit, the
* upper triangular part of the array A is overwritten by the
* upper triangular part of the updated matrix.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array A must contain the lower
* triangular part of the hermitian matrix and the strictly
* upper triangular part of A is not referenced. On exit, the
* lower triangular part of the array A is overwritten by the
* lower triangular part of the updated matrix.
* Note that the imaginary parts of the diagonal elements need
* not be set, they are assumed to be zero, and on exit they
* are set to zero.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, n ).
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP1, TEMP2
INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX, DBLE
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
$ .NOT.LSAME( UPLO, 'L' ) )THEN
INFO = 1
ELSE IF( N.LT.0 )THEN
INFO = 2
ELSE IF( INCX.EQ.0 )THEN
INFO = 5
ELSE IF( INCY.EQ.0 )THEN
INFO = 7
ELSE IF( LDA.LT.MAX( 1, N ) )THEN
INFO = 9
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZHER2 ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
$ RETURN
*
* Set up the start points in X and Y if the increments are not both
* unity.
*
IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( N - 1 )*INCX
END IF
IF( INCY.GT.0 )THEN
KY = 1
ELSE
KY = 1 - ( N - 1 )*INCY
END IF
JX = KX
JY = KY
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through the triangular part
* of A.
*
IF( LSAME( UPLO, 'U' ) )THEN
*
* Form A when A is stored in the upper triangle.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 20, J = 1, N
IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
TEMP1 = ALPHA*DCONJG( Y( J ) )
TEMP2 = DCONJG( ALPHA*X( J ) )
DO 10, I = 1, J - 1
A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
10 CONTINUE
A( J, J ) = DBLE( A( J, J ) ) +
$ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
ELSE
A( J, J ) = DBLE( A( J, J ) )
END IF
20 CONTINUE
ELSE
DO 40, J = 1, N
IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
TEMP1 = ALPHA*DCONJG( Y( JY ) )
TEMP2 = DCONJG( ALPHA*X( JX ) )
IX = KX
IY = KY
DO 30, I = 1, J - 1
A( I, J ) = A( I, J ) + X( IX )*TEMP1
$ + Y( IY )*TEMP2
IX = IX + INCX
IY = IY + INCY
30 CONTINUE
A( J, J ) = DBLE( A( J, J ) ) +
$ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
ELSE
A( J, J ) = DBLE( A( J, J ) )
END IF
JX = JX + INCX
JY = JY + INCY
40 CONTINUE
END IF
ELSE
*
* Form A when A is stored in the lower triangle.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 60, J = 1, N
IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
TEMP1 = ALPHA*DCONJG( Y( J ) )
TEMP2 = DCONJG( ALPHA*X( J ) )
A( J, J ) = DBLE( A( J, J ) ) +
$ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
DO 50, I = J + 1, N
A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
50 CONTINUE
ELSE
A( J, J ) = DBLE( A( J, J ) )
END IF
60 CONTINUE
ELSE
DO 80, J = 1, N
IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
TEMP1 = ALPHA*DCONJG( Y( JY ) )
TEMP2 = DCONJG( ALPHA*X( JX ) )
A( J, J ) = DBLE( A( J, J ) ) +
$ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
IX = JX
IY = JY
DO 70, I = J + 1, N
IX = IX + INCX
IY = IY + INCY
A( I, J ) = A( I, J ) + X( IX )*TEMP1
$ + Y( IY )*TEMP2
70 CONTINUE
ELSE
A( J, J ) = DBLE( A( J, J ) )
END IF
JX = JX + INCX
JY = JY + INCY
80 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZHER2 .
*
END
SUBROUTINE ZHER2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, BETA,
$ C, LDC )
* .. Scalar Arguments ..
CHARACTER TRANS, UPLO
INTEGER K, LDA, LDB, LDC, N
DOUBLE PRECISION BETA
COMPLEX*16 ALPHA
* ..
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
* ..
*
* Purpose
* =======
*
* ZHER2K performs one of the hermitian rank 2k operations
*
* C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + beta*C,
*
* or
*
* C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + beta*C,
*
* where alpha and beta are scalars with beta real, C is an n by n
* hermitian matrix and A and B are n by k matrices in the first case
* and k by n matrices in the second case.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the array C is to be referenced as
* follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of C
* is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of C
* is to be referenced.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' C := alpha*A*conjg( B' ) +
* conjg( alpha )*B*conjg( A' ) +
* beta*C.
*
* TRANS = 'C' or 'c' C := alpha*conjg( A' )*B +
* conjg( alpha )*conjg( B' )*A +
* beta*C.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix C. N must be
* at least zero.
* Unchanged on exit.
*
* K - INTEGER.
* On entry with TRANS = 'N' or 'n', K specifies the number
* of columns of the matrices A and B, and on entry with
* TRANS = 'C' or 'c', K specifies the number of rows of the
* matrices A and B. K must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
* k when TRANS = 'N' or 'n', and is n otherwise.
* Before entry with TRANS = 'N' or 'n', the leading n by k
* part of the array A must contain the matrix A, otherwise
* the leading k by n part of the array A must contain the
* matrix A.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When TRANS = 'N' or 'n'
* then LDA must be at least max( 1, n ), otherwise LDA must
* be at least max( 1, k ).
* Unchanged on exit.
*
* B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is
* k when TRANS = 'N' or 'n', and is n otherwise.
* Before entry with TRANS = 'N' or 'n', the leading n by k
* part of the array B must contain the matrix B, otherwise
* the leading k by n part of the array B must contain the
* matrix B.
* Unchanged on exit.
*
* LDB - INTEGER.
* On entry, LDB specifies the first dimension of B as declared
* in the calling (sub) program. When TRANS = 'N' or 'n'
* then LDB must be at least max( 1, n ), otherwise LDB must
* be at least max( 1, k ).
* Unchanged on exit.
*
* BETA - DOUBLE PRECISION .
* On entry, BETA specifies the scalar beta.
* Unchanged on exit.
*
* C - COMPLEX*16 array of DIMENSION ( LDC, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array C must contain the upper
* triangular part of the hermitian matrix and the strictly
* lower triangular part of C is not referenced. On exit, the
* upper triangular part of the array C is overwritten by the
* upper triangular part of the updated matrix.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array C must contain the lower
* triangular part of the hermitian matrix and the strictly
* upper triangular part of C is not referenced. On exit, the
* lower triangular part of the array C is overwritten by the
* lower triangular part of the updated matrix.
* Note that the imaginary parts of the diagonal elements need
* not be set, they are assumed to be zero, and on exit they
* are set to zero.
*
* LDC - INTEGER.
* On entry, LDC specifies the first dimension of C as declared
* in the calling (sub) program. LDC must be at least
* max( 1, n ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
* -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1.
* Ed Anderson, Cray Research Inc.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC DBLE, DCONJG, MAX
* ..
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I, INFO, J, L, NROWA
COMPLEX*16 TEMP1, TEMP2
* ..
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
IF( LSAME( TRANS, 'N' ) ) THEN
NROWA = N
ELSE
NROWA = K
END IF
UPPER = LSAME( UPLO, 'U' )
*
INFO = 0
IF( ( .NOT.UPPER ) .AND. ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
INFO = 1
ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ) .AND.
$ ( .NOT.LSAME( TRANS, 'C' ) ) ) THEN
INFO = 2
ELSE IF( N.LT.0 ) THEN
INFO = 3
ELSE IF( K.LT.0 ) THEN
INFO = 4
ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
INFO = 7
ELSE IF( LDB.LT.MAX( 1, NROWA ) ) THEN
INFO = 9
ELSE IF( LDC.LT.MAX( 1, N ) ) THEN
INFO = 12
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZHER2K', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
$ ( BETA.EQ.ONE ) ) )RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO ) THEN
IF( UPPER ) THEN
IF( BETA.EQ.DBLE( ZERO ) ) THEN
DO 20 J = 1, N
DO 10 I = 1, J
C( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40 J = 1, N
DO 30 I = 1, J - 1
C( I, J ) = BETA*C( I, J )
30 CONTINUE
C( J, J ) = BETA*DBLE( C( J, J ) )
40 CONTINUE
END IF
ELSE
IF( BETA.EQ.DBLE( ZERO ) ) THEN
DO 60 J = 1, N
DO 50 I = J, N
C( I, J ) = ZERO
50 CONTINUE
60 CONTINUE
ELSE
DO 80 J = 1, N
C( J, J ) = BETA*DBLE( C( J, J ) )
DO 70 I = J + 1, N
C( I, J ) = BETA*C( I, J )
70 CONTINUE
80 CONTINUE
END IF
END IF
RETURN
END IF
*
* Start the operations.
*
IF( LSAME( TRANS, 'N' ) ) THEN
*
* Form C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) +
* C.
*
IF( UPPER ) THEN
DO 130 J = 1, N
IF( BETA.EQ.DBLE( ZERO ) ) THEN
DO 90 I = 1, J
C( I, J ) = ZERO
90 CONTINUE
ELSE IF( BETA.NE.ONE ) THEN
DO 100 I = 1, J - 1
C( I, J ) = BETA*C( I, J )
100 CONTINUE
C( J, J ) = BETA*DBLE( C( J, J ) )
ELSE
C( J, J ) = DBLE( C( J, J ) )
END IF
DO 120 L = 1, K
IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
$ THEN
TEMP1 = ALPHA*DCONJG( B( J, L ) )
TEMP2 = DCONJG( ALPHA*A( J, L ) )
DO 110 I = 1, J - 1
C( I, J ) = C( I, J ) + A( I, L )*TEMP1 +
$ B( I, L )*TEMP2
110 CONTINUE
C( J, J ) = DBLE( C( J, J ) ) +
$ DBLE( A( J, L )*TEMP1+B( J, L )*TEMP2 )
END IF
120 CONTINUE
130 CONTINUE
ELSE
DO 180 J = 1, N
IF( BETA.EQ.DBLE( ZERO ) ) THEN
DO 140 I = J, N
C( I, J ) = ZERO
140 CONTINUE
ELSE IF( BETA.NE.ONE ) THEN
DO 150 I = J + 1, N
C( I, J ) = BETA*C( I, J )
150 CONTINUE
C( J, J ) = BETA*DBLE( C( J, J ) )
ELSE
C( J, J ) = DBLE( C( J, J ) )
END IF
DO 170 L = 1, K
IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
$ THEN
TEMP1 = ALPHA*DCONJG( B( J, L ) )
TEMP2 = DCONJG( ALPHA*A( J, L ) )
DO 160 I = J + 1, N
C( I, J ) = C( I, J ) + A( I, L )*TEMP1 +
$ B( I, L )*TEMP2
160 CONTINUE
C( J, J ) = DBLE( C( J, J ) ) +
$ DBLE( A( J, L )*TEMP1+B( J, L )*TEMP2 )
END IF
170 CONTINUE
180 CONTINUE
END IF
ELSE
*
* Form C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A +
* C.
*
IF( UPPER ) THEN
DO 210 J = 1, N
DO 200 I = 1, J
TEMP1 = ZERO
TEMP2 = ZERO
DO 190 L = 1, K
TEMP1 = TEMP1 + DCONJG( A( L, I ) )*B( L, J )
TEMP2 = TEMP2 + DCONJG( B( L, I ) )*A( L, J )
190 CONTINUE
IF( I.EQ.J ) THEN
IF( BETA.EQ.DBLE( ZERO ) ) THEN
C( J, J ) = DBLE( ALPHA*TEMP1+DCONJG( ALPHA )*
$ TEMP2 )
ELSE
C( J, J ) = BETA*DBLE( C( J, J ) ) +
$ DBLE( ALPHA*TEMP1+DCONJG( ALPHA )*
$ TEMP2 )
END IF
ELSE
IF( BETA.EQ.DBLE( ZERO ) ) THEN
C( I, J ) = ALPHA*TEMP1 + DCONJG( ALPHA )*TEMP2
ELSE
C( I, J ) = BETA*C( I, J ) + ALPHA*TEMP1 +
$ DCONJG( ALPHA )*TEMP2
END IF
END IF
200 CONTINUE
210 CONTINUE
ELSE
DO 240 J = 1, N
DO 230 I = J, N
TEMP1 = ZERO
TEMP2 = ZERO
DO 220 L = 1, K
TEMP1 = TEMP1 + DCONJG( A( L, I ) )*B( L, J )
TEMP2 = TEMP2 + DCONJG( B( L, I ) )*A( L, J )
220 CONTINUE
IF( I.EQ.J ) THEN
IF( BETA.EQ.DBLE( ZERO ) ) THEN
C( J, J ) = DBLE( ALPHA*TEMP1+DCONJG( ALPHA )*
$ TEMP2 )
ELSE
C( J, J ) = BETA*DBLE( C( J, J ) ) +
$ DBLE( ALPHA*TEMP1+DCONJG( ALPHA )*
$ TEMP2 )
END IF
ELSE
IF( BETA.EQ.DBLE( ZERO ) ) THEN
C( I, J ) = ALPHA*TEMP1 + DCONJG( ALPHA )*TEMP2
ELSE
C( I, J ) = BETA*C( I, J ) + ALPHA*TEMP1 +
$ DCONJG( ALPHA )*TEMP2
END IF
END IF
230 CONTINUE
240 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZHER2K.
*
END
subroutine zscal(n,za,zx,incx)
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
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
subroutine zswap (n,zx,incx,zy,incy)
c
c interchanges two vectors.
c jack dongarra, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
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
SUBROUTINE ZTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
$ B, LDB )
* .. Scalar Arguments ..
CHARACTER*1 SIDE, UPLO, TRANSA, DIAG
INTEGER M, N, LDA, LDB
COMPLEX*16 ALPHA
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), B( LDB, * )
* ..
*
* Purpose
* =======
*
* ZTRMM performs one of the matrix-matrix operations
*
* B := alpha*op( A )*B, or B := alpha*B*op( A )
*
* where alpha is a scalar, B is an m by n matrix, A is a unit, or
* non-unit, upper or lower triangular matrix and op( A ) is one of
*
* op( A ) = A or op( A ) = A' or op( A ) = conjg( A' ).
*
* Parameters
* ==========
*
* SIDE - CHARACTER*1.
* On entry, SIDE specifies whether op( A ) multiplies B from
* the left or right as follows:
*
* SIDE = 'L' or 'l' B := alpha*op( A )*B.
*
* SIDE = 'R' or 'r' B := alpha*B*op( A ).
*
* Unchanged on exit.
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix A is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANSA - CHARACTER*1.
* On entry, TRANSA specifies the form of op( A ) to be used in
* the matrix multiplication as follows:
*
* TRANSA = 'N' or 'n' op( A ) = A.
*
* TRANSA = 'T' or 't' op( A ) = A'.
*
* TRANSA = 'C' or 'c' op( A ) = conjg( A' ).
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit triangular
* as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of B. M must be at
* least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of B. N must be
* at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha. When alpha is
* zero then A is not referenced and B need not be set before
* entry.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, k ), where k is m
* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
* Before entry with UPLO = 'U' or 'u', the leading k by k
* upper triangular part of the array A must contain the upper
* triangular matrix and the strictly lower triangular part of
* A is not referenced.
* Before entry with UPLO = 'L' or 'l', the leading k by k
* lower triangular part of the array A must contain the lower
* triangular matrix and the strictly upper triangular part of
* A is not referenced.
* Note that when DIAG = 'U' or 'u', the diagonal elements of
* A are not referenced either, but are assumed to be unity.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When SIDE = 'L' or 'l' then
* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
* then LDA must be at least max( 1, n ).
* Unchanged on exit.
*
* B - COMPLEX*16 array of DIMENSION ( LDB, n ).
* Before entry, the leading m by n part of the array B must
* contain the matrix B, and on exit is overwritten by the
* transformed matrix.
*
* LDB - INTEGER.
* On entry, LDB specifies the first dimension of B as declared
* in the calling (sub) program. LDB must be at least
* max( 1, m ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX
* .. Local Scalars ..
LOGICAL LSIDE, NOCONJ, NOUNIT, UPPER
INTEGER I, INFO, J, K, NROWA
COMPLEX*16 TEMP
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
LSIDE = LSAME( SIDE , 'L' )
IF( LSIDE )THEN
NROWA = M
ELSE
NROWA = N
END IF
NOCONJ = LSAME( TRANSA, 'T' )
NOUNIT = LSAME( DIAG , 'N' )
UPPER = LSAME( UPLO , 'U' )
*
INFO = 0
IF( ( .NOT.LSIDE ).AND.
$ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.UPPER ).AND.
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
INFO = 2
ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
$ ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
$ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN
INFO = 3
ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND.
$ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN
INFO = 4
ELSE IF( M .LT.0 )THEN
INFO = 5
ELSE IF( N .LT.0 )THEN
INFO = 6
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 9
ELSE IF( LDB.LT.MAX( 1, M ) )THEN
INFO = 11
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZTRMM ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, M
B( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
RETURN
END IF
*
* Start the operations.
*
IF( LSIDE )THEN
IF( LSAME( TRANSA, 'N' ) )THEN
*
* Form B := alpha*A*B.
*
IF( UPPER )THEN
DO 50, J = 1, N
DO 40, K = 1, M
IF( B( K, J ).NE.ZERO )THEN
TEMP = ALPHA*B( K, J )
DO 30, I = 1, K - 1
B( I, J ) = B( I, J ) + TEMP*A( I, K )
30 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP*A( K, K )
B( K, J ) = TEMP
END IF
40 CONTINUE
50 CONTINUE
ELSE
DO 80, J = 1, N
DO 70 K = M, 1, -1
IF( B( K, J ).NE.ZERO )THEN
TEMP = ALPHA*B( K, J )
B( K, J ) = TEMP
IF( NOUNIT )
$ B( K, J ) = B( K, J )*A( K, K )
DO 60, I = K + 1, M
B( I, J ) = B( I, J ) + TEMP*A( I, K )
60 CONTINUE
END IF
70 CONTINUE
80 CONTINUE
END IF
ELSE
*
* Form B := alpha*A'*B or B := alpha*conjg( A' )*B.
*
IF( UPPER )THEN
DO 120, J = 1, N
DO 110, I = M, 1, -1
TEMP = B( I, J )
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*A( I, I )
DO 90, K = 1, I - 1
TEMP = TEMP + A( K, I )*B( K, J )
90 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( A( I, I ) )
DO 100, K = 1, I - 1
TEMP = TEMP + DCONJG( A( K, I ) )*B( K, J )
100 CONTINUE
END IF
B( I, J ) = ALPHA*TEMP
110 CONTINUE
120 CONTINUE
ELSE
DO 160, J = 1, N
DO 150, I = 1, M
TEMP = B( I, J )
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*A( I, I )
DO 130, K = I + 1, M
TEMP = TEMP + A( K, I )*B( K, J )
130 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( A( I, I ) )
DO 140, K = I + 1, M
TEMP = TEMP + DCONJG( A( K, I ) )*B( K, J )
140 CONTINUE
END IF
B( I, J ) = ALPHA*TEMP
150 CONTINUE
160 CONTINUE
END IF
END IF
ELSE
IF( LSAME( TRANSA, 'N' ) )THEN
*
* Form B := alpha*B*A.
*
IF( UPPER )THEN
DO 200, J = N, 1, -1
TEMP = ALPHA
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 170, I = 1, M
B( I, J ) = TEMP*B( I, J )
170 CONTINUE
DO 190, K = 1, J - 1
IF( A( K, J ).NE.ZERO )THEN
TEMP = ALPHA*A( K, J )
DO 180, I = 1, M
B( I, J ) = B( I, J ) + TEMP*B( I, K )
180 CONTINUE
END IF
190 CONTINUE
200 CONTINUE
ELSE
DO 240, J = 1, N
TEMP = ALPHA
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 210, I = 1, M
B( I, J ) = TEMP*B( I, J )
210 CONTINUE
DO 230, K = J + 1, N
IF( A( K, J ).NE.ZERO )THEN
TEMP = ALPHA*A( K, J )
DO 220, I = 1, M
B( I, J ) = B( I, J ) + TEMP*B( I, K )
220 CONTINUE
END IF
230 CONTINUE
240 CONTINUE
END IF
ELSE
*
* Form B := alpha*B*A' or B := alpha*B*conjg( A' ).
*
IF( UPPER )THEN
DO 280, K = 1, N
DO 260, J = 1, K - 1
IF( A( J, K ).NE.ZERO )THEN
IF( NOCONJ )THEN
TEMP = ALPHA*A( J, K )
ELSE
TEMP = ALPHA*DCONJG( A( J, K ) )
END IF
DO 250, I = 1, M
B( I, J ) = B( I, J ) + TEMP*B( I, K )
250 CONTINUE
END IF
260 CONTINUE
TEMP = ALPHA
IF( NOUNIT )THEN
IF( NOCONJ )THEN
TEMP = TEMP*A( K, K )
ELSE
TEMP = TEMP*DCONJG( A( K, K ) )
END IF
END IF
IF( TEMP.NE.ONE )THEN
DO 270, I = 1, M
B( I, K ) = TEMP*B( I, K )
270 CONTINUE
END IF
280 CONTINUE
ELSE
DO 320, K = N, 1, -1
DO 300, J = K + 1, N
IF( A( J, K ).NE.ZERO )THEN
IF( NOCONJ )THEN
TEMP = ALPHA*A( J, K )
ELSE
TEMP = ALPHA*DCONJG( A( J, K ) )
END IF
DO 290, I = 1, M
B( I, J ) = B( I, J ) + TEMP*B( I, K )
290 CONTINUE
END IF
300 CONTINUE
TEMP = ALPHA
IF( NOUNIT )THEN
IF( NOCONJ )THEN
TEMP = TEMP*A( K, K )
ELSE
TEMP = TEMP*DCONJG( A( K, K ) )
END IF
END IF
IF( TEMP.NE.ONE )THEN
DO 310, I = 1, M
B( I, K ) = TEMP*B( I, K )
310 CONTINUE
END IF
320 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ZTRMM .
*
END
SUBROUTINE ZTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, LDA, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * )
* ..
*
* Purpose
* =======
*
* ZTRMV performs one of the matrix-vector operations
*
* x := A*x, or x := A'*x, or x := conjg( A' )*x,
*
* where x is an n element vector and A is an n by n unit, or non-unit,
* upper or lower triangular matrix.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' x := A*x.
*
* TRANS = 'T' or 't' x := A'*x.
*
* TRANS = 'C' or 'c' x := conjg( A' )*x.
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit
* triangular as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array A must contain the upper
* triangular matrix and the strictly lower triangular part of
* A is not referenced.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array A must contain the lower
* triangular matrix and the strictly upper triangular part of
* A is not referenced.
* Note that when DIAG = 'U' or 'u', the diagonal elements of
* A are not referenced either, but are assumed to be unity.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, n ).
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x. On exit, X is overwritten with the
* tranformed vector x.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, J, JX, KX
LOGICAL NOCONJ, NOUNIT
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO , 'U' ).AND.
$ .NOT.LSAME( UPLO , 'L' ) )THEN
INFO = 1
ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
$ .NOT.LSAME( TRANS, 'T' ).AND.
$ .NOT.LSAME( TRANS, 'C' ) )THEN
INFO = 2
ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
$ .NOT.LSAME( DIAG , 'N' ) )THEN
INFO = 3
ELSE IF( N.LT.0 )THEN
INFO = 4
ELSE IF( LDA.LT.MAX( 1, N ) )THEN
INFO = 6
ELSE IF( INCX.EQ.0 )THEN
INFO = 8
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZTRMV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = LSAME( DIAG , 'N' )
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF( INCX.LE.0 )THEN
KX = 1 - ( N - 1 )*INCX
ELSE IF( INCX.NE.1 )THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form x := A*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
IF( INCX.EQ.1 )THEN
DO 20, J = 1, N
IF( X( J ).NE.ZERO )THEN
TEMP = X( J )
DO 10, I = 1, J - 1
X( I ) = X( I ) + TEMP*A( I, J )
10 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )*A( J, J )
END IF
20 CONTINUE
ELSE
JX = KX
DO 40, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = X( JX )
IX = KX
DO 30, I = 1, J - 1
X( IX ) = X( IX ) + TEMP*A( I, J )
IX = IX + INCX
30 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )*A( J, J )
END IF
JX = JX + INCX
40 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 60, J = N, 1, -1
IF( X( J ).NE.ZERO )THEN
TEMP = X( J )
DO 50, I = N, J + 1, -1
X( I ) = X( I ) + TEMP*A( I, J )
50 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )*A( J, J )
END IF
60 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 80, J = N, 1, -1
IF( X( JX ).NE.ZERO )THEN
TEMP = X( JX )
IX = KX
DO 70, I = N, J + 1, -1
X( IX ) = X( IX ) + TEMP*A( I, J )
IX = IX - INCX
70 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )*A( J, J )
END IF
JX = JX - INCX
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := A'*x or x := conjg( A' )*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
IF( INCX.EQ.1 )THEN
DO 110, J = N, 1, -1
TEMP = X( J )
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 90, I = J - 1, 1, -1
TEMP = TEMP + A( I, J )*X( I )
90 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( A( J, J ) )
DO 100, I = J - 1, 1, -1
TEMP = TEMP + DCONJG( A( I, J ) )*X( I )
100 CONTINUE
END IF
X( J ) = TEMP
110 CONTINUE
ELSE
JX = KX + ( N - 1 )*INCX
DO 140, J = N, 1, -1
TEMP = X( JX )
IX = JX
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 120, I = J - 1, 1, -1
IX = IX - INCX
TEMP = TEMP + A( I, J )*X( IX )
120 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( A( J, J ) )
DO 130, I = J - 1, 1, -1
IX = IX - INCX
TEMP = TEMP + DCONJG( A( I, J ) )*X( IX )
130 CONTINUE
END IF
X( JX ) = TEMP
JX = JX - INCX
140 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 170, J = 1, N
TEMP = X( J )
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 150, I = J + 1, N
TEMP = TEMP + A( I, J )*X( I )
150 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( A( J, J ) )
DO 160, I = J + 1, N
TEMP = TEMP + DCONJG( A( I, J ) )*X( I )
160 CONTINUE
END IF
X( J ) = TEMP
170 CONTINUE
ELSE
JX = KX
DO 200, J = 1, N
TEMP = X( JX )
IX = JX
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 180, I = J + 1, N
IX = IX + INCX
TEMP = TEMP + A( I, J )*X( IX )
180 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( A( J, J ) )
DO 190, I = J + 1, N
IX = IX + INCX
TEMP = TEMP + DCONJG( A( I, J ) )*X( IX )
190 CONTINUE
END IF
X( JX ) = TEMP
JX = JX + INCX
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ZTRMV .
*
END
SUBROUTINE ZTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
$ B, LDB )
* .. Scalar Arguments ..
CHARACTER*1 SIDE, UPLO, TRANSA, DIAG
INTEGER M, N, LDA, LDB
COMPLEX*16 ALPHA
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), B( LDB, * )
* ..
*
* Purpose
* =======
*
* ZTRSM solves one of the matrix equations
*
* op( A )*X = alpha*B, or X*op( A ) = alpha*B,
*
* where alpha is a scalar, X and B are m by n matrices, A is a unit, or
* non-unit, upper or lower triangular matrix and op( A ) is one of
*
* op( A ) = A or op( A ) = A' or op( A ) = conjg( A' ).
*
* The matrix X is overwritten on B.
*
* Parameters
* ==========
*
* SIDE - CHARACTER*1.
* On entry, SIDE specifies whether op( A ) appears on the left
* or right of X as follows:
*
* SIDE = 'L' or 'l' op( A )*X = alpha*B.
*
* SIDE = 'R' or 'r' X*op( A ) = alpha*B.
*
* Unchanged on exit.
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix A is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANSA - CHARACTER*1.
* On entry, TRANSA specifies the form of op( A ) to be used in
* the matrix multiplication as follows:
*
* TRANSA = 'N' or 'n' op( A ) = A.
*
* TRANSA = 'T' or 't' op( A ) = A'.
*
* TRANSA = 'C' or 'c' op( A ) = conjg( A' ).
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit triangular
* as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of B. M must be at
* least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of B. N must be
* at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha. When alpha is
* zero then A is not referenced and B need not be set before
* entry.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, k ), where k is m
* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
* Before entry with UPLO = 'U' or 'u', the leading k by k
* upper triangular part of the array A must contain the upper
* triangular matrix and the strictly lower triangular part of
* A is not referenced.
* Before entry with UPLO = 'L' or 'l', the leading k by k
* lower triangular part of the array A must contain the lower
* triangular matrix and the strictly upper triangular part of
* A is not referenced.
* Note that when DIAG = 'U' or 'u', the diagonal elements of
* A are not referenced either, but are assumed to be unity.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When SIDE = 'L' or 'l' then
* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
* then LDA must be at least max( 1, n ).
* Unchanged on exit.
*
* B - COMPLEX*16 array of DIMENSION ( LDB, n ).
* Before entry, the leading m by n part of the array B must
* contain the right-hand side matrix B, and on exit is
* overwritten by the solution matrix X.
*
* LDB - INTEGER.
* On entry, LDB specifies the first dimension of B as declared
* in the calling (sub) program. LDB must be at least
* max( 1, m ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX
* .. Local Scalars ..
LOGICAL LSIDE, NOCONJ, NOUNIT, UPPER
INTEGER I, INFO, J, K, NROWA
COMPLEX*16 TEMP
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
LSIDE = LSAME( SIDE , 'L' )
IF( LSIDE )THEN
NROWA = M
ELSE
NROWA = N
END IF
NOCONJ = LSAME( TRANSA, 'T' )
NOUNIT = LSAME( DIAG , 'N' )
UPPER = LSAME( UPLO , 'U' )
*
INFO = 0
IF( ( .NOT.LSIDE ).AND.
$ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.UPPER ).AND.
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
INFO = 2
ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
$ ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
$ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN
INFO = 3
ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND.
$ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN
INFO = 4
ELSE IF( M .LT.0 )THEN
INFO = 5
ELSE IF( N .LT.0 )THEN
INFO = 6
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 9
ELSE IF( LDB.LT.MAX( 1, M ) )THEN
INFO = 11
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZTRSM ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, M
B( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
RETURN
END IF
*
* Start the operations.
*
IF( LSIDE )THEN
IF( LSAME( TRANSA, 'N' ) )THEN
*
* Form B := alpha*inv( A )*B.
*
IF( UPPER )THEN
DO 60, J = 1, N
IF( ALPHA.NE.ONE )THEN
DO 30, I = 1, M
B( I, J ) = ALPHA*B( I, J )
30 CONTINUE
END IF
DO 50, K = M, 1, -1
IF( B( K, J ).NE.ZERO )THEN
IF( NOUNIT )
$ B( K, J ) = B( K, J )/A( K, K )
DO 40, I = 1, K - 1
B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
40 CONTINUE
END IF
50 CONTINUE
60 CONTINUE
ELSE
DO 100, J = 1, N
IF( ALPHA.NE.ONE )THEN
DO 70, I = 1, M
B( I, J ) = ALPHA*B( I, J )
70 CONTINUE
END IF
DO 90 K = 1, M
IF( B( K, J ).NE.ZERO )THEN
IF( NOUNIT )
$ B( K, J ) = B( K, J )/A( K, K )
DO 80, I = K + 1, M
B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
80 CONTINUE
END IF
90 CONTINUE
100 CONTINUE
END IF
ELSE
*
* Form B := alpha*inv( A' )*B
* or B := alpha*inv( conjg( A' ) )*B.
*
IF( UPPER )THEN
DO 140, J = 1, N
DO 130, I = 1, M
TEMP = ALPHA*B( I, J )
IF( NOCONJ )THEN
DO 110, K = 1, I - 1
TEMP = TEMP - A( K, I )*B( K, J )
110 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( I, I )
ELSE
DO 120, K = 1, I - 1
TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J )
120 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( A( I, I ) )
END IF
B( I, J ) = TEMP
130 CONTINUE
140 CONTINUE
ELSE
DO 180, J = 1, N
DO 170, I = M, 1, -1
TEMP = ALPHA*B( I, J )
IF( NOCONJ )THEN
DO 150, K = I + 1, M
TEMP = TEMP - A( K, I )*B( K, J )
150 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( I, I )
ELSE
DO 160, K = I + 1, M
TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J )
160 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( A( I, I ) )
END IF
B( I, J ) = TEMP
170 CONTINUE
180 CONTINUE
END IF
END IF
ELSE
IF( LSAME( TRANSA, 'N' ) )THEN
*
* Form B := alpha*B*inv( A ).
*
IF( UPPER )THEN
DO 230, J = 1, N
IF( ALPHA.NE.ONE )THEN
DO 190, I = 1, M
B( I, J ) = ALPHA*B( I, J )
190 CONTINUE
END IF
DO 210, K = 1, J - 1
IF( A( K, J ).NE.ZERO )THEN
DO 200, I = 1, M
B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
200 CONTINUE
END IF
210 CONTINUE
IF( NOUNIT )THEN
TEMP = ONE/A( J, J )
DO 220, I = 1, M
B( I, J ) = TEMP*B( I, J )
220 CONTINUE
END IF
230 CONTINUE
ELSE
DO 280, J = N, 1, -1
IF( ALPHA.NE.ONE )THEN
DO 240, I = 1, M
B( I, J ) = ALPHA*B( I, J )
240 CONTINUE
END IF
DO 260, K = J + 1, N
IF( A( K, J ).NE.ZERO )THEN
DO 250, I = 1, M
B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
250 CONTINUE
END IF
260 CONTINUE
IF( NOUNIT )THEN
TEMP = ONE/A( J, J )
DO 270, I = 1, M
B( I, J ) = TEMP*B( I, J )
270 CONTINUE
END IF
280 CONTINUE
END IF
ELSE
*
* Form B := alpha*B*inv( A' )
* or B := alpha*B*inv( conjg( A' ) ).
*
IF( UPPER )THEN
DO 330, K = N, 1, -1
IF( NOUNIT )THEN
IF( NOCONJ )THEN
TEMP = ONE/A( K, K )
ELSE
TEMP = ONE/DCONJG( A( K, K ) )
END IF
DO 290, I = 1, M
B( I, K ) = TEMP*B( I, K )
290 CONTINUE
END IF
DO 310, J = 1, K - 1
IF( A( J, K ).NE.ZERO )THEN
IF( NOCONJ )THEN
TEMP = A( J, K )
ELSE
TEMP = DCONJG( A( J, K ) )
END IF
DO 300, I = 1, M
B( I, J ) = B( I, J ) - TEMP*B( I, K )
300 CONTINUE
END IF
310 CONTINUE
IF( ALPHA.NE.ONE )THEN
DO 320, I = 1, M
B( I, K ) = ALPHA*B( I, K )
320 CONTINUE
END IF
330 CONTINUE
ELSE
DO 380, K = 1, N
IF( NOUNIT )THEN
IF( NOCONJ )THEN
TEMP = ONE/A( K, K )
ELSE
TEMP = ONE/DCONJG( A( K, K ) )
END IF
DO 340, I = 1, M
B( I, K ) = TEMP*B( I, K )
340 CONTINUE
END IF
DO 360, J = K + 1, N
IF( A( J, K ).NE.ZERO )THEN
IF( NOCONJ )THEN
TEMP = A( J, K )
ELSE
TEMP = DCONJG( A( J, K ) )
END IF
DO 350, I = 1, M
B( I, J ) = B( I, J ) - TEMP*B( I, K )
350 CONTINUE
END IF
360 CONTINUE
IF( ALPHA.NE.ONE )THEN
DO 370, I = 1, M
B( I, K ) = ALPHA*B( I, K )
370 CONTINUE
END IF
380 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ZTRSM .
*
END
SUBROUTINE ZTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, LDA, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * )
* ..
*
* Purpose
* =======
*
* ZTRSV solves one of the systems of equations
*
* A*x = b, or A'*x = b, or conjg( A' )*x = b,
*
* where b and x are n element vectors and A is an n by n unit, or
* non-unit, upper or lower triangular matrix.
*
* No test for singularity or near-singularity is included in this
* routine. Such tests must be performed before calling this routine.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the equations to be solved as
* follows:
*
* TRANS = 'N' or 'n' A*x = b.
*
* TRANS = 'T' or 't' A'*x = b.
*
* TRANS = 'C' or 'c' conjg( A' )*x = b.
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit
* triangular as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array A must contain the upper
* triangular matrix and the strictly lower triangular part of
* A is not referenced.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array A must contain the lower
* triangular matrix and the strictly upper triangular part of
* A is not referenced.
* Note that when DIAG = 'U' or 'u', the diagonal elements of
* A are not referenced either, but are assumed to be unity.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, n ).
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element right-hand side vector b. On exit, X is overwritten
* with the solution vector x.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, J, JX, KX
LOGICAL NOCONJ, NOUNIT
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO , 'U' ).AND.
$ .NOT.LSAME( UPLO , 'L' ) )THEN
INFO = 1
ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
$ .NOT.LSAME( TRANS, 'T' ).AND.
$ .NOT.LSAME( TRANS, 'C' ) )THEN
INFO = 2
ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
$ .NOT.LSAME( DIAG , 'N' ) )THEN
INFO = 3
ELSE IF( N.LT.0 )THEN
INFO = 4
ELSE IF( LDA.LT.MAX( 1, N ) )THEN
INFO = 6
ELSE IF( INCX.EQ.0 )THEN
INFO = 8
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZTRSV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = LSAME( DIAG , 'N' )
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF( INCX.LE.0 )THEN
KX = 1 - ( N - 1 )*INCX
ELSE IF( INCX.NE.1 )THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form x := inv( A )*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
IF( INCX.EQ.1 )THEN
DO 20, J = N, 1, -1
IF( X( J ).NE.ZERO )THEN
IF( NOUNIT )
$ X( J ) = X( J )/A( J, J )
TEMP = X( J )
DO 10, I = J - 1, 1, -1
X( I ) = X( I ) - TEMP*A( I, J )
10 CONTINUE
END IF
20 CONTINUE
ELSE
JX = KX + ( N - 1 )*INCX
DO 40, J = N, 1, -1
IF( X( JX ).NE.ZERO )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )/A( J, J )
TEMP = X( JX )
IX = JX
DO 30, I = J - 1, 1, -1
IX = IX - INCX
X( IX ) = X( IX ) - TEMP*A( I, J )
30 CONTINUE
END IF
JX = JX - INCX
40 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 60, J = 1, N
IF( X( J ).NE.ZERO )THEN
IF( NOUNIT )
$ X( J ) = X( J )/A( J, J )
TEMP = X( J )
DO 50, I = J + 1, N
X( I ) = X( I ) - TEMP*A( I, J )
50 CONTINUE
END IF
60 CONTINUE
ELSE
JX = KX
DO 80, J = 1, N
IF( X( JX ).NE.ZERO )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )/A( J, J )
TEMP = X( JX )
IX = JX
DO 70, I = J + 1, N
IX = IX + INCX
X( IX ) = X( IX ) - TEMP*A( I, J )
70 CONTINUE
END IF
JX = JX + INCX
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := inv( A' )*x or x := inv( conjg( A' ) )*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
IF( INCX.EQ.1 )THEN
DO 110, J = 1, N
TEMP = X( J )
IF( NOCONJ )THEN
DO 90, I = 1, J - 1
TEMP = TEMP - A( I, J )*X( I )
90 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( J, J )
ELSE
DO 100, I = 1, J - 1
TEMP = TEMP - DCONJG( A( I, J ) )*X( I )
100 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( A( J, J ) )
END IF
X( J ) = TEMP
110 CONTINUE
ELSE
JX = KX
DO 140, J = 1, N
IX = KX
TEMP = X( JX )
IF( NOCONJ )THEN
DO 120, I = 1, J - 1
TEMP = TEMP - A( I, J )*X( IX )
IX = IX + INCX
120 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( J, J )
ELSE
DO 130, I = 1, J - 1
TEMP = TEMP - DCONJG( A( I, J ) )*X( IX )
IX = IX + INCX
130 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( A( J, J ) )
END IF
X( JX ) = TEMP
JX = JX + INCX
140 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 170, J = N, 1, -1
TEMP = X( J )
IF( NOCONJ )THEN
DO 150, I = N, J + 1, -1
TEMP = TEMP - A( I, J )*X( I )
150 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( J, J )
ELSE
DO 160, I = N, J + 1, -1
TEMP = TEMP - DCONJG( A( I, J ) )*X( I )
160 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( A( J, J ) )
END IF
X( J ) = TEMP
170 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 200, J = N, 1, -1
IX = KX
TEMP = X( JX )
IF( NOCONJ )THEN
DO 180, I = N, J + 1, -1
TEMP = TEMP - A( I, J )*X( IX )
IX = IX - INCX
180 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( J, J )
ELSE
DO 190, I = N, J + 1, -1
TEMP = TEMP - DCONJG( A( I, J ) )*X( IX )
IX = IX - INCX
190 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( A( J, J ) )
END IF
X( JX ) = TEMP
JX = JX - INCX
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ZTRSV .
*
END
subroutine zdrot (n,zx,incx,zy,incy,c,s)
c
c applies a plane rotation, where the cos and sin (c and s) are
c double precision and the vectors zx and zy are double complex.
c jack dongarra, linpack, 3/11/78.
c
double complex zx(1),zy(1),ztemp
double precision 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
ztemp = c*zx(ix) + s*zy(iy)
zy(iy) = c*zy(iy) - s*zx(ix)
zx(ix) = ztemp
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
ztemp = c*zx(i) + s*zy(i)
zy(i) = c*zy(i) - s*zx(i)
zx(i) = ztemp
30 continue
return
end
SUBROUTINE ZGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
$ BETA, Y, INCY )
* .. Scalar Arguments ..
COMPLEX*16 ALPHA, BETA
INTEGER INCX, INCY, KL, KU, LDA, M, N
CHARACTER*1 TRANS
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * ), Y( * )
* ..
*
* Purpose
* =======
*
* ZGBMV performs one of the matrix-vector operations
*
* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, or
*
* y := alpha*conjg( A' )*x + beta*y,
*
* where alpha and beta are scalars, x and y are vectors and A is an
* m by n band matrix, with kl sub-diagonals and ku super-diagonals.
*
* Parameters
* ==========
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
*
* TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
*
* TRANS = 'C' or 'c' y := alpha*conjg( A' )*x + beta*y.
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of the matrix A.
* M must be at least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* KL - INTEGER.
* On entry, KL specifies the number of sub-diagonals of the
* matrix A. KL must satisfy 0 .le. KL.
* Unchanged on exit.
*
* KU - INTEGER.
* On entry, KU specifies the number of super-diagonals of the
* matrix A. KU must satisfy 0 .le. KU.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry, the leading ( kl + ku + 1 ) by n part of the
* array A must contain the matrix of coefficients, supplied
* column by column, with the leading diagonal of the matrix in
* row ( ku + 1 ) of the array, the first super-diagonal
* starting at position 2 in row ku, the first sub-diagonal
* starting at position 1 in row ( ku + 2 ), and so on.
* Elements in the array A that do not correspond to elements
* in the band matrix (such as the top left ku by ku triangle)
* are not referenced.
* The following program segment will transfer a band matrix
* from conventional full matrix storage to band storage:
*
* DO 20, J = 1, N
* K = KU + 1 - J
* DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
* A( K + I, J ) = matrix( I, J )
* 10 CONTINUE
* 20 CONTINUE
*
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* ( kl + ku + 1 ).
* Unchanged on exit.
*
* X - COMPLEX*16 array of DIMENSION at least
* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
* and at least
* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
* Before entry, the incremented array X must contain the
* vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* BETA - COMPLEX*16 .
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then Y need not be set on input.
* Unchanged on exit.
*
* Y - COMPLEX*16 array of DIMENSION at least
* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
* and at least
* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
* Before entry, the incremented array Y must contain the
* vector y. On exit, Y is overwritten by the updated vector y.
*
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, IY, J, JX, JY, K, KUP1, KX, KY,
$ LENX, LENY
LOGICAL NOCONJ
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( TRANS, 'N' ).AND.
$ .NOT.LSAME( TRANS, 'T' ).AND.
$ .NOT.LSAME( TRANS, 'C' ) )THEN
INFO = 1
ELSE IF( M.LT.0 )THEN
INFO = 2
ELSE IF( N.LT.0 )THEN
INFO = 3
ELSE IF( KL.LT.0 )THEN
INFO = 4
ELSE IF( KU.LT.0 )THEN
INFO = 5
ELSE IF( LDA.LT.( KL + KU + 1 ) )THEN
INFO = 8
ELSE IF( INCX.EQ.0 )THEN
INFO = 10
ELSE IF( INCY.EQ.0 )THEN
INFO = 13
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZGBMV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
$ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
*
* Set LENX and LENY, the lengths of the vectors x and y, and set
* up the start points in X and Y.
*
IF( LSAME( TRANS, 'N' ) )THEN
LENX = N
LENY = M
ELSE
LENX = M
LENY = N
END IF
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( LENX - 1 )*INCX
END IF
IF( INCY.GT.0 )THEN
KY = 1
ELSE
KY = 1 - ( LENY - 1 )*INCY
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through the band part of A.
*
* First form y := beta*y.
*
IF( BETA.NE.ONE )THEN
IF( INCY.EQ.1 )THEN
IF( BETA.EQ.ZERO )THEN
DO 10, I = 1, LENY
Y( I ) = ZERO
10 CONTINUE
ELSE
DO 20, I = 1, LENY
Y( I ) = BETA*Y( I )
20 CONTINUE
END IF
ELSE
IY = KY
IF( BETA.EQ.ZERO )THEN
DO 30, I = 1, LENY
Y( IY ) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40, I = 1, LENY
Y( IY ) = BETA*Y( IY )
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF( ALPHA.EQ.ZERO )
$ RETURN
KUP1 = KU + 1
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form y := alpha*A*x + y.
*
JX = KX
IF( INCY.EQ.1 )THEN
DO 60, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*X( JX )
K = KUP1 - J
DO 50, I = MAX( 1, J - KU ), MIN( M, J + KL )
Y( I ) = Y( I ) + TEMP*A( K + I, J )
50 CONTINUE
END IF
JX = JX + INCX
60 CONTINUE
ELSE
DO 80, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*X( JX )
IY = KY
K = KUP1 - J
DO 70, I = MAX( 1, J - KU ), MIN( M, J + KL )
Y( IY ) = Y( IY ) + TEMP*A( K + I, J )
IY = IY + INCY
70 CONTINUE
END IF
JX = JX + INCX
IF( J.GT.KU )
$ KY = KY + INCY
80 CONTINUE
END IF
ELSE
*
* Form y := alpha*A'*x + y or y := alpha*conjg( A' )*x + y.
*
JY = KY
IF( INCX.EQ.1 )THEN
DO 110, J = 1, N
TEMP = ZERO
K = KUP1 - J
IF( NOCONJ )THEN
DO 90, I = MAX( 1, J - KU ), MIN( M, J + KL )
TEMP = TEMP + A( K + I, J )*X( I )
90 CONTINUE
ELSE
DO 100, I = MAX( 1, J - KU ), MIN( M, J + KL )
TEMP = TEMP + DCONJG( A( K + I, J ) )*X( I )
100 CONTINUE
END IF
Y( JY ) = Y( JY ) + ALPHA*TEMP
JY = JY + INCY
110 CONTINUE
ELSE
DO 140, J = 1, N
TEMP = ZERO
IX = KX
K = KUP1 - J
IF( NOCONJ )THEN
DO 120, I = MAX( 1, J - KU ), MIN( M, J + KL )
TEMP = TEMP + A( K + I, J )*X( IX )
IX = IX + INCX
120 CONTINUE
ELSE
DO 130, I = MAX( 1, J - KU ), MIN( M, J + KL )
TEMP = TEMP + DCONJG( A( K + I, J ) )*X( IX )
IX = IX + INCX
130 CONTINUE
END IF
Y( JY ) = Y( JY ) + ALPHA*TEMP
JY = JY + INCY
IF( J.GT.KU )
$ KX = KX + INCX
140 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZGBMV .
*
END
SUBROUTINE ZGERU ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
* .. Scalar Arguments ..
COMPLEX*16 ALPHA
INTEGER INCX, INCY, LDA, M, N
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * ), Y( * )
* ..
*
* Purpose
* =======
*
* ZGERU performs the rank 1 operation
*
* A := alpha*x*y' + A,
*
* where alpha is a scalar, x is an m element vector, y is an n element
* vector and A is an m by n matrix.
*
* Parameters
* ==========
*
* M - INTEGER.
* On entry, M specifies the number of rows of the matrix A.
* M must be at least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( m - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the m
* element vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* Y - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCY ) ).
* Before entry, the incremented array Y must contain the n
* element vector y.
* Unchanged on exit.
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry, the leading m by n part of the array A must
* contain the matrix of coefficients. On exit, A is
* overwritten by the updated matrix.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, m ).
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, J, JY, KX
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( M.LT.0 )THEN
INFO = 1
ELSE IF( N.LT.0 )THEN
INFO = 2
ELSE IF( INCX.EQ.0 )THEN
INFO = 5
ELSE IF( INCY.EQ.0 )THEN
INFO = 7
ELSE IF( LDA.LT.MAX( 1, M ) )THEN
INFO = 9
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZGERU ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
$ RETURN
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF( INCY.GT.0 )THEN
JY = 1
ELSE
JY = 1 - ( N - 1 )*INCY
END IF
IF( INCX.EQ.1 )THEN
DO 20, J = 1, N
IF( Y( JY ).NE.ZERO )THEN
TEMP = ALPHA*Y( JY )
DO 10, I = 1, M
A( I, J ) = A( I, J ) + X( I )*TEMP
10 CONTINUE
END IF
JY = JY + INCY
20 CONTINUE
ELSE
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( M - 1 )*INCX
END IF
DO 40, J = 1, N
IF( Y( JY ).NE.ZERO )THEN
TEMP = ALPHA*Y( JY )
IX = KX
DO 30, I = 1, M
A( I, J ) = A( I, J ) + X( IX )*TEMP
IX = IX + INCX
30 CONTINUE
END IF
JY = JY + INCY
40 CONTINUE
END IF
*
RETURN
*
* End of ZGERU .
*
END
SUBROUTINE ZHBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX,
$ BETA, Y, INCY )
* .. Scalar Arguments ..
COMPLEX*16 ALPHA, BETA
INTEGER INCX, INCY, K, LDA, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * ), Y( * )
* ..
*
* Purpose
* =======
*
* ZHBMV performs the matrix-vector operation
*
* y := alpha*A*x + beta*y,
*
* where alpha and beta are scalars, x and y are n element vectors and
* A is an n by n hermitian band matrix, with k super-diagonals.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the band matrix A is being supplied as
* follows:
*
* UPLO = 'U' or 'u' The upper triangular part of A is
* being supplied.
*
* UPLO = 'L' or 'l' The lower triangular part of A is
* being supplied.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* K - INTEGER.
* On entry, K specifies the number of super-diagonals of the
* matrix A. K must satisfy 0 .le. K.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
* by n part of the array A must contain the upper triangular
* band part of the hermitian matrix, supplied column by
* column, with the leading diagonal of the matrix in row
* ( k + 1 ) of the array, the first super-diagonal starting at
* position 2 in row k, and so on. The top left k by k triangle
* of the array A is not referenced.
* The following program segment will transfer the upper
* triangular part of a hermitian band matrix from conventional
* full matrix storage to band storage:
*
* DO 20, J = 1, N
* M = K + 1 - J
* DO 10, I = MAX( 1, J - K ), J
* A( M + I, J ) = matrix( I, J )
* 10 CONTINUE
* 20 CONTINUE
*
* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
* by n part of the array A must contain the lower triangular
* band part of the hermitian matrix, supplied column by
* column, with the leading diagonal of the matrix in row 1 of
* the array, the first sub-diagonal starting at position 1 in
* row 2, and so on. The bottom right k by k triangle of the
* array A is not referenced.
* The following program segment will transfer the lower
* triangular part of a hermitian band matrix from conventional
* full matrix storage to band storage:
*
* DO 20, J = 1, N
* M = 1 - J
* DO 10, I = J, MIN( N, J + K )
* A( M + I, J ) = matrix( I, J )
* 10 CONTINUE
* 20 CONTINUE
*
* Note that the imaginary parts of the diagonal elements need
* not be set and are assumed to be zero.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* ( k + 1 ).
* Unchanged on exit.
*
* X - COMPLEX*16 array of DIMENSION at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the
* vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* BETA - COMPLEX*16 .
* On entry, BETA specifies the scalar beta.
* Unchanged on exit.
*
* Y - COMPLEX*16 array of DIMENSION at least
* ( 1 + ( n - 1 )*abs( INCY ) ).
* Before entry, the incremented array Y must contain the
* vector y. On exit, Y is overwritten by the updated vector y.
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP1, TEMP2
INTEGER I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX, MIN, DBLE
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
$ .NOT.LSAME( UPLO, 'L' ) )THEN
INFO = 1
ELSE IF( N.LT.0 )THEN
INFO = 2
ELSE IF( K.LT.0 )THEN
INFO = 3
ELSE IF( LDA.LT.( K + 1 ) )THEN
INFO = 6
ELSE IF( INCX.EQ.0 )THEN
INFO = 8
ELSE IF( INCY.EQ.0 )THEN
INFO = 11
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZHBMV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* Set up the start points in X and Y.
*
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( N - 1 )*INCX
END IF
IF( INCY.GT.0 )THEN
KY = 1
ELSE
KY = 1 - ( N - 1 )*INCY
END IF
*
* Start the operations. In this version the elements of the array A
* are accessed sequentially with one pass through A.
*
* First form y := beta*y.
*
IF( BETA.NE.ONE )THEN
IF( INCY.EQ.1 )THEN
IF( BETA.EQ.ZERO )THEN
DO 10, I = 1, N
Y( I ) = ZERO
10 CONTINUE
ELSE
DO 20, I = 1, N
Y( I ) = BETA*Y( I )
20 CONTINUE
END IF
ELSE
IY = KY
IF( BETA.EQ.ZERO )THEN
DO 30, I = 1, N
Y( IY ) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40, I = 1, N
Y( IY ) = BETA*Y( IY )
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF( ALPHA.EQ.ZERO )
$ RETURN
IF( LSAME( UPLO, 'U' ) )THEN
*
* Form y when upper triangle of A is stored.
*
KPLUS1 = K + 1
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 60, J = 1, N
TEMP1 = ALPHA*X( J )
TEMP2 = ZERO
L = KPLUS1 - J
DO 50, I = MAX( 1, J - K ), J - 1
Y( I ) = Y( I ) + TEMP1*A( L + I, J )
TEMP2 = TEMP2 + DCONJG( A( L + I, J ) )*X( I )
50 CONTINUE
Y( J ) = Y( J ) + TEMP1*DBLE( A( KPLUS1, J ) )
$ + ALPHA*TEMP2
60 CONTINUE
ELSE
JX = KX
JY = KY
DO 80, J = 1, N
TEMP1 = ALPHA*X( JX )
TEMP2 = ZERO
IX = KX
IY = KY
L = KPLUS1 - J
DO 70, I = MAX( 1, J - K ), J - 1
Y( IY ) = Y( IY ) + TEMP1*A( L + I, J )
TEMP2 = TEMP2 + DCONJG( A( L + I, J ) )*X( IX )
IX = IX + INCX
IY = IY + INCY
70 CONTINUE
Y( JY ) = Y( JY ) + TEMP1*DBLE( A( KPLUS1, J ) )
$ + ALPHA*TEMP2
JX = JX + INCX
JY = JY + INCY
IF( J.GT.K )THEN
KX = KX + INCX
KY = KY + INCY
END IF
80 CONTINUE
END IF
ELSE
*
* Form y when lower triangle of A is stored.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 100, J = 1, N
TEMP1 = ALPHA*X( J )
TEMP2 = ZERO
Y( J ) = Y( J ) + TEMP1*DBLE( A( 1, J ) )
L = 1 - J
DO 90, I = J + 1, MIN( N, J + K )
Y( I ) = Y( I ) + TEMP1*A( L + I, J )
TEMP2 = TEMP2 + DCONJG( A( L + I, J ) )*X( I )
90 CONTINUE
Y( J ) = Y( J ) + ALPHA*TEMP2
100 CONTINUE
ELSE
JX = KX
JY = KY
DO 120, J = 1, N
TEMP1 = ALPHA*X( JX )
TEMP2 = ZERO
Y( JY ) = Y( JY ) + TEMP1*DBLE( A( 1, J ) )
L = 1 - J
IX = JX
IY = JY
DO 110, I = J + 1, MIN( N, J + K )
IX = IX + INCX
IY = IY + INCY
Y( IY ) = Y( IY ) + TEMP1*A( L + I, J )
TEMP2 = TEMP2 + DCONJG( A( L + I, J ) )*X( IX )
110 CONTINUE
Y( JY ) = Y( JY ) + ALPHA*TEMP2
JX = JX + INCX
JY = JY + INCY
120 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZHBMV .
*
END
SUBROUTINE ZHEMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB,
$ BETA, C, LDC )
* .. Scalar Arguments ..
CHARACTER*1 SIDE, UPLO
INTEGER M, N, LDA, LDB, LDC
COMPLEX*16 ALPHA, BETA
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
* ..
*
* Purpose
* =======
*
* ZHEMM performs one of the matrix-matrix operations
*
* C := alpha*A*B + beta*C,
*
* or
*
* C := alpha*B*A + beta*C,
*
* where alpha and beta are scalars, A is an hermitian matrix and B and
* C are m by n matrices.
*
* Parameters
* ==========
*
* SIDE - CHARACTER*1.
* On entry, SIDE specifies whether the hermitian matrix A
* appears on the left or right in the operation as follows:
*
* SIDE = 'L' or 'l' C := alpha*A*B + beta*C,
*
* SIDE = 'R' or 'r' C := alpha*B*A + beta*C,
*
* Unchanged on exit.
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the hermitian matrix A is to be
* referenced as follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of the
* hermitian matrix is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of the
* hermitian matrix is to be referenced.
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of the matrix C.
* M must be at least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of the matrix C.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
* m when SIDE = 'L' or 'l' and is n otherwise.
* Before entry with SIDE = 'L' or 'l', the m by m part of
* the array A must contain the hermitian matrix, such that
* when UPLO = 'U' or 'u', the leading m by m upper triangular
* part of the array A must contain the upper triangular part
* of the hermitian matrix and the strictly lower triangular
* part of A is not referenced, and when UPLO = 'L' or 'l',
* the leading m by m lower triangular part of the array A
* must contain the lower triangular part of the hermitian
* matrix and the strictly upper triangular part of A is not
* referenced.
* Before entry with SIDE = 'R' or 'r', the n by n part of
* the array A must contain the hermitian matrix, such that
* when UPLO = 'U' or 'u', the leading n by n upper triangular
* part of the array A must contain the upper triangular part
* of the hermitian matrix and the strictly lower triangular
* part of A is not referenced, and when UPLO = 'L' or 'l',
* the leading n by n lower triangular part of the array A
* must contain the lower triangular part of the hermitian
* matrix and the strictly upper triangular part of A is not
* referenced.
* Note that the imaginary parts of the diagonal elements need
* not be set, they are assumed to be zero.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When SIDE = 'L' or 'l' then
* LDA must be at least max( 1, m ), otherwise LDA must be at
* least max( 1, n ).
* Unchanged on exit.
*
* B - COMPLEX*16 array of DIMENSION ( LDB, n ).
* Before entry, the leading m by n part of the array B must
* contain the matrix B.
* Unchanged on exit.
*
* LDB - INTEGER.
* On entry, LDB specifies the first dimension of B as declared
* in the calling (sub) program. LDB must be at least
* max( 1, m ).
* Unchanged on exit.
*
* BETA - COMPLEX*16 .
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then C need not be set on input.
* Unchanged on exit.
*
* C - COMPLEX*16 array of DIMENSION ( LDC, n ).
* Before entry, the leading m by n part of the array C must
* contain the matrix C, except when beta is zero, in which
* case C need not be set on entry.
* On exit, the array C is overwritten by the m by n updated
* matrix.
*
* LDC - INTEGER.
* On entry, LDC specifies the first dimension of C as declared
* in the calling (sub) program. LDC must be at least
* max( 1, m ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX, DBLE
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I, INFO, J, K, NROWA
COMPLEX*16 TEMP1, TEMP2
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* ..
* .. Executable Statements ..
*
* Set NROWA as the number of rows of A.
*
IF( LSAME( SIDE, 'L' ) )THEN
NROWA = M
ELSE
NROWA = N
END IF
UPPER = LSAME( UPLO, 'U' )
*
* Test the input parameters.
*
INFO = 0
IF( ( .NOT.LSAME( SIDE, 'L' ) ).AND.
$ ( .NOT.LSAME( SIDE, 'R' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.UPPER ).AND.
$ ( .NOT.LSAME( UPLO, 'L' ) ) )THEN
INFO = 2
ELSE IF( M .LT.0 )THEN
INFO = 3
ELSE IF( N .LT.0 )THEN
INFO = 4
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 7
ELSE IF( LDB.LT.MAX( 1, M ) )THEN
INFO = 9
ELSE IF( LDC.LT.MAX( 1, M ) )THEN
INFO = 12
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZHEMM ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
$ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
IF( BETA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, M
C( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40, J = 1, N
DO 30, I = 1, M
C( I, J ) = BETA*C( I, J )
30 CONTINUE
40 CONTINUE
END IF
RETURN
END IF
*
* Start the operations.
*
IF( LSAME( SIDE, 'L' ) )THEN
*
* Form C := alpha*A*B + beta*C.
*
IF( UPPER )THEN
DO 70, J = 1, N
DO 60, I = 1, M
TEMP1 = ALPHA*B( I, J )
TEMP2 = ZERO
DO 50, K = 1, I - 1
C( K, J ) = C( K, J ) + TEMP1*A( K, I )
TEMP2 = TEMP2 +
$ B( K, J )*DCONJG( A( K, I ) )
50 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = TEMP1*DBLE( A( I, I ) ) +
$ ALPHA*TEMP2
ELSE
C( I, J ) = BETA *C( I, J ) +
$ TEMP1*DBLE( A( I, I ) ) +
$ ALPHA*TEMP2
END IF
60 CONTINUE
70 CONTINUE
ELSE
DO 100, J = 1, N
DO 90, I = M, 1, -1
TEMP1 = ALPHA*B( I, J )
TEMP2 = ZERO
DO 80, K = I + 1, M
C( K, J ) = C( K, J ) + TEMP1*A( K, I )
TEMP2 = TEMP2 +
$ B( K, J )*DCONJG( A( K, I ) )
80 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = TEMP1*DBLE( A( I, I ) ) +
$ ALPHA*TEMP2
ELSE
C( I, J ) = BETA *C( I, J ) +
$ TEMP1*DBLE( A( I, I ) ) +
$ ALPHA*TEMP2
END IF
90 CONTINUE
100 CONTINUE
END IF
ELSE
*
* Form C := alpha*B*A + beta*C.
*
DO 170, J = 1, N
TEMP1 = ALPHA*DBLE( A( J, J ) )
IF( BETA.EQ.ZERO )THEN
DO 110, I = 1, M
C( I, J ) = TEMP1*B( I, J )
110 CONTINUE
ELSE
DO 120, I = 1, M
C( I, J ) = BETA*C( I, J ) + TEMP1*B( I, J )
120 CONTINUE
END IF
DO 140, K = 1, J - 1
IF( UPPER )THEN
TEMP1 = ALPHA*A( K, J )
ELSE
TEMP1 = ALPHA*DCONJG( A( J, K ) )
END IF
DO 130, I = 1, M
C( I, J ) = C( I, J ) + TEMP1*B( I, K )
130 CONTINUE
140 CONTINUE
DO 160, K = J + 1, N
IF( UPPER )THEN
TEMP1 = ALPHA*DCONJG( A( J, K ) )
ELSE
TEMP1 = ALPHA*A( K, J )
END IF
DO 150, I = 1, M
C( I, J ) = C( I, J ) + TEMP1*B( I, K )
150 CONTINUE
160 CONTINUE
170 CONTINUE
END IF
*
RETURN
*
* End of ZHEMM .
*
END
SUBROUTINE ZHER ( UPLO, N, ALPHA, X, INCX, A, LDA )
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA
INTEGER INCX, LDA, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * )
* ..
*
* Purpose
* =======
*
* ZHER performs the hermitian rank 1 operation
*
* A := alpha*x*conjg( x' ) + A,
*
* where alpha is a real scalar, x is an n element vector and A is an
* n by n hermitian matrix.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the array A is to be referenced as
* follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of A
* is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of A
* is to be referenced.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - DOUBLE PRECISION.
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array A must contain the upper
* triangular part of the hermitian matrix and the strictly
* lower triangular part of A is not referenced. On exit, the
* upper triangular part of the array A is overwritten by the
* upper triangular part of the updated matrix.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array A must contain the lower
* triangular part of the hermitian matrix and the strictly
* upper triangular part of A is not referenced. On exit, the
* lower triangular part of the array A is overwritten by the
* lower triangular part of the updated matrix.
* Note that the imaginary parts of the diagonal elements need
* not be set, they are assumed to be zero, and on exit they
* are set to zero.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, n ).
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, J, JX, KX
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX, DBLE
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
$ .NOT.LSAME( UPLO, 'L' ) )THEN
INFO = 1
ELSE IF( N.LT.0 )THEN
INFO = 2
ELSE IF( INCX.EQ.0 )THEN
INFO = 5
ELSE IF( LDA.LT.MAX( 1, N ) )THEN
INFO = 7
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZHER ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ALPHA.EQ.DBLE( ZERO ) ) )
$ RETURN
*
* Set the start point in X if the increment is not unity.
*
IF( INCX.LE.0 )THEN
KX = 1 - ( N - 1 )*INCX
ELSE IF( INCX.NE.1 )THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through the triangular part
* of A.
*
IF( LSAME( UPLO, 'U' ) )THEN
*
* Form A when A is stored in upper triangle.
*
IF( INCX.EQ.1 )THEN
DO 20, J = 1, N
IF( X( J ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( X( J ) )
DO 10, I = 1, J - 1
A( I, J ) = A( I, J ) + X( I )*TEMP
10 CONTINUE
A( J, J ) = DBLE( A( J, J ) ) + DBLE( X( J )*TEMP )
ELSE
A( J, J ) = DBLE( A( J, J ) )
END IF
20 CONTINUE
ELSE
JX = KX
DO 40, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( X( JX ) )
IX = KX
DO 30, I = 1, J - 1
A( I, J ) = A( I, J ) + X( IX )*TEMP
IX = IX + INCX
30 CONTINUE
A( J, J ) = DBLE( A( J, J ) ) + DBLE( X( JX )*TEMP )
ELSE
A( J, J ) = DBLE( A( J, J ) )
END IF
JX = JX + INCX
40 CONTINUE
END IF
ELSE
*
* Form A when A is stored in lower triangle.
*
IF( INCX.EQ.1 )THEN
DO 60, J = 1, N
IF( X( J ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( X( J ) )
A( J, J ) = DBLE( A( J, J ) ) + DBLE( TEMP*X( J ) )
DO 50, I = J + 1, N
A( I, J ) = A( I, J ) + X( I )*TEMP
50 CONTINUE
ELSE
A( J, J ) = DBLE( A( J, J ) )
END IF
60 CONTINUE
ELSE
JX = KX
DO 80, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( X( JX ) )
A( J, J ) = DBLE( A( J, J ) ) + DBLE( TEMP*X( JX ) )
IX = JX
DO 70, I = J + 1, N
IX = IX + INCX
A( I, J ) = A( I, J ) + X( IX )*TEMP
70 CONTINUE
ELSE
A( J, J ) = DBLE( A( J, J ) )
END IF
JX = JX + INCX
80 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZHER .
*
END
SUBROUTINE ZHERK( UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC )
* .. Scalar Arguments ..
CHARACTER TRANS, UPLO
INTEGER K, LDA, LDC, N
DOUBLE PRECISION ALPHA, BETA
* ..
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), C( LDC, * )
* ..
*
* Purpose
* =======
*
* ZHERK performs one of the hermitian rank k operations
*
* C := alpha*A*conjg( A' ) + beta*C,
*
* or
*
* C := alpha*conjg( A' )*A + beta*C,
*
* where alpha and beta are real scalars, C is an n by n hermitian
* matrix and A is an n by k matrix in the first case and a k by n
* matrix in the second case.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the array C is to be referenced as
* follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of C
* is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of C
* is to be referenced.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' C := alpha*A*conjg( A' ) + beta*C.
*
* TRANS = 'C' or 'c' C := alpha*conjg( A' )*A + beta*C.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix C. N must be
* at least zero.
* Unchanged on exit.
*
* K - INTEGER.
* On entry with TRANS = 'N' or 'n', K specifies the number
* of columns of the matrix A, and on entry with
* TRANS = 'C' or 'c', K specifies the number of rows of the
* matrix A. K must be at least zero.
* Unchanged on exit.
*
* ALPHA - DOUBLE PRECISION .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
* k when TRANS = 'N' or 'n', and is n otherwise.
* Before entry with TRANS = 'N' or 'n', the leading n by k
* part of the array A must contain the matrix A, otherwise
* the leading k by n part of the array A must contain the
* matrix A.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When TRANS = 'N' or 'n'
* then LDA must be at least max( 1, n ), otherwise LDA must
* be at least max( 1, k ).
* Unchanged on exit.
*
* BETA - DOUBLE PRECISION.
* On entry, BETA specifies the scalar beta.
* Unchanged on exit.
*
* C - COMPLEX*16 array of DIMENSION ( LDC, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array C must contain the upper
* triangular part of the hermitian matrix and the strictly
* lower triangular part of C is not referenced. On exit, the
* upper triangular part of the array C is overwritten by the
* upper triangular part of the updated matrix.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array C must contain the lower
* triangular part of the hermitian matrix and the strictly
* upper triangular part of C is not referenced. On exit, the
* lower triangular part of the array C is overwritten by the
* lower triangular part of the updated matrix.
* Note that the imaginary parts of the diagonal elements need
* not be set, they are assumed to be zero, and on exit they
* are set to zero.
*
* LDC - INTEGER.
* On entry, LDC specifies the first dimension of C as declared
* in the calling (sub) program. LDC must be at least
* max( 1, n ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
* -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1.
* Ed Anderson, Cray Research Inc.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC DBLE, DCMPLX, DCONJG, MAX
* ..
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I, INFO, J, L, NROWA
DOUBLE PRECISION RTEMP
COMPLEX*16 TEMP
* ..
* .. Parameters ..
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
IF( LSAME( TRANS, 'N' ) ) THEN
NROWA = N
ELSE
NROWA = K
END IF
UPPER = LSAME( UPLO, 'U' )
*
INFO = 0
IF( ( .NOT.UPPER ) .AND. ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
INFO = 1
ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ) .AND.
$ ( .NOT.LSAME( TRANS, 'C' ) ) ) THEN
INFO = 2
ELSE IF( N.LT.0 ) THEN
INFO = 3
ELSE IF( K.LT.0 ) THEN
INFO = 4
ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN
INFO = 7
ELSE IF( LDC.LT.MAX( 1, N ) ) THEN
INFO = 10
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZHERK ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND.
$ ( BETA.EQ.ONE ) ) )RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO ) THEN
IF( UPPER ) THEN
IF( BETA.EQ.ZERO ) THEN
DO 20 J = 1, N
DO 10 I = 1, J
C( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40 J = 1, N
DO 30 I = 1, J - 1
C( I, J ) = BETA*C( I, J )
30 CONTINUE
C( J, J ) = BETA*DBLE( C( J, J ) )
40 CONTINUE
END IF
ELSE
IF( BETA.EQ.ZERO ) THEN
DO 60 J = 1, N
DO 50 I = J, N
C( I, J ) = ZERO
50 CONTINUE
60 CONTINUE
ELSE
DO 80 J = 1, N
C( J, J ) = BETA*DBLE( C( J, J ) )
DO 70 I = J + 1, N
C( I, J ) = BETA*C( I, J )
70 CONTINUE
80 CONTINUE
END IF
END IF
RETURN
END IF
*
* Start the operations.
*
IF( LSAME( TRANS, 'N' ) ) THEN
*
* Form C := alpha*A*conjg( A' ) + beta*C.
*
IF( UPPER ) THEN
DO 130 J = 1, N
IF( BETA.EQ.ZERO ) THEN
DO 90 I = 1, J
C( I, J ) = ZERO
90 CONTINUE
ELSE IF( BETA.NE.ONE ) THEN
DO 100 I = 1, J - 1
C( I, J ) = BETA*C( I, J )
100 CONTINUE
C( J, J ) = BETA*DBLE( C( J, J ) )
ELSE
C( J, J ) = DBLE( C( J, J ) )
END IF
DO 120 L = 1, K
IF( A( J, L ).NE.DCMPLX( ZERO ) ) THEN
TEMP = ALPHA*DCONJG( A( J, L ) )
DO 110 I = 1, J - 1
C( I, J ) = C( I, J ) + TEMP*A( I, L )
110 CONTINUE
C( J, J ) = DBLE( C( J, J ) ) +
$ DBLE( TEMP*A( I, L ) )
END IF
120 CONTINUE
130 CONTINUE
ELSE
DO 180 J = 1, N
IF( BETA.EQ.ZERO ) THEN
DO 140 I = J, N
C( I, J ) = ZERO
140 CONTINUE
ELSE IF( BETA.NE.ONE ) THEN
C( J, J ) = BETA*DBLE( C( J, J ) )
DO 150 I = J + 1, N
C( I, J ) = BETA*C( I, J )
150 CONTINUE
ELSE
C( J, J ) = DBLE( C( J, J ) )
END IF
DO 170 L = 1, K
IF( A( J, L ).NE.DCMPLX( ZERO ) ) THEN
TEMP = ALPHA*DCONJG( A( J, L ) )
C( J, J ) = DBLE( C( J, J ) ) +
$ DBLE( TEMP*A( J, L ) )
DO 160 I = J + 1, N
C( I, J ) = C( I, J ) + TEMP*A( I, L )
160 CONTINUE
END IF
170 CONTINUE
180 CONTINUE
END IF
ELSE
*
* Form C := alpha*conjg( A' )*A + beta*C.
*
IF( UPPER ) THEN
DO 220 J = 1, N
DO 200 I = 1, J - 1
TEMP = ZERO
DO 190 L = 1, K
TEMP = TEMP + DCONJG( A( L, I ) )*A( L, J )
190 CONTINUE
IF( BETA.EQ.ZERO ) THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
200 CONTINUE
RTEMP = ZERO
DO 210 L = 1, K
RTEMP = RTEMP + DCONJG( A( L, J ) )*A( L, J )
210 CONTINUE
IF( BETA.EQ.ZERO ) THEN
C( J, J ) = ALPHA*RTEMP
ELSE
C( J, J ) = ALPHA*RTEMP + BETA*DBLE( C( J, J ) )
END IF
220 CONTINUE
ELSE
DO 260 J = 1, N
RTEMP = ZERO
DO 230 L = 1, K
RTEMP = RTEMP + DCONJG( A( L, J ) )*A( L, J )
230 CONTINUE
IF( BETA.EQ.ZERO ) THEN
C( J, J ) = ALPHA*RTEMP
ELSE
C( J, J ) = ALPHA*RTEMP + BETA*DBLE( C( J, J ) )
END IF
DO 250 I = J + 1, N
TEMP = ZERO
DO 240 L = 1, K
TEMP = TEMP + DCONJG( A( L, I ) )*A( L, J )
240 CONTINUE
IF( BETA.EQ.ZERO ) THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
250 CONTINUE
260 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZHERK .
*
END
SUBROUTINE ZHPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
* .. Scalar Arguments ..
COMPLEX*16 ALPHA, BETA
INTEGER INCX, INCY, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 AP( * ), X( * ), Y( * )
* ..
*
* Purpose
* =======
*
* ZHPMV performs the matrix-vector operation
*
* y := alpha*A*x + beta*y,
*
* where alpha and beta are scalars, x and y are n element vectors and
* A is an n by n hermitian matrix, supplied in packed form.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the matrix A is supplied in the packed
* array AP as follows:
*
* UPLO = 'U' or 'u' The upper triangular part of A is
* supplied in AP.
*
* UPLO = 'L' or 'l' The lower triangular part of A is
* supplied in AP.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* AP - COMPLEX*16 array of DIMENSION at least
* ( ( n*( n + 1 ) )/2 ).
* Before entry with UPLO = 'U' or 'u', the array AP must
* contain the upper triangular part of the hermitian matrix
* packed sequentially, column by column, so that AP( 1 )
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
* and a( 2, 2 ) respectively, and so on.
* Before entry with UPLO = 'L' or 'l', the array AP must
* contain the lower triangular part of the hermitian matrix
* packed sequentially, column by column, so that AP( 1 )
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
* and a( 3, 1 ) respectively, and so on.
* Note that the imaginary parts of the diagonal elements need
* not be set and are assumed to be zero.
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* BETA - COMPLEX*16 .
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then Y need not be set on input.
* Unchanged on exit.
*
* Y - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCY ) ).
* Before entry, the incremented array Y must contain the n
* element vector y. On exit, Y is overwritten by the updated
* vector y.
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP1, TEMP2
INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, DBLE
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
$ .NOT.LSAME( UPLO, 'L' ) )THEN
INFO = 1
ELSE IF( N.LT.0 )THEN
INFO = 2
ELSE IF( INCX.EQ.0 )THEN
INFO = 6
ELSE IF( INCY.EQ.0 )THEN
INFO = 9
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZHPMV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* Set up the start points in X and Y.
*
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( N - 1 )*INCX
END IF
IF( INCY.GT.0 )THEN
KY = 1
ELSE
KY = 1 - ( N - 1 )*INCY
END IF
*
* Start the operations. In this version the elements of the array AP
* are accessed sequentially with one pass through AP.
*
* First form y := beta*y.
*
IF( BETA.NE.ONE )THEN
IF( INCY.EQ.1 )THEN
IF( BETA.EQ.ZERO )THEN
DO 10, I = 1, N
Y( I ) = ZERO
10 CONTINUE
ELSE
DO 20, I = 1, N
Y( I ) = BETA*Y( I )
20 CONTINUE
END IF
ELSE
IY = KY
IF( BETA.EQ.ZERO )THEN
DO 30, I = 1, N
Y( IY ) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40, I = 1, N
Y( IY ) = BETA*Y( IY )
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF( ALPHA.EQ.ZERO )
$ RETURN
KK = 1
IF( LSAME( UPLO, 'U' ) )THEN
*
* Form y when AP contains the upper triangle.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 60, J = 1, N
TEMP1 = ALPHA*X( J )
TEMP2 = ZERO
K = KK
DO 50, I = 1, J - 1
Y( I ) = Y( I ) + TEMP1*AP( K )
TEMP2 = TEMP2 + DCONJG( AP( K ) )*X( I )
K = K + 1
50 CONTINUE
Y( J ) = Y( J ) + TEMP1*DBLE( AP( KK + J - 1 ) )
$ + ALPHA*TEMP2
KK = KK + J
60 CONTINUE
ELSE
JX = KX
JY = KY
DO 80, J = 1, N
TEMP1 = ALPHA*X( JX )
TEMP2 = ZERO
IX = KX
IY = KY
DO 70, K = KK, KK + J - 2
Y( IY ) = Y( IY ) + TEMP1*AP( K )
TEMP2 = TEMP2 + DCONJG( AP( K ) )*X( IX )
IX = IX + INCX
IY = IY + INCY
70 CONTINUE
Y( JY ) = Y( JY ) + TEMP1*DBLE( AP( KK + J - 1 ) )
$ + ALPHA*TEMP2
JX = JX + INCX
JY = JY + INCY
KK = KK + J
80 CONTINUE
END IF
ELSE
*
* Form y when AP contains the lower triangle.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 100, J = 1, N
TEMP1 = ALPHA*X( J )
TEMP2 = ZERO
Y( J ) = Y( J ) + TEMP1*DBLE( AP( KK ) )
K = KK + 1
DO 90, I = J + 1, N
Y( I ) = Y( I ) + TEMP1*AP( K )
TEMP2 = TEMP2 + DCONJG( AP( K ) )*X( I )
K = K + 1
90 CONTINUE
Y( J ) = Y( J ) + ALPHA*TEMP2
KK = KK + ( N - J + 1 )
100 CONTINUE
ELSE
JX = KX
JY = KY
DO 120, J = 1, N
TEMP1 = ALPHA*X( JX )
TEMP2 = ZERO
Y( JY ) = Y( JY ) + TEMP1*DBLE( AP( KK ) )
IX = JX
IY = JY
DO 110, K = KK + 1, KK + N - J
IX = IX + INCX
IY = IY + INCY
Y( IY ) = Y( IY ) + TEMP1*AP( K )
TEMP2 = TEMP2 + DCONJG( AP( K ) )*X( IX )
110 CONTINUE
Y( JY ) = Y( JY ) + ALPHA*TEMP2
JX = JX + INCX
JY = JY + INCY
KK = KK + ( N - J + 1 )
120 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZHPMV .
*
END
SUBROUTINE ZHPR ( UPLO, N, ALPHA, X, INCX, AP )
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA
INTEGER INCX, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 AP( * ), X( * )
* ..
*
* Purpose
* =======
*
* ZHPR performs the hermitian rank 1 operation
*
* A := alpha*x*conjg( x' ) + A,
*
* where alpha is a real scalar, x is an n element vector and A is an
* n by n hermitian matrix, supplied in packed form.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the matrix A is supplied in the packed
* array AP as follows:
*
* UPLO = 'U' or 'u' The upper triangular part of A is
* supplied in AP.
*
* UPLO = 'L' or 'l' The lower triangular part of A is
* supplied in AP.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - DOUBLE PRECISION.
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* AP - COMPLEX*16 array of DIMENSION at least
* ( ( n*( n + 1 ) )/2 ).
* Before entry with UPLO = 'U' or 'u', the array AP must
* contain the upper triangular part of the hermitian matrix
* packed sequentially, column by column, so that AP( 1 )
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
* and a( 2, 2 ) respectively, and so on. On exit, the array
* AP is overwritten by the upper triangular part of the
* updated matrix.
* Before entry with UPLO = 'L' or 'l', the array AP must
* contain the lower triangular part of the hermitian matrix
* packed sequentially, column by column, so that AP( 1 )
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
* and a( 3, 1 ) respectively, and so on. On exit, the array
* AP is overwritten by the lower triangular part of the
* updated matrix.
* Note that the imaginary parts of the diagonal elements need
* not be set, they are assumed to be zero, and on exit they
* are set to zero.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, J, JX, K, KK, KX
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, DBLE
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
$ .NOT.LSAME( UPLO, 'L' ) )THEN
INFO = 1
ELSE IF( N.LT.0 )THEN
INFO = 2
ELSE IF( INCX.EQ.0 )THEN
INFO = 5
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZHPR ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ALPHA.EQ.DBLE( ZERO ) ) )
$ RETURN
*
* Set the start point in X if the increment is not unity.
*
IF( INCX.LE.0 )THEN
KX = 1 - ( N - 1 )*INCX
ELSE IF( INCX.NE.1 )THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of the array AP
* are accessed sequentially with one pass through AP.
*
KK = 1
IF( LSAME( UPLO, 'U' ) )THEN
*
* Form A when upper triangle is stored in AP.
*
IF( INCX.EQ.1 )THEN
DO 20, J = 1, N
IF( X( J ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( X( J ) )
K = KK
DO 10, I = 1, J - 1
AP( K ) = AP( K ) + X( I )*TEMP
K = K + 1
10 CONTINUE
AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
$ + DBLE( X( J )*TEMP )
ELSE
AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
END IF
KK = KK + J
20 CONTINUE
ELSE
JX = KX
DO 40, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( X( JX ) )
IX = KX
DO 30, K = KK, KK + J - 2
AP( K ) = AP( K ) + X( IX )*TEMP
IX = IX + INCX
30 CONTINUE
AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
$ + DBLE( X( JX )*TEMP )
ELSE
AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
END IF
JX = JX + INCX
KK = KK + J
40 CONTINUE
END IF
ELSE
*
* Form A when lower triangle is stored in AP.
*
IF( INCX.EQ.1 )THEN
DO 60, J = 1, N
IF( X( J ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( X( J ) )
AP( KK ) = DBLE( AP( KK ) ) + DBLE( TEMP*X( J ) )
K = KK + 1
DO 50, I = J + 1, N
AP( K ) = AP( K ) + X( I )*TEMP
K = K + 1
50 CONTINUE
ELSE
AP( KK ) = DBLE( AP( KK ) )
END IF
KK = KK + N - J + 1
60 CONTINUE
ELSE
JX = KX
DO 80, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( X( JX ) )
AP( KK ) = DBLE( AP( KK ) ) + DBLE( TEMP*X( JX ) )
IX = JX
DO 70, K = KK + 1, KK + N - J
IX = IX + INCX
AP( K ) = AP( K ) + X( IX )*TEMP
70 CONTINUE
ELSE
AP( KK ) = DBLE( AP( KK ) )
END IF
JX = JX + INCX
KK = KK + N - J + 1
80 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZHPR .
*
END
SUBROUTINE ZHPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP )
* .. Scalar Arguments ..
COMPLEX*16 ALPHA
INTEGER INCX, INCY, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 AP( * ), X( * ), Y( * )
* ..
*
* Purpose
* =======
*
* ZHPR2 performs the hermitian rank 2 operation
*
* A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
*
* where alpha is a scalar, x and y are n element vectors and A is an
* n by n hermitian matrix, supplied in packed form.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the matrix A is supplied in the packed
* array AP as follows:
*
* UPLO = 'U' or 'u' The upper triangular part of A is
* supplied in AP.
*
* UPLO = 'L' or 'l' The lower triangular part of A is
* supplied in AP.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* Y - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCY ) ).
* Before entry, the incremented array Y must contain the n
* element vector y.
* Unchanged on exit.
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
* AP - COMPLEX*16 array of DIMENSION at least
* ( ( n*( n + 1 ) )/2 ).
* Before entry with UPLO = 'U' or 'u', the array AP must
* contain the upper triangular part of the hermitian matrix
* packed sequentially, column by column, so that AP( 1 )
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
* and a( 2, 2 ) respectively, and so on. On exit, the array
* AP is overwritten by the upper triangular part of the
* updated matrix.
* Before entry with UPLO = 'L' or 'l', the array AP must
* contain the lower triangular part of the hermitian matrix
* packed sequentially, column by column, so that AP( 1 )
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
* and a( 3, 1 ) respectively, and so on. On exit, the array
* AP is overwritten by the lower triangular part of the
* updated matrix.
* Note that the imaginary parts of the diagonal elements need
* not be set, they are assumed to be zero, and on exit they
* are set to zero.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP1, TEMP2
INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, DBLE
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
$ .NOT.LSAME( UPLO, 'L' ) )THEN
INFO = 1
ELSE IF( N.LT.0 )THEN
INFO = 2
ELSE IF( INCX.EQ.0 )THEN
INFO = 5
ELSE IF( INCY.EQ.0 )THEN
INFO = 7
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZHPR2 ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
$ RETURN
*
* Set up the start points in X and Y if the increments are not both
* unity.
*
IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( N - 1 )*INCX
END IF
IF( INCY.GT.0 )THEN
KY = 1
ELSE
KY = 1 - ( N - 1 )*INCY
END IF
JX = KX
JY = KY
END IF
*
* Start the operations. In this version the elements of the array AP
* are accessed sequentially with one pass through AP.
*
KK = 1
IF( LSAME( UPLO, 'U' ) )THEN
*
* Form A when upper triangle is stored in AP.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 20, J = 1, N
IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
TEMP1 = ALPHA*DCONJG( Y( J ) )
TEMP2 = DCONJG( ALPHA*X( J ) )
K = KK
DO 10, I = 1, J - 1
AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2
K = K + 1
10 CONTINUE
AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) +
$ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
ELSE
AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
END IF
KK = KK + J
20 CONTINUE
ELSE
DO 40, J = 1, N
IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
TEMP1 = ALPHA*DCONJG( Y( JY ) )
TEMP2 = DCONJG( ALPHA*X( JX ) )
IX = KX
IY = KY
DO 30, K = KK, KK + J - 2
AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2
IX = IX + INCX
IY = IY + INCY
30 CONTINUE
AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) +
$ DBLE( X( JX )*TEMP1 +
$ Y( JY )*TEMP2 )
ELSE
AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
END IF
JX = JX + INCX
JY = JY + INCY
KK = KK + J
40 CONTINUE
END IF
ELSE
*
* Form A when lower triangle is stored in AP.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 60, J = 1, N
IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
TEMP1 = ALPHA*DCONJG( Y( J ) )
TEMP2 = DCONJG( ALPHA*X( J ) )
AP( KK ) = DBLE( AP( KK ) ) +
$ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
K = KK + 1
DO 50, I = J + 1, N
AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2
K = K + 1
50 CONTINUE
ELSE
AP( KK ) = DBLE( AP( KK ) )
END IF
KK = KK + N - J + 1
60 CONTINUE
ELSE
DO 80, J = 1, N
IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
TEMP1 = ALPHA*DCONJG( Y( JY ) )
TEMP2 = DCONJG( ALPHA*X( JX ) )
AP( KK ) = DBLE( AP( KK ) ) +
$ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
IX = JX
IY = JY
DO 70, K = KK + 1, KK + N - J
IX = IX + INCX
IY = IY + INCY
AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2
70 CONTINUE
ELSE
AP( KK ) = DBLE( AP( KK ) )
END IF
JX = JX + INCX
JY = JY + INCY
KK = KK + N - J + 1
80 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZHPR2 .
*
END
subroutine zrotg(ca,cb,c,s)
double complex ca,cb,s
double precision c
double precision norm,scale
double complex alpha
intrinsic dconjg, dcmplx
if (cdabs(ca) .ne. 0.0d0) go to 10
c = 0.0d0
s = (1.0d0,0.0d0)
ca = cb
go to 20
10 continue
scale = cdabs(ca) + cdabs(cb)
norm = scale*dsqrt((cdabs(ca/dcmplx(scale,0.0d0)))**2 +
* (cdabs(cb/dcmplx(scale,0.0d0)))**2)
alpha = ca /cdabs(ca)
c = cdabs(ca) / norm
s = alpha * dconjg(cb) / norm
ca = alpha * norm
20 continue
return
end
SUBROUTINE ZSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB,
$ BETA, C, LDC )
* .. Scalar Arguments ..
CHARACTER*1 SIDE, UPLO
INTEGER M, N, LDA, LDB, LDC
COMPLEX*16 ALPHA, BETA
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
* ..
*
* Purpose
* =======
*
* ZSYMM performs one of the matrix-matrix operations
*
* C := alpha*A*B + beta*C,
*
* or
*
* C := alpha*B*A + beta*C,
*
* where alpha and beta are scalars, A is a symmetric matrix and B and
* C are m by n matrices.
*
* Parameters
* ==========
*
* SIDE - CHARACTER*1.
* On entry, SIDE specifies whether the symmetric matrix A
* appears on the left or right in the operation as follows:
*
* SIDE = 'L' or 'l' C := alpha*A*B + beta*C,
*
* SIDE = 'R' or 'r' C := alpha*B*A + beta*C,
*
* Unchanged on exit.
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the symmetric matrix A is to be
* referenced as follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of the
* symmetric matrix is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of the
* symmetric matrix is to be referenced.
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of the matrix C.
* M must be at least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of the matrix C.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
* m when SIDE = 'L' or 'l' and is n otherwise.
* Before entry with SIDE = 'L' or 'l', the m by m part of
* the array A must contain the symmetric matrix, such that
* when UPLO = 'U' or 'u', the leading m by m upper triangular
* part of the array A must contain the upper triangular part
* of the symmetric matrix and the strictly lower triangular
* part of A is not referenced, and when UPLO = 'L' or 'l',
* the leading m by m lower triangular part of the array A
* must contain the lower triangular part of the symmetric
* matrix and the strictly upper triangular part of A is not
* referenced.
* Before entry with SIDE = 'R' or 'r', the n by n part of
* the array A must contain the symmetric matrix, such that
* when UPLO = 'U' or 'u', the leading n by n upper triangular
* part of the array A must contain the upper triangular part
* of the symmetric matrix and the strictly lower triangular
* part of A is not referenced, and when UPLO = 'L' or 'l',
* the leading n by n lower triangular part of the array A
* must contain the lower triangular part of the symmetric
* matrix and the strictly upper triangular part of A is not
* referenced.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When SIDE = 'L' or 'l' then
* LDA must be at least max( 1, m ), otherwise LDA must be at
* least max( 1, n ).
* Unchanged on exit.
*
* B - COMPLEX*16 array of DIMENSION ( LDB, n ).
* Before entry, the leading m by n part of the array B must
* contain the matrix B.
* Unchanged on exit.
*
* LDB - INTEGER.
* On entry, LDB specifies the first dimension of B as declared
* in the calling (sub) program. LDB must be at least
* max( 1, m ).
* Unchanged on exit.
*
* BETA - COMPLEX*16 .
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then C need not be set on input.
* Unchanged on exit.
*
* C - COMPLEX*16 array of DIMENSION ( LDC, n ).
* Before entry, the leading m by n part of the array C must
* contain the matrix C, except when beta is zero, in which
* case C need not be set on entry.
* On exit, the array C is overwritten by the m by n updated
* matrix.
*
* LDC - INTEGER.
* On entry, LDC specifies the first dimension of C as declared
* in the calling (sub) program. LDC must be at least
* max( 1, m ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I, INFO, J, K, NROWA
COMPLEX*16 TEMP1, TEMP2
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* ..
* .. Executable Statements ..
*
* Set NROWA as the number of rows of A.
*
IF( LSAME( SIDE, 'L' ) )THEN
NROWA = M
ELSE
NROWA = N
END IF
UPPER = LSAME( UPLO, 'U' )
*
* Test the input parameters.
*
INFO = 0
IF( ( .NOT.LSAME( SIDE, 'L' ) ).AND.
$ ( .NOT.LSAME( SIDE, 'R' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.UPPER ).AND.
$ ( .NOT.LSAME( UPLO, 'L' ) ) )THEN
INFO = 2
ELSE IF( M .LT.0 )THEN
INFO = 3
ELSE IF( N .LT.0 )THEN
INFO = 4
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 7
ELSE IF( LDB.LT.MAX( 1, M ) )THEN
INFO = 9
ELSE IF( LDC.LT.MAX( 1, M ) )THEN
INFO = 12
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZSYMM ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
$ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
IF( BETA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, M
C( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40, J = 1, N
DO 30, I = 1, M
C( I, J ) = BETA*C( I, J )
30 CONTINUE
40 CONTINUE
END IF
RETURN
END IF
*
* Start the operations.
*
IF( LSAME( SIDE, 'L' ) )THEN
*
* Form C := alpha*A*B + beta*C.
*
IF( UPPER )THEN
DO 70, J = 1, N
DO 60, I = 1, M
TEMP1 = ALPHA*B( I, J )
TEMP2 = ZERO
DO 50, K = 1, I - 1
C( K, J ) = C( K, J ) + TEMP1 *A( K, I )
TEMP2 = TEMP2 + B( K, J )*A( K, I )
50 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2
ELSE
C( I, J ) = BETA *C( I, J ) +
$ TEMP1*A( I, I ) + ALPHA*TEMP2
END IF
60 CONTINUE
70 CONTINUE
ELSE
DO 100, J = 1, N
DO 90, I = M, 1, -1
TEMP1 = ALPHA*B( I, J )
TEMP2 = ZERO
DO 80, K = I + 1, M
C( K, J ) = C( K, J ) + TEMP1 *A( K, I )
TEMP2 = TEMP2 + B( K, J )*A( K, I )
80 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2
ELSE
C( I, J ) = BETA *C( I, J ) +
$ TEMP1*A( I, I ) + ALPHA*TEMP2
END IF
90 CONTINUE
100 CONTINUE
END IF
ELSE
*
* Form C := alpha*B*A + beta*C.
*
DO 170, J = 1, N
TEMP1 = ALPHA*A( J, J )
IF( BETA.EQ.ZERO )THEN
DO 110, I = 1, M
C( I, J ) = TEMP1*B( I, J )
110 CONTINUE
ELSE
DO 120, I = 1, M
C( I, J ) = BETA*C( I, J ) + TEMP1*B( I, J )
120 CONTINUE
END IF
DO 140, K = 1, J - 1
IF( UPPER )THEN
TEMP1 = ALPHA*A( K, J )
ELSE
TEMP1 = ALPHA*A( J, K )
END IF
DO 130, I = 1, M
C( I, J ) = C( I, J ) + TEMP1*B( I, K )
130 CONTINUE
140 CONTINUE
DO 160, K = J + 1, N
IF( UPPER )THEN
TEMP1 = ALPHA*A( J, K )
ELSE
TEMP1 = ALPHA*A( K, J )
END IF
DO 150, I = 1, M
C( I, J ) = C( I, J ) + TEMP1*B( I, K )
150 CONTINUE
160 CONTINUE
170 CONTINUE
END IF
*
RETURN
*
* End of ZSYMM .
*
END
SUBROUTINE ZSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB,
$ BETA, C, LDC )
* .. Scalar Arguments ..
CHARACTER*1 UPLO, TRANS
INTEGER N, K, LDA, LDB, LDC
COMPLEX*16 ALPHA, BETA
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
* ..
*
* Purpose
* =======
*
* ZSYR2K performs one of the symmetric rank 2k operations
*
* C := alpha*A*B' + alpha*B*A' + beta*C,
*
* or
*
* C := alpha*A'*B + alpha*B'*A + beta*C,
*
* where alpha and beta are scalars, C is an n by n symmetric matrix
* and A and B are n by k matrices in the first case and k by n
* matrices in the second case.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the array C is to be referenced as
* follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of C
* is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of C
* is to be referenced.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' C := alpha*A*B' + alpha*B*A' +
* beta*C.
*
* TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A +
* beta*C.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix C. N must be
* at least zero.
* Unchanged on exit.
*
* K - INTEGER.
* On entry with TRANS = 'N' or 'n', K specifies the number
* of columns of the matrices A and B, and on entry with
* TRANS = 'T' or 't', K specifies the number of rows of the
* matrices A and B. K must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
* k when TRANS = 'N' or 'n', and is n otherwise.
* Before entry with TRANS = 'N' or 'n', the leading n by k
* part of the array A must contain the matrix A, otherwise
* the leading k by n part of the array A must contain the
* matrix A.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When TRANS = 'N' or 'n'
* then LDA must be at least max( 1, n ), otherwise LDA must
* be at least max( 1, k ).
* Unchanged on exit.
*
* B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is
* k when TRANS = 'N' or 'n', and is n otherwise.
* Before entry with TRANS = 'N' or 'n', the leading n by k
* part of the array B must contain the matrix B, otherwise
* the leading k by n part of the array B must contain the
* matrix B.
* Unchanged on exit.
*
* LDB - INTEGER.
* On entry, LDB specifies the first dimension of B as declared
* in the calling (sub) program. When TRANS = 'N' or 'n'
* then LDB must be at least max( 1, n ), otherwise LDB must
* be at least max( 1, k ).
* Unchanged on exit.
*
* BETA - COMPLEX*16 .
* On entry, BETA specifies the scalar beta.
* Unchanged on exit.
*
* C - COMPLEX*16 array of DIMENSION ( LDC, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array C must contain the upper
* triangular part of the symmetric matrix and the strictly
* lower triangular part of C is not referenced. On exit, the
* upper triangular part of the array C is overwritten by the
* upper triangular part of the updated matrix.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array C must contain the lower
* triangular part of the symmetric matrix and the strictly
* upper triangular part of C is not referenced. On exit, the
* lower triangular part of the array C is overwritten by the
* lower triangular part of the updated matrix.
*
* LDC - INTEGER.
* On entry, LDC specifies the first dimension of C as declared
* in the calling (sub) program. LDC must be at least
* max( 1, n ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I, INFO, J, L, NROWA
COMPLEX*16 TEMP1, TEMP2
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
IF( LSAME( TRANS, 'N' ) )THEN
NROWA = N
ELSE
NROWA = K
END IF
UPPER = LSAME( UPLO, 'U' )
*
INFO = 0
IF( ( .NOT.UPPER ).AND.
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
$ ( .NOT.LSAME( TRANS, 'T' ) ) )THEN
INFO = 2
ELSE IF( N .LT.0 )THEN
INFO = 3
ELSE IF( K .LT.0 )THEN
INFO = 4
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 7
ELSE IF( LDB.LT.MAX( 1, NROWA ) )THEN
INFO = 9
ELSE IF( LDC.LT.MAX( 1, N ) )THEN
INFO = 12
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZSYR2K', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.
$ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
IF( UPPER )THEN
IF( BETA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, J
C( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40, J = 1, N
DO 30, I = 1, J
C( I, J ) = BETA*C( I, J )
30 CONTINUE
40 CONTINUE
END IF
ELSE
IF( BETA.EQ.ZERO )THEN
DO 60, J = 1, N
DO 50, I = J, N
C( I, J ) = ZERO
50 CONTINUE
60 CONTINUE
ELSE
DO 80, J = 1, N
DO 70, I = J, N
C( I, J ) = BETA*C( I, J )
70 CONTINUE
80 CONTINUE
END IF
END IF
RETURN
END IF
*
* Start the operations.
*
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form C := alpha*A*B' + alpha*B*A' + C.
*
IF( UPPER )THEN
DO 130, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 90, I = 1, J
C( I, J ) = ZERO
90 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 100, I = 1, J
C( I, J ) = BETA*C( I, J )
100 CONTINUE
END IF
DO 120, L = 1, K
IF( ( A( J, L ).NE.ZERO ).OR.
$ ( B( J, L ).NE.ZERO ) )THEN
TEMP1 = ALPHA*B( J, L )
TEMP2 = ALPHA*A( J, L )
DO 110, I = 1, J
C( I, J ) = C( I, J ) + A( I, L )*TEMP1 +
$ B( I, L )*TEMP2
110 CONTINUE
END IF
120 CONTINUE
130 CONTINUE
ELSE
DO 180, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 140, I = J, N
C( I, J ) = ZERO
140 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 150, I = J, N
C( I, J ) = BETA*C( I, J )
150 CONTINUE
END IF
DO 170, L = 1, K
IF( ( A( J, L ).NE.ZERO ).OR.
$ ( B( J, L ).NE.ZERO ) )THEN
TEMP1 = ALPHA*B( J, L )
TEMP2 = ALPHA*A( J, L )
DO 160, I = J, N
C( I, J ) = C( I, J ) + A( I, L )*TEMP1 +
$ B( I, L )*TEMP2
160 CONTINUE
END IF
170 CONTINUE
180 CONTINUE
END IF
ELSE
*
* Form C := alpha*A'*B + alpha*B'*A + C.
*
IF( UPPER )THEN
DO 210, J = 1, N
DO 200, I = 1, J
TEMP1 = ZERO
TEMP2 = ZERO
DO 190, L = 1, K
TEMP1 = TEMP1 + A( L, I )*B( L, J )
TEMP2 = TEMP2 + B( L, I )*A( L, J )
190 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2
ELSE
C( I, J ) = BETA *C( I, J ) +
$ ALPHA*TEMP1 + ALPHA*TEMP2
END IF
200 CONTINUE
210 CONTINUE
ELSE
DO 240, J = 1, N
DO 230, I = J, N
TEMP1 = ZERO
TEMP2 = ZERO
DO 220, L = 1, K
TEMP1 = TEMP1 + A( L, I )*B( L, J )
TEMP2 = TEMP2 + B( L, I )*A( L, J )
220 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2
ELSE
C( I, J ) = BETA *C( I, J ) +
$ ALPHA*TEMP1 + ALPHA*TEMP2
END IF
230 CONTINUE
240 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZSYR2K.
*
END
SUBROUTINE ZSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA,
$ BETA, C, LDC )
* .. Scalar Arguments ..
CHARACTER*1 UPLO, TRANS
INTEGER N, K, LDA, LDC
COMPLEX*16 ALPHA, BETA
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), C( LDC, * )
* ..
*
* Purpose
* =======
*
* ZSYRK performs one of the symmetric rank k operations
*
* C := alpha*A*A' + beta*C,
*
* or
*
* C := alpha*A'*A + beta*C,
*
* where alpha and beta are scalars, C is an n by n symmetric matrix
* and A is an n by k matrix in the first case and a k by n matrix
* in the second case.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the array C is to be referenced as
* follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of C
* is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of C
* is to be referenced.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
*
* TRANS = 'T' or 't' C := alpha*A'*A + beta*C.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix C. N must be
* at least zero.
* Unchanged on exit.
*
* K - INTEGER.
* On entry with TRANS = 'N' or 'n', K specifies the number
* of columns of the matrix A, and on entry with
* TRANS = 'T' or 't', K specifies the number of rows of the
* matrix A. K must be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
* k when TRANS = 'N' or 'n', and is n otherwise.
* Before entry with TRANS = 'N' or 'n', the leading n by k
* part of the array A must contain the matrix A, otherwise
* the leading k by n part of the array A must contain the
* matrix A.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When TRANS = 'N' or 'n'
* then LDA must be at least max( 1, n ), otherwise LDA must
* be at least max( 1, k ).
* Unchanged on exit.
*
* BETA - COMPLEX*16 .
* On entry, BETA specifies the scalar beta.
* Unchanged on exit.
*
* C - COMPLEX*16 array of DIMENSION ( LDC, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array C must contain the upper
* triangular part of the symmetric matrix and the strictly
* lower triangular part of C is not referenced. On exit, the
* upper triangular part of the array C is overwritten by the
* upper triangular part of the updated matrix.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array C must contain the lower
* triangular part of the symmetric matrix and the strictly
* upper triangular part of C is not referenced. On exit, the
* lower triangular part of the array C is overwritten by the
* lower triangular part of the updated matrix.
*
* LDC - INTEGER.
* On entry, LDC specifies the first dimension of C as declared
* in the calling (sub) program. LDC must be at least
* max( 1, n ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I, INFO, J, L, NROWA
COMPLEX*16 TEMP
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
IF( LSAME( TRANS, 'N' ) )THEN
NROWA = N
ELSE
NROWA = K
END IF
UPPER = LSAME( UPLO, 'U' )
*
INFO = 0
IF( ( .NOT.UPPER ).AND.
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
$ ( .NOT.LSAME( TRANS, 'T' ) ) )THEN
INFO = 2
ELSE IF( N .LT.0 )THEN
INFO = 3
ELSE IF( K .LT.0 )THEN
INFO = 4
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 7
ELSE IF( LDC.LT.MAX( 1, N ) )THEN
INFO = 10
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZSYRK ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.
$ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
IF( UPPER )THEN
IF( BETA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, J
C( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40, J = 1, N
DO 30, I = 1, J
C( I, J ) = BETA*C( I, J )
30 CONTINUE
40 CONTINUE
END IF
ELSE
IF( BETA.EQ.ZERO )THEN
DO 60, J = 1, N
DO 50, I = J, N
C( I, J ) = ZERO
50 CONTINUE
60 CONTINUE
ELSE
DO 80, J = 1, N
DO 70, I = J, N
C( I, J ) = BETA*C( I, J )
70 CONTINUE
80 CONTINUE
END IF
END IF
RETURN
END IF
*
* Start the operations.
*
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form C := alpha*A*A' + beta*C.
*
IF( UPPER )THEN
DO 130, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 90, I = 1, J
C( I, J ) = ZERO
90 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 100, I = 1, J
C( I, J ) = BETA*C( I, J )
100 CONTINUE
END IF
DO 120, L = 1, K
IF( A( J, L ).NE.ZERO )THEN
TEMP = ALPHA*A( J, L )
DO 110, I = 1, J
C( I, J ) = C( I, J ) + TEMP*A( I, L )
110 CONTINUE
END IF
120 CONTINUE
130 CONTINUE
ELSE
DO 180, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 140, I = J, N
C( I, J ) = ZERO
140 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 150, I = J, N
C( I, J ) = BETA*C( I, J )
150 CONTINUE
END IF
DO 170, L = 1, K
IF( A( J, L ).NE.ZERO )THEN
TEMP = ALPHA*A( J, L )
DO 160, I = J, N
C( I, J ) = C( I, J ) + TEMP*A( I, L )
160 CONTINUE
END IF
170 CONTINUE
180 CONTINUE
END IF
ELSE
*
* Form C := alpha*A'*A + beta*C.
*
IF( UPPER )THEN
DO 210, J = 1, N
DO 200, I = 1, J
TEMP = ZERO
DO 190, L = 1, K
TEMP = TEMP + A( L, I )*A( L, J )
190 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
200 CONTINUE
210 CONTINUE
ELSE
DO 240, J = 1, N
DO 230, I = J, N
TEMP = ZERO
DO 220, L = 1, K
TEMP = TEMP + A( L, I )*A( L, J )
220 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
230 CONTINUE
240 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZSYRK .
*
END
SUBROUTINE ZTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, K, LDA, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * )
* ..
*
* Purpose
* =======
*
* ZTBMV performs one of the matrix-vector operations
*
* x := A*x, or x := A'*x, or x := conjg( A' )*x,
*
* where x is an n element vector and A is an n by n unit, or non-unit,
* upper or lower triangular band matrix, with ( k + 1 ) diagonals.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' x := A*x.
*
* TRANS = 'T' or 't' x := A'*x.
*
* TRANS = 'C' or 'c' x := conjg( A' )*x.
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit
* triangular as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* K - INTEGER.
* On entry with UPLO = 'U' or 'u', K specifies the number of
* super-diagonals of the matrix A.
* On entry with UPLO = 'L' or 'l', K specifies the number of
* sub-diagonals of the matrix A.
* K must satisfy 0 .le. K.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
* by n part of the array A must contain the upper triangular
* band part of the matrix of coefficients, supplied column by
* column, with the leading diagonal of the matrix in row
* ( k + 1 ) of the array, the first super-diagonal starting at
* position 2 in row k, and so on. The top left k by k triangle
* of the array A is not referenced.
* The following program segment will transfer an upper
* triangular band matrix from conventional full matrix storage
* to band storage:
*
* DO 20, J = 1, N
* M = K + 1 - J
* DO 10, I = MAX( 1, J - K ), J
* A( M + I, J ) = matrix( I, J )
* 10 CONTINUE
* 20 CONTINUE
*
* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
* by n part of the array A must contain the lower triangular
* band part of the matrix of coefficients, supplied column by
* column, with the leading diagonal of the matrix in row 1 of
* the array, the first sub-diagonal starting at position 1 in
* row 2, and so on. The bottom right k by k triangle of the
* array A is not referenced.
* The following program segment will transfer a lower
* triangular band matrix from conventional full matrix storage
* to band storage:
*
* DO 20, J = 1, N
* M = 1 - J
* DO 10, I = J, MIN( N, J + K )
* A( M + I, J ) = matrix( I, J )
* 10 CONTINUE
* 20 CONTINUE
*
* Note that when DIAG = 'U' or 'u' the elements of the array A
* corresponding to the diagonal elements of the matrix are not
* referenced, but are assumed to be unity.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* ( k + 1 ).
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x. On exit, X is overwritten with the
* tranformed vector x.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L
LOGICAL NOCONJ, NOUNIT
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO , 'U' ).AND.
$ .NOT.LSAME( UPLO , 'L' ) )THEN
INFO = 1
ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
$ .NOT.LSAME( TRANS, 'T' ).AND.
$ .NOT.LSAME( TRANS, 'C' ) )THEN
INFO = 2
ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
$ .NOT.LSAME( DIAG , 'N' ) )THEN
INFO = 3
ELSE IF( N.LT.0 )THEN
INFO = 4
ELSE IF( K.LT.0 )THEN
INFO = 5
ELSE IF( LDA.LT.( K + 1 ) )THEN
INFO = 7
ELSE IF( INCX.EQ.0 )THEN
INFO = 9
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZTBMV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = LSAME( DIAG , 'N' )
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF( INCX.LE.0 )THEN
KX = 1 - ( N - 1 )*INCX
ELSE IF( INCX.NE.1 )THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form x := A*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
KPLUS1 = K + 1
IF( INCX.EQ.1 )THEN
DO 20, J = 1, N
IF( X( J ).NE.ZERO )THEN
TEMP = X( J )
L = KPLUS1 - J
DO 10, I = MAX( 1, J - K ), J - 1
X( I ) = X( I ) + TEMP*A( L + I, J )
10 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )*A( KPLUS1, J )
END IF
20 CONTINUE
ELSE
JX = KX
DO 40, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = X( JX )
IX = KX
L = KPLUS1 - J
DO 30, I = MAX( 1, J - K ), J - 1
X( IX ) = X( IX ) + TEMP*A( L + I, J )
IX = IX + INCX
30 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )*A( KPLUS1, J )
END IF
JX = JX + INCX
IF( J.GT.K )
$ KX = KX + INCX
40 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 60, J = N, 1, -1
IF( X( J ).NE.ZERO )THEN
TEMP = X( J )
L = 1 - J
DO 50, I = MIN( N, J + K ), J + 1, -1
X( I ) = X( I ) + TEMP*A( L + I, J )
50 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )*A( 1, J )
END IF
60 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 80, J = N, 1, -1
IF( X( JX ).NE.ZERO )THEN
TEMP = X( JX )
IX = KX
L = 1 - J
DO 70, I = MIN( N, J + K ), J + 1, -1
X( IX ) = X( IX ) + TEMP*A( L + I, J )
IX = IX - INCX
70 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )*A( 1, J )
END IF
JX = JX - INCX
IF( ( N - J ).GE.K )
$ KX = KX - INCX
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := A'*x or x := conjg( A' )*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
KPLUS1 = K + 1
IF( INCX.EQ.1 )THEN
DO 110, J = N, 1, -1
TEMP = X( J )
L = KPLUS1 - J
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*A( KPLUS1, J )
DO 90, I = J - 1, MAX( 1, J - K ), -1
TEMP = TEMP + A( L + I, J )*X( I )
90 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( A( KPLUS1, J ) )
DO 100, I = J - 1, MAX( 1, J - K ), -1
TEMP = TEMP + DCONJG( A( L + I, J ) )*X( I )
100 CONTINUE
END IF
X( J ) = TEMP
110 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 140, J = N, 1, -1
TEMP = X( JX )
KX = KX - INCX
IX = KX
L = KPLUS1 - J
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*A( KPLUS1, J )
DO 120, I = J - 1, MAX( 1, J - K ), -1
TEMP = TEMP + A( L + I, J )*X( IX )
IX = IX - INCX
120 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( A( KPLUS1, J ) )
DO 130, I = J - 1, MAX( 1, J - K ), -1
TEMP = TEMP + DCONJG( A( L + I, J ) )*X( IX )
IX = IX - INCX
130 CONTINUE
END IF
X( JX ) = TEMP
JX = JX - INCX
140 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 170, J = 1, N
TEMP = X( J )
L = 1 - J
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*A( 1, J )
DO 150, I = J + 1, MIN( N, J + K )
TEMP = TEMP + A( L + I, J )*X( I )
150 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( A( 1, J ) )
DO 160, I = J + 1, MIN( N, J + K )
TEMP = TEMP + DCONJG( A( L + I, J ) )*X( I )
160 CONTINUE
END IF
X( J ) = TEMP
170 CONTINUE
ELSE
JX = KX
DO 200, J = 1, N
TEMP = X( JX )
KX = KX + INCX
IX = KX
L = 1 - J
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*A( 1, J )
DO 180, I = J + 1, MIN( N, J + K )
TEMP = TEMP + A( L + I, J )*X( IX )
IX = IX + INCX
180 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( A( 1, J ) )
DO 190, I = J + 1, MIN( N, J + K )
TEMP = TEMP + DCONJG( A( L + I, J ) )*X( IX )
IX = IX + INCX
190 CONTINUE
END IF
X( JX ) = TEMP
JX = JX + INCX
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ZTBMV .
*
END
SUBROUTINE ZTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, K, LDA, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), X( * )
* ..
*
* Purpose
* =======
*
* ZTBSV solves one of the systems of equations
*
* A*x = b, or A'*x = b, or conjg( A' )*x = b,
*
* where b and x are n element vectors and A is an n by n unit, or
* non-unit, upper or lower triangular band matrix, with ( k + 1 )
* diagonals.
*
* No test for singularity or near-singularity is included in this
* routine. Such tests must be performed before calling this routine.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the equations to be solved as
* follows:
*
* TRANS = 'N' or 'n' A*x = b.
*
* TRANS = 'T' or 't' A'*x = b.
*
* TRANS = 'C' or 'c' conjg( A' )*x = b.
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit
* triangular as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* K - INTEGER.
* On entry with UPLO = 'U' or 'u', K specifies the number of
* super-diagonals of the matrix A.
* On entry with UPLO = 'L' or 'l', K specifies the number of
* sub-diagonals of the matrix A.
* K must satisfy 0 .le. K.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, n ).
* Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
* by n part of the array A must contain the upper triangular
* band part of the matrix of coefficients, supplied column by
* column, with the leading diagonal of the matrix in row
* ( k + 1 ) of the array, the first super-diagonal starting at
* position 2 in row k, and so on. The top left k by k triangle
* of the array A is not referenced.
* The following program segment will transfer an upper
* triangular band matrix from conventional full matrix storage
* to band storage:
*
* DO 20, J = 1, N
* M = K + 1 - J
* DO 10, I = MAX( 1, J - K ), J
* A( M + I, J ) = matrix( I, J )
* 10 CONTINUE
* 20 CONTINUE
*
* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
* by n part of the array A must contain the lower triangular
* band part of the matrix of coefficients, supplied column by
* column, with the leading diagonal of the matrix in row 1 of
* the array, the first sub-diagonal starting at position 1 in
* row 2, and so on. The bottom right k by k triangle of the
* array A is not referenced.
* The following program segment will transfer a lower
* triangular band matrix from conventional full matrix storage
* to band storage:
*
* DO 20, J = 1, N
* M = 1 - J
* DO 10, I = J, MIN( N, J + K )
* A( M + I, J ) = matrix( I, J )
* 10 CONTINUE
* 20 CONTINUE
*
* Note that when DIAG = 'U' or 'u' the elements of the array A
* corresponding to the diagonal elements of the matrix are not
* referenced, but are assumed to be unity.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* ( k + 1 ).
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element right-hand side vector b. On exit, X is overwritten
* with the solution vector x.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L
LOGICAL NOCONJ, NOUNIT
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO , 'U' ).AND.
$ .NOT.LSAME( UPLO , 'L' ) )THEN
INFO = 1
ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
$ .NOT.LSAME( TRANS, 'T' ).AND.
$ .NOT.LSAME( TRANS, 'C' ) )THEN
INFO = 2
ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
$ .NOT.LSAME( DIAG , 'N' ) )THEN
INFO = 3
ELSE IF( N.LT.0 )THEN
INFO = 4
ELSE IF( K.LT.0 )THEN
INFO = 5
ELSE IF( LDA.LT.( K + 1 ) )THEN
INFO = 7
ELSE IF( INCX.EQ.0 )THEN
INFO = 9
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZTBSV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = LSAME( DIAG , 'N' )
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF( INCX.LE.0 )THEN
KX = 1 - ( N - 1 )*INCX
ELSE IF( INCX.NE.1 )THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed by sequentially with one pass through A.
*
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form x := inv( A )*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
KPLUS1 = K + 1
IF( INCX.EQ.1 )THEN
DO 20, J = N, 1, -1
IF( X( J ).NE.ZERO )THEN
L = KPLUS1 - J
IF( NOUNIT )
$ X( J ) = X( J )/A( KPLUS1, J )
TEMP = X( J )
DO 10, I = J - 1, MAX( 1, J - K ), -1
X( I ) = X( I ) - TEMP*A( L + I, J )
10 CONTINUE
END IF
20 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 40, J = N, 1, -1
KX = KX - INCX
IF( X( JX ).NE.ZERO )THEN
IX = KX
L = KPLUS1 - J
IF( NOUNIT )
$ X( JX ) = X( JX )/A( KPLUS1, J )
TEMP = X( JX )
DO 30, I = J - 1, MAX( 1, J - K ), -1
X( IX ) = X( IX ) - TEMP*A( L + I, J )
IX = IX - INCX
30 CONTINUE
END IF
JX = JX - INCX
40 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 60, J = 1, N
IF( X( J ).NE.ZERO )THEN
L = 1 - J
IF( NOUNIT )
$ X( J ) = X( J )/A( 1, J )
TEMP = X( J )
DO 50, I = J + 1, MIN( N, J + K )
X( I ) = X( I ) - TEMP*A( L + I, J )
50 CONTINUE
END IF
60 CONTINUE
ELSE
JX = KX
DO 80, J = 1, N
KX = KX + INCX
IF( X( JX ).NE.ZERO )THEN
IX = KX
L = 1 - J
IF( NOUNIT )
$ X( JX ) = X( JX )/A( 1, J )
TEMP = X( JX )
DO 70, I = J + 1, MIN( N, J + K )
X( IX ) = X( IX ) - TEMP*A( L + I, J )
IX = IX + INCX
70 CONTINUE
END IF
JX = JX + INCX
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := inv( A' )*x or x := inv( conjg( A') )*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
KPLUS1 = K + 1
IF( INCX.EQ.1 )THEN
DO 110, J = 1, N
TEMP = X( J )
L = KPLUS1 - J
IF( NOCONJ )THEN
DO 90, I = MAX( 1, J - K ), J - 1
TEMP = TEMP - A( L + I, J )*X( I )
90 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( KPLUS1, J )
ELSE
DO 100, I = MAX( 1, J - K ), J - 1
TEMP = TEMP - DCONJG( A( L + I, J ) )*X( I )
100 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( A( KPLUS1, J ) )
END IF
X( J ) = TEMP
110 CONTINUE
ELSE
JX = KX
DO 140, J = 1, N
TEMP = X( JX )
IX = KX
L = KPLUS1 - J
IF( NOCONJ )THEN
DO 120, I = MAX( 1, J - K ), J - 1
TEMP = TEMP - A( L + I, J )*X( IX )
IX = IX + INCX
120 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( KPLUS1, J )
ELSE
DO 130, I = MAX( 1, J - K ), J - 1
TEMP = TEMP - DCONJG( A( L + I, J ) )*X( IX )
IX = IX + INCX
130 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( A( KPLUS1, J ) )
END IF
X( JX ) = TEMP
JX = JX + INCX
IF( J.GT.K )
$ KX = KX + INCX
140 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 170, J = N, 1, -1
TEMP = X( J )
L = 1 - J
IF( NOCONJ )THEN
DO 150, I = MIN( N, J + K ), J + 1, -1
TEMP = TEMP - A( L + I, J )*X( I )
150 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( 1, J )
ELSE
DO 160, I = MIN( N, J + K ), J + 1, -1
TEMP = TEMP - DCONJG( A( L + I, J ) )*X( I )
160 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( A( 1, J ) )
END IF
X( J ) = TEMP
170 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 200, J = N, 1, -1
TEMP = X( JX )
IX = KX
L = 1 - J
IF( NOCONJ )THEN
DO 180, I = MIN( N, J + K ), J + 1, -1
TEMP = TEMP - A( L + I, J )*X( IX )
IX = IX - INCX
180 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( 1, J )
ELSE
DO 190, I = MIN( N, J + K ), J + 1, -1
TEMP = TEMP - DCONJG( A( L + I, J ) )*X( IX )
IX = IX - INCX
190 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( A( 1, J ) )
END IF
X( JX ) = TEMP
JX = JX - INCX
IF( ( N - J ).GE.K )
$ KX = KX - INCX
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ZTBSV .
*
END
SUBROUTINE ZTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
COMPLEX*16 AP( * ), X( * )
* ..
*
* Purpose
* =======
*
* ZTPMV performs one of the matrix-vector operations
*
* x := A*x, or x := A'*x, or x := conjg( A' )*x,
*
* where x is an n element vector and A is an n by n unit, or non-unit,
* upper or lower triangular matrix, supplied in packed form.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' x := A*x.
*
* TRANS = 'T' or 't' x := A'*x.
*
* TRANS = 'C' or 'c' x := conjg( A' )*x.
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit
* triangular as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* AP - COMPLEX*16 array of DIMENSION at least
* ( ( n*( n + 1 ) )/2 ).
* Before entry with UPLO = 'U' or 'u', the array AP must
* contain the upper triangular matrix packed sequentially,
* column by column, so that AP( 1 ) contains a( 1, 1 ),
* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
* respectively, and so on.
* Before entry with UPLO = 'L' or 'l', the array AP must
* contain the lower triangular matrix packed sequentially,
* column by column, so that AP( 1 ) contains a( 1, 1 ),
* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
* respectively, and so on.
* Note that when DIAG = 'U' or 'u', the diagonal elements of
* A are not referenced, but are assumed to be unity.
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x. On exit, X is overwritten with the
* tranformed vector x.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, J, JX, K, KK, KX
LOGICAL NOCONJ, NOUNIT
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO , 'U' ).AND.
$ .NOT.LSAME( UPLO , 'L' ) )THEN
INFO = 1
ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
$ .NOT.LSAME( TRANS, 'T' ).AND.
$ .NOT.LSAME( TRANS, 'C' ) )THEN
INFO = 2
ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
$ .NOT.LSAME( DIAG , 'N' ) )THEN
INFO = 3
ELSE IF( N.LT.0 )THEN
INFO = 4
ELSE IF( INCX.EQ.0 )THEN
INFO = 7
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZTPMV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = LSAME( DIAG , 'N' )
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF( INCX.LE.0 )THEN
KX = 1 - ( N - 1 )*INCX
ELSE IF( INCX.NE.1 )THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of AP are
* accessed sequentially with one pass through AP.
*
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form x:= A*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
KK = 1
IF( INCX.EQ.1 )THEN
DO 20, J = 1, N
IF( X( J ).NE.ZERO )THEN
TEMP = X( J )
K = KK
DO 10, I = 1, J - 1
X( I ) = X( I ) + TEMP*AP( K )
K = K + 1
10 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )*AP( KK + J - 1 )
END IF
KK = KK + J
20 CONTINUE
ELSE
JX = KX
DO 40, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = X( JX )
IX = KX
DO 30, K = KK, KK + J - 2
X( IX ) = X( IX ) + TEMP*AP( K )
IX = IX + INCX
30 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )*AP( KK + J - 1 )
END IF
JX = JX + INCX
KK = KK + J
40 CONTINUE
END IF
ELSE
KK = ( N*( N + 1 ) )/2
IF( INCX.EQ.1 )THEN
DO 60, J = N, 1, -1
IF( X( J ).NE.ZERO )THEN
TEMP = X( J )
K = KK
DO 50, I = N, J + 1, -1
X( I ) = X( I ) + TEMP*AP( K )
K = K - 1
50 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )*AP( KK - N + J )
END IF
KK = KK - ( N - J + 1 )
60 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 80, J = N, 1, -1
IF( X( JX ).NE.ZERO )THEN
TEMP = X( JX )
IX = KX
DO 70, K = KK, KK - ( N - ( J + 1 ) ), -1
X( IX ) = X( IX ) + TEMP*AP( K )
IX = IX - INCX
70 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )*AP( KK - N + J )
END IF
JX = JX - INCX
KK = KK - ( N - J + 1 )
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := A'*x or x := conjg( A' )*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
KK = ( N*( N + 1 ) )/2
IF( INCX.EQ.1 )THEN
DO 110, J = N, 1, -1
TEMP = X( J )
K = KK - 1
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*AP( KK )
DO 90, I = J - 1, 1, -1
TEMP = TEMP + AP( K )*X( I )
K = K - 1
90 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( AP( KK ) )
DO 100, I = J - 1, 1, -1
TEMP = TEMP + DCONJG( AP( K ) )*X( I )
K = K - 1
100 CONTINUE
END IF
X( J ) = TEMP
KK = KK - J
110 CONTINUE
ELSE
JX = KX + ( N - 1 )*INCX
DO 140, J = N, 1, -1
TEMP = X( JX )
IX = JX
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*AP( KK )
DO 120, K = KK - 1, KK - J + 1, -1
IX = IX - INCX
TEMP = TEMP + AP( K )*X( IX )
120 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( AP( KK ) )
DO 130, K = KK - 1, KK - J + 1, -1
IX = IX - INCX
TEMP = TEMP + DCONJG( AP( K ) )*X( IX )
130 CONTINUE
END IF
X( JX ) = TEMP
JX = JX - INCX
KK = KK - J
140 CONTINUE
END IF
ELSE
KK = 1
IF( INCX.EQ.1 )THEN
DO 170, J = 1, N
TEMP = X( J )
K = KK + 1
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*AP( KK )
DO 150, I = J + 1, N
TEMP = TEMP + AP( K )*X( I )
K = K + 1
150 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( AP( KK ) )
DO 160, I = J + 1, N
TEMP = TEMP + DCONJG( AP( K ) )*X( I )
K = K + 1
160 CONTINUE
END IF
X( J ) = TEMP
KK = KK + ( N - J + 1 )
170 CONTINUE
ELSE
JX = KX
DO 200, J = 1, N
TEMP = X( JX )
IX = JX
IF( NOCONJ )THEN
IF( NOUNIT )
$ TEMP = TEMP*AP( KK )
DO 180, K = KK + 1, KK + N - J
IX = IX + INCX
TEMP = TEMP + AP( K )*X( IX )
180 CONTINUE
ELSE
IF( NOUNIT )
$ TEMP = TEMP*DCONJG( AP( KK ) )
DO 190, K = KK + 1, KK + N - J
IX = IX + INCX
TEMP = TEMP + DCONJG( AP( K ) )*X( IX )
190 CONTINUE
END IF
X( JX ) = TEMP
JX = JX + INCX
KK = KK + ( N - J + 1 )
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ZTPMV .
*
END
SUBROUTINE ZTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
COMPLEX*16 AP( * ), X( * )
* ..
*
* Purpose
* =======
*
* ZTPSV solves one of the systems of equations
*
* A*x = b, or A'*x = b, or conjg( A' )*x = b,
*
* where b and x are n element vectors and A is an n by n unit, or
* non-unit, upper or lower triangular matrix, supplied in packed form.
*
* No test for singularity or near-singularity is included in this
* routine. Such tests must be performed before calling this routine.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the equations to be solved as
* follows:
*
* TRANS = 'N' or 'n' A*x = b.
*
* TRANS = 'T' or 't' A'*x = b.
*
* TRANS = 'C' or 'c' conjg( A' )*x = b.
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit
* triangular as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* AP - COMPLEX*16 array of DIMENSION at least
* ( ( n*( n + 1 ) )/2 ).
* Before entry with UPLO = 'U' or 'u', the array AP must
* contain the upper triangular matrix packed sequentially,
* column by column, so that AP( 1 ) contains a( 1, 1 ),
* AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
* respectively, and so on.
* Before entry with UPLO = 'L' or 'l', the array AP must
* contain the lower triangular matrix packed sequentially,
* column by column, so that AP( 1 ) contains a( 1, 1 ),
* AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
* respectively, and so on.
* Note that when DIAG = 'U' or 'u', the diagonal elements of
* A are not referenced, but are assumed to be unity.
* Unchanged on exit.
*
* X - COMPLEX*16 array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element right-hand side vector b. On exit, X is overwritten
* with the solution vector x.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
COMPLEX*16 TEMP
INTEGER I, INFO, IX, J, JX, K, KK, KX
LOGICAL NOCONJ, NOUNIT
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO , 'U' ).AND.
$ .NOT.LSAME( UPLO , 'L' ) )THEN
INFO = 1
ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
$ .NOT.LSAME( TRANS, 'T' ).AND.
$ .NOT.LSAME( TRANS, 'C' ) )THEN
INFO = 2
ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
$ .NOT.LSAME( DIAG , 'N' ) )THEN
INFO = 3
ELSE IF( N.LT.0 )THEN
INFO = 4
ELSE IF( INCX.EQ.0 )THEN
INFO = 7
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZTPSV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = LSAME( DIAG , 'N' )
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF( INCX.LE.0 )THEN
KX = 1 - ( N - 1 )*INCX
ELSE IF( INCX.NE.1 )THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of AP are
* accessed sequentially with one pass through AP.
*
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form x := inv( A )*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
KK = ( N*( N + 1 ) )/2
IF( INCX.EQ.1 )THEN
DO 20, J = N, 1, -1
IF( X( J ).NE.ZERO )THEN
IF( NOUNIT )
$ X( J ) = X( J )/AP( KK )
TEMP = X( J )
K = KK - 1
DO 10, I = J - 1, 1, -1
X( I ) = X( I ) - TEMP*AP( K )
K = K - 1
10 CONTINUE
END IF
KK = KK - J
20 CONTINUE
ELSE
JX = KX + ( N - 1 )*INCX
DO 40, J = N, 1, -1
IF( X( JX ).NE.ZERO )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )/AP( KK )
TEMP = X( JX )
IX = JX
DO 30, K = KK - 1, KK - J + 1, -1
IX = IX - INCX
X( IX ) = X( IX ) - TEMP*AP( K )
30 CONTINUE
END IF
JX = JX - INCX
KK = KK - J
40 CONTINUE
END IF
ELSE
KK = 1
IF( INCX.EQ.1 )THEN
DO 60, J = 1, N
IF( X( J ).NE.ZERO )THEN
IF( NOUNIT )
$ X( J ) = X( J )/AP( KK )
TEMP = X( J )
K = KK + 1
DO 50, I = J + 1, N
X( I ) = X( I ) - TEMP*AP( K )
K = K + 1
50 CONTINUE
END IF
KK = KK + ( N - J + 1 )
60 CONTINUE
ELSE
JX = KX
DO 80, J = 1, N
IF( X( JX ).NE.ZERO )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )/AP( KK )
TEMP = X( JX )
IX = JX
DO 70, K = KK + 1, KK + N - J
IX = IX + INCX
X( IX ) = X( IX ) - TEMP*AP( K )
70 CONTINUE
END IF
JX = JX + INCX
KK = KK + ( N - J + 1 )
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := inv( A' )*x or x := inv( conjg( A' ) )*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
KK = 1
IF( INCX.EQ.1 )THEN
DO 110, J = 1, N
TEMP = X( J )
K = KK
IF( NOCONJ )THEN
DO 90, I = 1, J - 1
TEMP = TEMP - AP( K )*X( I )
K = K + 1
90 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/AP( KK + J - 1 )
ELSE
DO 100, I = 1, J - 1
TEMP = TEMP - DCONJG( AP( K ) )*X( I )
K = K + 1
100 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( AP( KK + J - 1 ) )
END IF
X( J ) = TEMP
KK = KK + J
110 CONTINUE
ELSE
JX = KX
DO 140, J = 1, N
TEMP = X( JX )
IX = KX
IF( NOCONJ )THEN
DO 120, K = KK, KK + J - 2
TEMP = TEMP - AP( K )*X( IX )
IX = IX + INCX
120 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/AP( KK + J - 1 )
ELSE
DO 130, K = KK, KK + J - 2
TEMP = TEMP - DCONJG( AP( K ) )*X( IX )
IX = IX + INCX
130 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( AP( KK + J - 1 ) )
END IF
X( JX ) = TEMP
JX = JX + INCX
KK = KK + J
140 CONTINUE
END IF
ELSE
KK = ( N*( N + 1 ) )/2
IF( INCX.EQ.1 )THEN
DO 170, J = N, 1, -1
TEMP = X( J )
K = KK
IF( NOCONJ )THEN
DO 150, I = N, J + 1, -1
TEMP = TEMP - AP( K )*X( I )
K = K - 1
150 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/AP( KK - N + J )
ELSE
DO 160, I = N, J + 1, -1
TEMP = TEMP - DCONJG( AP( K ) )*X( I )
K = K - 1
160 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( AP( KK - N + J ) )
END IF
X( J ) = TEMP
KK = KK - ( N - J + 1 )
170 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 200, J = N, 1, -1
TEMP = X( JX )
IX = KX
IF( NOCONJ )THEN
DO 180, K = KK, KK - ( N - ( J + 1 ) ), -1
TEMP = TEMP - AP( K )*X( IX )
IX = IX - INCX
180 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/AP( KK - N + J )
ELSE
DO 190, K = KK, KK - ( N - ( J + 1 ) ), -1
TEMP = TEMP - DCONJG( AP( K ) )*X( IX )
IX = IX - INCX
190 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/DCONJG( AP( KK - N + J ) )
END IF
X( JX ) = TEMP
JX = JX - INCX
KK = KK - ( N - J + 1 )
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ZTPSV .
*
END
SUBROUTINE ZGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
$ BETA, C, LDC )
* .. Scalar Arguments ..
CHARACTER*1 TRANSA, TRANSB
INTEGER M, N, K, LDA, LDB, LDC
COMPLEX*16 ALPHA, BETA
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
* ..
*
* Purpose
* =======
*
* ZGEMM performs one of the matrix-matrix operations
*
* C := alpha*op( A )*op( B ) + beta*C,
*
* where op( X ) is one of
*
* op( X ) = X or op( X ) = X' or op( X ) = conjg( X' ),
*
* alpha and beta are scalars, and A, B and C are matrices, with op( A )
* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
*
* Parameters
* ==========
*
* TRANSA - CHARACTER*1.
* On entry, TRANSA specifies the form of op( A ) to be used in
* the matrix multiplication as follows:
*
* TRANSA = 'N' or 'n', op( A ) = A.
*
* TRANSA = 'T' or 't', op( A ) = A'.
*
* TRANSA = 'C' or 'c', op( A ) = conjg( A' ).
*
* Unchanged on exit.
*
* TRANSB - CHARACTER*1.
* On entry, TRANSB specifies the form of op( B ) to be used in
* the matrix multiplication as follows:
*
* TRANSB = 'N' or 'n', op( B ) = B.
*
* TRANSB = 'T' or 't', op( B ) = B'.
*
* TRANSB = 'C' or 'c', op( B ) = conjg( B' ).
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of the matrix
* op( A ) and of the matrix C. M must be at least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of the matrix
* op( B ) and the number of columns of the matrix C. N must be
* at least zero.
* Unchanged on exit.
*
* K - INTEGER.
* On entry, K specifies the number of columns of the matrix
* op( A ) and the number of rows of the matrix op( B ). K must
* be at least zero.
* Unchanged on exit.
*
* ALPHA - COMPLEX*16 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
* k when TRANSA = 'N' or 'n', and is m otherwise.
* Before entry with TRANSA = 'N' or 'n', the leading m by k
* part of the array A must contain the matrix A, otherwise
* the leading k by m part of the array A must contain the
* matrix A.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When TRANSA = 'N' or 'n' then
* LDA must be at least max( 1, m ), otherwise LDA must be at
* least max( 1, k ).
* Unchanged on exit.
*
* B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is
* n when TRANSB = 'N' or 'n', and is k otherwise.
* Before entry with TRANSB = 'N' or 'n', the leading k by n
* part of the array B must contain the matrix B, otherwise
* the leading n by k part of the array B must contain the
* matrix B.
* Unchanged on exit.
*
* LDB - INTEGER.
* On entry, LDB specifies the first dimension of B as declared
* in the calling (sub) program. When TRANSB = 'N' or 'n' then
* LDB must be at least max( 1, k ), otherwise LDB must be at
* least max( 1, n ).
* Unchanged on exit.
*
* BETA - COMPLEX*16 .
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then C need not be set on input.
* Unchanged on exit.
*
* C - COMPLEX*16 array of DIMENSION ( LDC, n ).
* Before entry, the leading m by n part of the array C must
* contain the matrix C, except when beta is zero, in which
* case C need not be set on entry.
* On exit, the array C is overwritten by the m by n matrix
* ( alpha*op( A )*op( B ) + beta*C ).
*
* LDC - INTEGER.
* On entry, LDC specifies the first dimension of C as declared
* in the calling (sub) program. LDC must be at least
* max( 1, m ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, MAX
* .. Local Scalars ..
LOGICAL CONJA, CONJB, NOTA, NOTB
INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB
COMPLEX*16 TEMP
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* ..
* .. Executable Statements ..
*
* Set NOTA and NOTB as true if A and B respectively are not
* conjugated or transposed, set CONJA and CONJB as true if A and
* B respectively are to be transposed but not conjugated and set
* NROWA, NCOLA and NROWB as the number of rows and columns of A
* and the number of rows of B respectively.
*
NOTA = LSAME( TRANSA, 'N' )
NOTB = LSAME( TRANSB, 'N' )
CONJA = LSAME( TRANSA, 'C' )
CONJB = LSAME( TRANSB, 'C' )
IF( NOTA )THEN
NROWA = M
NCOLA = K
ELSE
NROWA = K
NCOLA = M
END IF
IF( NOTB )THEN
NROWB = K
ELSE
NROWB = N
END IF
*
* Test the input parameters.
*
INFO = 0
IF( ( .NOT.NOTA ).AND.
$ ( .NOT.CONJA ).AND.
$ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.NOTB ).AND.
$ ( .NOT.CONJB ).AND.
$ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN
INFO = 2
ELSE IF( M .LT.0 )THEN
INFO = 3
ELSE IF( N .LT.0 )THEN
INFO = 4
ELSE IF( K .LT.0 )THEN
INFO = 5
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 8
ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
INFO = 10
ELSE IF( LDC.LT.MAX( 1, M ) )THEN
INFO = 13
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'ZGEMM ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
$ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
IF( BETA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, M
C( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40, J = 1, N
DO 30, I = 1, M
C( I, J ) = BETA*C( I, J )
30 CONTINUE
40 CONTINUE
END IF
RETURN
END IF
*
* Start the operations.
*
IF( NOTB )THEN
IF( NOTA )THEN
*
* Form C := alpha*A*B + beta*C.
*
DO 90, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 50, I = 1, M
C( I, J ) = ZERO
50 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 60, I = 1, M
C( I, J ) = BETA*C( I, J )
60 CONTINUE
END IF
DO 80, L = 1, K
IF( B( L, J ).NE.ZERO )THEN
TEMP = ALPHA*B( L, J )
DO 70, I = 1, M
C( I, J ) = C( I, J ) + TEMP*A( I, L )
70 CONTINUE
END IF
80 CONTINUE
90 CONTINUE
ELSE IF( CONJA )THEN
*
* Form C := alpha*conjg( A' )*B + beta*C.
*
DO 120, J = 1, N
DO 110, I = 1, M
TEMP = ZERO
DO 100, L = 1, K
TEMP = TEMP + DCONJG( A( L, I ) )*B( L, J )
100 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
110 CONTINUE
120 CONTINUE
ELSE
*
* Form C := alpha*A'*B + beta*C
*
DO 150, J = 1, N
DO 140, I = 1, M
TEMP = ZERO
DO 130, L = 1, K
TEMP = TEMP + A( L, I )*B( L, J )
130 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
140 CONTINUE
150 CONTINUE
END IF
ELSE IF( NOTA )THEN
IF( CONJB )THEN
*
* Form C := alpha*A*conjg( B' ) + beta*C.
*
DO 200, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 160, I = 1, M
C( I, J ) = ZERO
160 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 170, I = 1, M
C( I, J ) = BETA*C( I, J )
170 CONTINUE
END IF
DO 190, L = 1, K
IF( B( J, L ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( B( J, L ) )
DO 180, I = 1, M
C( I, J ) = C( I, J ) + TEMP*A( I, L )
180 CONTINUE
END IF
190 CONTINUE
200 CONTINUE
ELSE
*
* Form C := alpha*A*B' + beta*C
*
DO 250, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 210, I = 1, M
C( I, J ) = ZERO
210 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 220, I = 1, M
C( I, J ) = BETA*C( I, J )
220 CONTINUE
END IF
DO 240, L = 1, K
IF( B( J, L ).NE.ZERO )THEN
TEMP = ALPHA*B( J, L )
DO 230, I = 1, M
C( I, J ) = C( I, J ) + TEMP*A( I, L )
230 CONTINUE
END IF
240 CONTINUE
250 CONTINUE
END IF
ELSE IF( CONJA )THEN
IF( CONJB )THEN
*
* Form C := alpha*conjg( A' )*conjg( B' ) + beta*C.
*
DO 280, J = 1, N
DO 270, I = 1, M
TEMP = ZERO
DO 260, L = 1, K
TEMP = TEMP +
$ DCONJG( A( L, I ) )*DCONJG( B( J, L ) )
260 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
270 CONTINUE
280 CONTINUE
ELSE
*
* Form C := alpha*conjg( A' )*B' + beta*C
*
DO 310, J = 1, N
DO 300, I = 1, M
TEMP = ZERO
DO 290, L = 1, K
TEMP = TEMP + DCONJG( A( L, I ) )*B( J, L )
290 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
300 CONTINUE
310 CONTINUE
END IF
ELSE
IF( CONJB )THEN
*
* Form C := alpha*A'*conjg( B' ) + beta*C
*
DO 340, J = 1, N
DO 330, I = 1, M
TEMP = ZERO
DO 320, L = 1, K
TEMP = TEMP + A( L, I )*DCONJG( B( J, L ) )
320 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
330 CONTINUE
340 CONTINUE
ELSE
*
* Form C := alpha*A'*B' + beta*C
*
DO 370, J = 1, N
DO 360, I = 1, M
TEMP = ZERO
DO 350, L = 1, K
TEMP = TEMP + A( L, I )*B( J, L )
350 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
360 CONTINUE
370 CONTINUE
END IF
END IF
*
RETURN
*
* End of ZGEMM .
*
END