*
************************************************************************
*
* File of the EXTENDED-COMPLEX Level-2 BLAS.
* ===================================================
*
* SUBROUTINE ECGEMV( TRANS, M, N, ALPHA, A, LDA, X, INCX,
* $ BETA, Y, INCY )
*
* SUBROUTINE ECGBMV( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
* $ BETA, Y, INCY )
*
* SUBROUTINE ECHEMV( UPLO, N, ALPHA, A, LDA, X, INCX,
* $ BETA, Y, INCY )
*
* SUBROUTINE ECHBMV( UPLO, N, K, ALPHA, A, LDA, X, INCX,
* $ BETA, Y, INCY )
*
* SUBROUTINE ECHPMV( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
*
* SUBROUTINE ECTRMV( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
*
* SUBROUTINE ECTBMV( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
*
* SUBROUTINE ECTPMV( UPLO, TRANS, DIAG, N, AP, X, INCX )
*
* SUBROUTINE ECTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
*
* SUBROUTINE ECTBSV( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
*
* SUBROUTINE ECTPSV( UPLO, TRANS, DIAG, N, AP, X, INCX )
*
* SUBROUTINE ECGERU( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
*
* SUBROUTINE ECGERC( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
*
* SUBROUTINE ECHER ( UPLO, N, ALPHA, X, INCX, A, LDA )
*
* SUBROUTINE ECHPR ( UPLO, N, ALPHA, X, INCX, AP )
*
* SUBROUTINE ECHER2( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
*
* SUBROUTINE ECHPR2( UPLO, N, ALPHA, X, INCX, Y, INCY, AP )
*
* See:
*
* Dongarra J. J., Du Croz J. J., Hammarling S. and Hanson R. J..
* A proposal for an extended set of Fortran Basic Linear Algebra
* Subprograms.
*
* Technical Memorandum No.41 (revision 1), Mathematics and
* Computer Science Division, Argonne National Laboratory, 9700
* South Cass Avenue, Argonne, Illinois 60439, US.
*
* Or
*
* NAG Technical Report TR3/86, Numerical Algorithms Group Ltd.,
* NAG Central Office, 256 Banbury Road, Oxford OX2 7DE, UK, and
* Numerical Algorithms Group Inc., 1101 31st Street, Suite 100,
* Downers Grove, Illinois 60515-1263, USA.
*
************************************************************************
*
SUBROUTINE ECGEMV( TRANS, M, N, ALPHA, A, LDA, X, INCX,
$ BETA, Y, INCY )
* .. Scalar Arguments ..
COMPLEX ALPHA, BETA
INTEGER INCX, INCY, LDA, M, N
CHARACTER*1 TRANS
* .. Array Arguments ..
COMPLEX*16 Y( * )
COMPLEX A( LDA, * ), X( * )
* ..
*
* Purpose
* =======
*
* ECGEMV 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. Additional precision arithmetic is used in the
* computation.
*
* 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 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX 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 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 .
* 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. At least double precision arithmetic is
* used in the computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX ONE
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
COMPLEX ZERO
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+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 DCMPLX, CONJG, 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( 'ECGEMV', 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.
*
IF( LSAME( TRANS, 'N' ) )THEN
LENX = N
LENY = M
ELSE
LENX = M
LENY = N
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 and set up the start points in X and Y if
* the increments are not both unity.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
IF( BETA.NE.ONE )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
END IF
ELSE
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
IF( BETA.NE.ONE )THEN
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.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 60, J = 1, N
IF( X( J ).NE.ZERO )THEN
TEMP = DCMPLX( ALPHA )*X( J )
DO 50, I = 1, M
Y( I ) = Y( 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
TEMP = DCMPLX( 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.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 110, J = 1, N
TEMP = ZERO
IF( NOCONJ )THEN
DO 90, I = 1, M
TEMP = TEMP + A( I, J )*DCMPLX( X( I ) )
90 CONTINUE
ELSE
DO 100, I = 1, M
TEMP = TEMP + CONJG( A( I, J ) )*DCMPLX( X( I ) )
100 CONTINUE
END IF
Y( J ) = Y( J ) + ALPHA*TEMP
110 CONTINUE
ELSE
JY = KY
DO 140, J = 1, N
TEMP = ZERO
IX = KX
IF( NOCONJ )THEN
DO 120, I = 1, M
TEMP = TEMP + A( I, J )*DCMPLX( X( IX ) )
IX = IX + INCX
120 CONTINUE
ELSE
DO 130, I = 1, M
TEMP = TEMP + CONJG( A( I, J ) )*DCMPLX( 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 ECGEMV.
*
END
*
************************************************************************
*
SUBROUTINE ECGERC( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
* .. Scalar Arguments ..
COMPLEX ALPHA
INTEGER INCX, INCY, LDA, M, N
* .. Array Arguments ..
COMPLEX*16 X( * ), Y( * )
COMPLEX A( LDA, * )
* ..
*
* Purpose
* =======
*
* ECGERC 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. Additional precision is used in the
* computation.
*
* 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 .
* 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 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. At least double precision
* arithmetic is used in the computation of A.
*
* 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 20-July-1986.
* 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 Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC DCONJG, CMPLX, 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( 'ECGERC', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.CMPLX( ZERO ) ) )
$ RETURN
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 20, J = 1, N
IF( Y( J ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( Y( J ) )
DO 10, I = 1, M
A( I, J ) = A( I, J ) + X( I )*TEMP
10 CONTINUE
END IF
20 CONTINUE
ELSE
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( M - 1 )*INCX
END IF
IF( INCY.GT.0 )THEN
JY = 1
ELSE
JY = 1 - ( N - 1 )*INCY
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 ECGERC.
*
END
*
************************************************************************
*
SUBROUTINE ECGERU( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
* .. Scalar Arguments ..
COMPLEX ALPHA
INTEGER INCX, INCY, LDA, M, N
* .. Array Arguments ..
COMPLEX*16 X( * ), Y( * )
COMPLEX A( LDA, * )
* ..
*
* Purpose
* =======
*
* ECGERU 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. Additional precision is used in the
* computation.
*
* 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 .
* 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 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. At least double precision
* arithmetic is used in the computation of A.
*
* 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 20-July-1986.
* 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 Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC CMPLX, 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( 'ECGERU', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.CMPLX( ZERO ) ) )
$ RETURN
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 20, J = 1, N
IF( Y( J ).NE.ZERO )THEN
TEMP = ALPHA*Y( J )
DO 10, I = 1, M
A( I, J ) = A( I, J ) + X( I )*TEMP
10 CONTINUE
END IF
20 CONTINUE
ELSE
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( M - 1 )*INCX
END IF
IF( INCY.GT.0 )THEN
JY = 1
ELSE
JY = 1 - ( N - 1 )*INCY
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 ECGERU.
*
END
*
************************************************************************
*
SUBROUTINE ECGBMV( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
$ BETA, Y, INCY )
* .. Scalar Arguments ..
COMPLEX ALPHA, BETA
INTEGER INCX, INCY, KL, KU, LDA, M, N
CHARACTER*1 TRANS
* .. Array Arguments ..
COMPLEX*16 Y( * )
COMPLEX A( LDA, * ), X( * )
* ..
*
* Purpose
* =======
*
* ECGBMV 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.
* Additional precision arithmetic is used in the computation.
*
* 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 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX 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 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 .
* 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.
* At least double precision arithmetic is used in the
* computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX ONE
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
COMPLEX ZERO
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+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 DCMPLX, MAX, MIN, CONJG
* ..
* .. 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( 'ECGBMV', 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.
*
IF( LSAME( TRANS, 'N' ) )THEN
LENX = N
LENY = M
ELSE
LENX = M
LENY = N
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 and set up the start points in X and Y
* if the increments are not both unity.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
IF( BETA.NE.ONE )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
END IF
ELSE
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
IF( BETA.NE.ONE )THEN
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.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
DO 60, J = 1, N
IF( X( J ).NE.ZERO )THEN
TEMP = DCMPLX( ALPHA )*X( J )
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
60 CONTINUE
ELSE
JX = KX
DO 80, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = DCMPLX( 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.
*
IF( ( INCX.EQ.1 ).AND.( INCY.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 )*DCMPLX( X( I ) )
90 CONTINUE
ELSE
DO 100, I = MAX( 1, J - KU ), MIN( M, J + KL )
TEMP = TEMP +
$ CONJG( A( K + I, J ) )*DCMPLX( X( I ) )
100 CONTINUE
END IF
Y( J ) = Y( J ) + ALPHA*TEMP
110 CONTINUE
ELSE
JY = KY
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 )*DCMPLX( X( IX ) )
IX = IX + INCX
120 CONTINUE
ELSE
DO 130, I = MAX( 1, J - KU ), MIN( M, J + KL )
TEMP = TEMP +
$ CONJG( A( K + I, J ) )*DCMPLX( 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 ECGBMV.
*
END
*
************************************************************************
*
SUBROUTINE ECHEMV( UPLO, N, ALPHA, A, LDA, X, INCX,
$ BETA, Y, INCY )
* .. Scalar Arguments ..
COMPLEX ALPHA, BETA
INTEGER INCX, INCY, LDA, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 Y( * )
COMPLEX A( LDA, * ), X( * )
* ..
*
* Purpose
* =======
*
* ECHEMV 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. Additional precision arithmetic is
* used in the computation.
*
* 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 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX 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 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 .
* 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. At least double precision arithmetic is used in
* the computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX ONE
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
COMPLEX ZERO
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+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 DCMPLX, CONJG, MAX, REAL
* ..
* .. 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( 'ECHEMV', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* 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 and set up the start points in X and Y if
* the increments are not both unity.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
IF( BETA.NE.ONE )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
END IF
ELSE
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
IF( BETA.NE.ONE )THEN
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 = DCMPLX( ALPHA )*X( J )
TEMP2 = ZERO
DO 50, I = 1, J - 1
Y( I ) = Y( I ) + TEMP1*A( I, J )
TEMP2 = TEMP2 + CONJG( A( I, J ) )*DCMPLX( X( I ) )
50 CONTINUE
Y( J ) = Y( J ) + TEMP1*REAL( A( J, J ) ) + ALPHA*TEMP2
60 CONTINUE
ELSE
IX = KX - INCX
DO 80, J = 1, N
TEMP1 = DCMPLX( ALPHA )*X( IX + INCX )
TEMP2 = ZERO
IX = KX
IY = KY
DO 70, I = 1, J - 1
Y( IY ) = Y( IY ) + TEMP1*A( I, J )
TEMP2 = TEMP2 +
$ CONJG( A( I, J ) )*DCMPLX( X( IX ) )
IX = IX + INCX
IY = IY + INCY
70 CONTINUE
Y( IY ) = Y( IY ) + TEMP1*REAL( A( J, J ) ) + ALPHA*TEMP2
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 = DCMPLX( ALPHA )*X( J )
TEMP2 = ZERO
Y( J ) = Y( J ) + TEMP1*REAL( A( J, J ) )
DO 90, I = J + 1, N
Y( I ) = Y( I ) + TEMP1*A( I, J )
TEMP2 = TEMP2 + CONJG( A( I, J ) )*DCMPLX( X( I ) )
90 CONTINUE
Y( J ) = Y( J ) + ALPHA*TEMP2
100 CONTINUE
ELSE
JX = KX
JY = KY
DO 120, J = 1, N
TEMP1 = DCMPLX( ALPHA )*X( JX )
TEMP2 = ZERO
Y( JY ) = Y( JY ) + TEMP1*REAL( 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 +
$ CONJG( A( I, J ) )*DCMPLX( 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 ECHEMV.
*
END
*
************************************************************************
*
SUBROUTINE ECHER ( UPLO, N, ALPHA, X, INCX, A, LDA )
* .. Scalar Arguments ..
REAL ALPHA
INTEGER INCX, LDA, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 X( * )
COMPLEX A( LDA, * )
* ..
*
* Purpose
* =======
*
* ECHER 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. Additional precision arithmetic is used in
* the computation.
*
* 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 - REAL .
* 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 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. At least double precision arithmetic is
* used in the computation of A.
*
* 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 20-July-1986.
* 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, DBLE, MAX, REAL
* ..
* .. 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( 'ECHER ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ALPHA.EQ.REAL( 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 ) = REAL( A( J, J ) ) + DBLE( X( J )*TEMP )
ELSE
A( J, J ) = REAL( 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 ) = REAL( A( J, J ) ) + DBLE( X( JX )*TEMP )
ELSE
A( J, J ) = REAL( 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 ) = REAL( 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 ) = REAL( 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 ) = REAL( 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 ) = REAL( A( J, J ) )
END IF
JX = JX + INCX
80 CONTINUE
END IF
END IF
*
RETURN
*
* End of ECHER .
*
END
*
************************************************************************
*
SUBROUTINE ECHER2( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
* .. Scalar Arguments ..
COMPLEX ALPHA
INTEGER INCX, INCY, LDA, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 X( * ), Y( * )
COMPLEX A( LDA, * )
* ..
*
* Purpose
* =======
*
* ECHER2 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. Additional precision arithmetic is used in the
* computation.
*
* 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 .
* 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 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. At least double precision arithmetic is
* used in the computation of A.
*
* 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 20-July-1986.
* 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, DBLE, CMPLX, MAX, REAL
* ..
* .. 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( 'ECHER2', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ALPHA.EQ.CMPLX( 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
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 ) = REAL( A( J, J ) ) +
$ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
ELSE
A( J, J ) = REAL( A( J, J ) )
END IF
20 CONTINUE
ELSE
JX = KX
JY = KY
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 ) = REAL( A( J, J ) ) +
$ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
ELSE
A( J, J ) = REAL( 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 ) = REAL( 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 ) = REAL( A( J, J ) )
END IF
60 CONTINUE
ELSE
JX = KX
JY = KY
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 ) = REAL( 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 ) = REAL( A( J, J ) )
END IF
JX = JX + INCX
JY = JY + INCY
80 CONTINUE
END IF
END IF
*
RETURN
*
* End of ECHER2.
*
END
*
************************************************************************
*
SUBROUTINE ECHPMV( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
* .. Scalar Arguments ..
COMPLEX ALPHA, BETA
INTEGER INCX, INCY, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 Y( * )
COMPLEX AP( * ), X( * )
* ..
*
* Purpose
* =======
*
* ECHPMV 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. Additional precision arithmetic is
* used in the computation.
*
* 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 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* AP - COMPLEX 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 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 .
* 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. At least double precision arithmetic is used in
* the computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX ONE
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
COMPLEX ZERO
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+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 DCMPLX, CONJG, REAL
* ..
* .. 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( 'ECHPMV', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* 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 and set up the start points in X and Y if
* the increments are not both unity.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
IF( BETA.NE.ONE )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
END IF
ELSE
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
IF( BETA.NE.ONE )THEN
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
K = 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 = DCMPLX( ALPHA )*X( J )
TEMP2 = ZERO
DO 50, I = 1, J - 1
Y( I ) = Y( I ) + TEMP1*AP( K )
TEMP2 = TEMP2 + CONJG( AP( K ) )*DCMPLX( X( I ) )
K = K + 1
50 CONTINUE
Y( J ) = Y( J ) + TEMP1*REAL( AP( K ) ) + ALPHA*TEMP2
K = K + 1
60 CONTINUE
ELSE
IX = KX - INCX
DO 80, J = 1, N
TEMP1 = DCMPLX( ALPHA )*X( IX + INCX )
TEMP2 = ZERO
IX = KX
IY = KY
KK = K
DO 70, K = KK, KK + J - 2
Y( IY ) = Y( IY ) + TEMP1*AP( K )
TEMP2 = TEMP2 + CONJG( AP( K ) )*DCMPLX( X( IX ) )
IX = IX + INCX
IY = IY + INCY
70 CONTINUE
Y( IY ) = Y( IY ) + TEMP1*REAL( AP( K ) ) + ALPHA*TEMP2
K = K + 1
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 = DCMPLX( ALPHA )*X( J )
TEMP2 = ZERO
Y( J ) = Y( J ) + TEMP1*REAL( AP( K ) )
K = K + 1
DO 90, I = J + 1, N
Y( I ) = Y( I ) + TEMP1*AP( K )
TEMP2 = TEMP2 + CONJG( AP( K ) )*DCMPLX( X( I ) )
K = K + 1
90 CONTINUE
Y( J ) = Y( J ) + ALPHA*TEMP2
100 CONTINUE
ELSE
JX = KX
JY = KY
DO 120, J = 1, N
TEMP1 = DCMPLX( ALPHA )*X( JX )
TEMP2 = ZERO
Y( JY ) = Y( JY ) + TEMP1*REAL( AP( K ) )
IX = JX
IY = JY
KK = K + 1
DO 110, K = KK, KK + N - ( J + 1 )
IX = IX + INCX
IY = IY + INCY
Y( IY ) = Y( IY ) + TEMP1*AP( K )
TEMP2 = TEMP2 + CONJG( AP( K ) )*DCMPLX( 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 ECHPMV.
*
END
*
************************************************************************
*
SUBROUTINE ECHPR ( UPLO, N, ALPHA, X, INCX, AP )
* .. Scalar Arguments ..
REAL ALPHA
INTEGER INCX, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 X( * )
COMPLEX AP( * )
* ..
*
* Purpose
* =======
*
* ECHPR 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. Additional precision arithmetic is used in
* the computation.
*
* 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 - REAL .
* 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 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. At least double precision arithmetic is
* used in the computation of A.
*
*
* Level 2 Blas routine.
*
* -- Written on 20-July-1986.
* 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, REAL
* ..
* .. 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( 'ECHPR ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ALPHA.EQ.REAL( 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.
*
K = 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 ) )
DO 10, I = 1, J - 1
AP( K ) = AP( K ) + X( I )*TEMP
K = K + 1
10 CONTINUE
AP( K ) = REAL( AP( K ) ) + DBLE( X( J )*TEMP )
ELSE
K = K + J - 1
AP( K ) = REAL( AP( K ) )
END IF
K = K + 1
20 CONTINUE
ELSE
JX = KX
DO 40, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( X( JX ) )
IX = KX
KK = K
DO 30, K = KK, KK + J - 2
AP( K ) = AP( K ) + X( IX )*TEMP
IX = IX + INCX
30 CONTINUE
AP( K ) = REAL( AP( K ) ) + DBLE( X( JX )*TEMP )
ELSE
K = K + J - 1
AP( K ) = REAL( AP( K ) )
END IF
K = K + 1
JX = JX + INCX
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( K ) = REAL( AP( K ) ) + DBLE( TEMP*X( J ) )
K = K + 1
DO 50, I = J + 1, N
AP( K ) = AP( K ) + X( I )*TEMP
K = K + 1
50 CONTINUE
ELSE
AP( K ) = REAL( AP( K ) )
K = K + N - J + 1
END IF
60 CONTINUE
ELSE
JX = KX
DO 80, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*DCONJG( X( JX ) )
AP( K ) = REAL( AP( K ) ) + DBLE( TEMP*X( JX ) )
IX = JX
KK = K + 1
DO 70, K = KK, KK + N - J - 1
IX = IX + INCX
AP( K ) = AP( K ) + X( IX )*TEMP
70 CONTINUE
ELSE
AP( K ) = REAL( AP( K ) )
K = K + N - J + 1
END IF
JX = JX + INCX
80 CONTINUE
END IF
END IF
*
RETURN
*
* End of ECHPR .
*
END
*
************************************************************************
*
SUBROUTINE ECHPR2( UPLO, N, ALPHA, X, INCX, Y, INCY, AP )
* .. Scalar Arguments ..
COMPLEX ALPHA
INTEGER INCX, INCY, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 X( * ), Y( * )
COMPLEX AP( * )
* ..
*
* Purpose
* =======
*
* ECHPR2 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. Additional precision arithmetic is used in
* the computation.
*
* 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 .
* 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 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. At least double precision arithmetic is
* used in the computation of A.
*
*
* Level 2 Blas routine.
*
* -- Written on 20-July-1986.
* 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, CMPLX, REAL
* ..
* .. 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( 'ECHPR2', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ALPHA.EQ.CMPLX( 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
END IF
*
* Start the operations. In this version the elements of the array AP
* are accessed sequentially with one pass through AP.
*
K = 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 ) )
DO 10, I = 1, J - 1
AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2
K = K + 1
10 CONTINUE
AP( K ) = REAL( AP( K ) ) +
$ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
ELSE
K = K + J - 1
AP( K ) = REAL( AP( K ) )
END IF
K = K + 1
20 CONTINUE
ELSE
JX = KX
JY = KY
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
KK = K
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( K ) = REAL( AP( K ) ) +
$ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
ELSE
K = K + J - 1
AP( K ) = REAL( AP( K ) )
END IF
K = K + 1
JX = JX + INCX
JY = JY + INCY
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( K ) = REAL( AP( K ) ) +
$ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
K = K + 1
DO 50, I = J + 1, N
AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2
K = K + 1
50 CONTINUE
ELSE
AP( K ) = REAL( AP( K ) )
K = K + N - J + 1
END IF
60 CONTINUE
ELSE
JX = KX
JY = KY
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( K ) = REAL( AP( K ) ) +
$ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
IX = JX
IY = JY
KK = K + 1
DO 70, K = KK, KK + N - J - 1
IX = IX + INCX
IY = IY + INCY
AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2
70 CONTINUE
ELSE
AP( K ) = REAL( AP( K ) )
K = K + N - J + 1
END IF
JX = JX + INCX
JY = JY + INCY
80 CONTINUE
END IF
END IF
*
RETURN
*
* End of ECHPR2.
*
END
*
************************************************************************
*
SUBROUTINE ECHBMV( UPLO, N, K, ALPHA, A, LDA, X, INCX,
$ BETA, Y, INCY )
* .. Scalar Arguments ..
COMPLEX ALPHA, BETA
INTEGER INCX, INCY, K, LDA, N
CHARACTER*1 UPLO
* .. Array Arguments ..
COMPLEX*16 Y( * )
COMPLEX A( LDA, * ), X( * )
* ..
*
* Purpose
* =======
*
* ECHBMV 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.
* Additional precision arithmetic is used in the computation.
*
* 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 .
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - COMPLEX 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 hemitian 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 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 .
* 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.
* At least double precision arithmetic is used in the
* computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX ONE
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
COMPLEX ZERO
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+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 DCMPLX, MAX, MIN, CONJG, REAL
* ..
* .. 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( 'ECHBMV', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* 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 and set up the start points in X and Y if
* the increments are not both unity.
*
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
IF( BETA.NE.ONE )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
END IF
ELSE
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
IF( BETA.NE.ONE )THEN
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 = DCMPLX( 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 +
$ CONJG( A( L + I, J ) )*DCMPLX( X( I ) )
50 CONTINUE
Y( J ) = Y( J ) + TEMP1*REAL( A( KPLUS1, J ) )
$ + ALPHA*TEMP2
60 CONTINUE
ELSE
IX = KX - INCX
DO 80, J = 1, N
TEMP1 = DCMPLX( ALPHA )*X( IX + INCX )
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 +
$ CONJG( A( L + I, J ) )*DCMPLX( X( IX ) )
IX = IX + INCX
IY = IY + INCY
70 CONTINUE
Y( IY ) = Y( IY ) + TEMP1*REAL( A( KPLUS1, J ) )
$ + ALPHA*TEMP2
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 = DCMPLX( ALPHA )*X( J )
TEMP2 = ZERO
Y( J ) = Y( J ) + TEMP1*REAL( 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 +
$ CONJG( A( L + I, J ) )*DCMPLX( X( I ) )
90 CONTINUE
Y( J ) = Y( J ) + ALPHA*TEMP2
100 CONTINUE
ELSE
JX = KX
JY = KY
DO 120, J = 1, N
TEMP1 = DCMPLX( ALPHA )*X( JX )
TEMP2 = ZERO
Y( JY ) = Y( JY ) + TEMP1*REAL( 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 +
$ CONJG( A( L + I, J ) )*DCMPLX( 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 ECHBMV.
*
END
*
************************************************************************
*
SUBROUTINE ECTRMV( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, LDA, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
COMPLEX*16 X( * )
COMPLEX A( LDA, * )
* ..
*
* Purpose
* =======
*
* ECTRMV performs one of the matrix-vector operations
*
* x := A*x, or x := A'*x, or x := conjg( A' )*x,
*
* where x is n element vector and A is an n by n unit, or non-unit,
* upper or lower triangular matrix. Additional precision arithmetic is
* used in the computation.
*
* 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 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. At least double precision arithmetic is
* used in the computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
INTEGER I, INFO, IX, J, JX, KX
LOGICAL NOCONJ, NOUNIT
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC CONJG, 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( 'ECTRMV', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = ( DIAG .EQ.'N' ).OR.( DIAG .EQ.'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
DO 10, I = 1, J - 1
X( I ) = X( I ) + X( J )*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
IX = KX
DO 30, I = 1, J - 1
X( IX ) = X( IX ) + X( JX )*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
DO 50, I = N, J + 1, -1
X( I ) = X( I ) + X( J )*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
IX = KX
DO 70, I = N, J + 1, -1
X( IX ) = X( IX ) + X( JX )*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
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( J ) = X( J )*A( J, J )
DO 90, I = J - 1, 1, -1
X( J ) = X( J ) + A( I, J )*X( I )
90 CONTINUE
ELSE
IF( NOUNIT )
$ X( J ) = X( J )*CONJG( A( J, J ) )
DO 100, I = J - 1, 1, -1
X( J ) = X( J ) + CONJG( A( I, J ) )*X( I )
100 CONTINUE
END IF
110 CONTINUE
ELSE
JX = KX + ( N - 1 )*INCX
DO 140, J = N, 1, -1
IX = JX
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )*A( J, J )
DO 120, I = J - 1, 1, -1
IX = IX - INCX
X( JX ) = X( JX ) + A( I, J )*X( IX )
120 CONTINUE
ELSE
IF( NOUNIT )
$ X( JX ) = X( JX )*CONJG( A( J, J ) )
DO 130, I = J - 1, 1, -1
IX = IX - INCX
X( JX ) = X( JX ) + CONJG( A( I, J ) )*X( IX )
130 CONTINUE
END IF
JX = JX - INCX
140 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 170, J = 1, N
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( J ) = X( J )*A( J, J )
DO 150, I = J + 1, N
X( J ) = X( J ) + A( I, J )*X( I )
150 CONTINUE
ELSE
IF( NOUNIT )
$ X( J ) = X( J )*CONJG( A( J, J ) )
DO 160, I = J + 1, N
X( J ) = X( J ) + CONJG( A( I, J ) )*X( I )
160 CONTINUE
END IF
170 CONTINUE
ELSE
JX = KX
DO 200, J = 1, N
IX = JX
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )*A( J, J )
DO 180, I = J + 1, N
IX = IX + INCX
X( JX ) = X( JX ) + A( I, J )*X( IX )
180 CONTINUE
ELSE
IF( NOUNIT )
$ X( JX ) = X( JX )*CONJG( A( J, J ) )
DO 190, I = J + 1, N
IX = IX + INCX
X( JX ) = X( JX ) + CONJG( A( I, J ) )*X( IX )
190 CONTINUE
END IF
JX = JX + INCX
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ECTRMV.
*
END
*
************************************************************************
*
SUBROUTINE ECTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, LDA, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
COMPLEX*16 X( * )
COMPLEX A( LDA, * )
* ..
*
* Purpose
* =======
*
* ECTRSV 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. Additional precision
* arithmetic is used in the computation.
*
* 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 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. At least double precision
* arithmetic is used in the computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
INTEGER I, INFO, IX, J, JX, KX
LOGICAL NOCONJ, NOUNIT
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC CONJG, 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( 'ECTRSV', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = ( DIAG .EQ.'N' ).OR.( DIAG .EQ.'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 )
DO 10, I = J - 1, 1, -1
X( I ) = X( I ) - X( J )*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 )
IX = JX
DO 30, I = J - 1, 1, -1
IX = IX - INCX
X( IX ) = X( IX ) - X( JX )*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 )
DO 50, I = J + 1, N
X( I ) = X( I ) - X( J )*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 )
IX = JX
DO 70, I = J + 1, N
IX = IX + INCX
X( IX ) = X( IX ) - X( JX )*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
IF( NOCONJ )THEN
DO 90, I = 1, J - 1
X( J ) = X( J ) - A( I, J )*X( I )
90 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/A( J, J )
ELSE
DO 100, I = 1, J - 1
X( J ) = X( J ) - CONJG( A( I, J ) )*X( I )
100 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/CONJG( A( J, J ) )
END IF
110 CONTINUE
ELSE
JX = KX
DO 140, J = 1, N
IX = KX
IF( NOCONJ )THEN
DO 120, I = 1, J - 1
X( JX ) = X( JX ) - A( I, J )*X( IX )
IX = IX + INCX
120 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/A( J, J )
ELSE
DO 130, I = 1, J - 1
X( JX ) = X( JX ) - CONJG( A( I, J ) )*X( IX )
IX = IX + INCX
130 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/CONJG( A( J, J ) )
END IF
JX = JX + INCX
140 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 170, J = N, 1, -1
IF( NOCONJ )THEN
DO 150, I = N, J + 1, -1
X( J ) = X( J ) - A( I, J )*X( I )
150 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/A( J, J )
ELSE
DO 160, I = N, J + 1, -1
X( J ) = X( J ) - CONJG( A( I, J ) )*X( I )
160 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/CONJG( A( J, J ) )
END IF
170 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 200, J = N, 1, -1
IX = KX
IF( NOCONJ )THEN
DO 180, I = N, J + 1, -1
X( JX ) = X( JX ) - A( I, J )*X( IX )
IX = IX - INCX
180 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/A( J, J )
ELSE
DO 190, I = N, J + 1, -1
X( JX ) = X( JX ) - CONJG( A( I, J ) )*X( IX )
IX = IX - INCX
190 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/CONJG( A( J, J ) )
END IF
JX = JX - INCX
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ECTRSV.
*
END
*
************************************************************************
*
SUBROUTINE ECTPMV( UPLO, TRANS, DIAG, N, AP, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
COMPLEX*16 X( * )
COMPLEX AP( * )
* ..
*
* Purpose
* =======
*
* ECTPMV performs one of the matrix-vector operations
*
* x := A*x, or x := A'*x, or x := conjg( A' )*x,
*
* where x is n element vector and A is an n by n unit, or non-unit,
* upper or lower triangular matrix. Additional precision arithmetic is
* used in the computation.
*
* 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 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. At least double precision arithmetic is
* used in the computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
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 CONJG
* ..
* .. 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( 'ECTPMV', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = ( DIAG .EQ.'N' ).OR.( DIAG .EQ.'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
K = 1
IF( INCX.EQ.1 )THEN
DO 20, J = 1, N
IF( X( J ).NE.ZERO )THEN
DO 10, I = 1, J - 1
X( I ) = X( I ) + X( J )*AP( K )
K = K + 1
10 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )*AP( K )
K = K + 1
ELSE
K = K + J
END IF
20 CONTINUE
ELSE
JX = KX
DO 40, J = 1, N
IF( X( JX ).NE.ZERO )THEN
IX = KX
KK = K
DO 30, K = KK, KK + J - 2
X( IX ) = X( IX ) + X( JX )*AP( K )
IX = IX + INCX
30 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )*AP( K )
K = K + 1
ELSE
K = K + J
END IF
JX = JX + INCX
40 CONTINUE
END IF
ELSE
K = ( N*( N + 1 ) )/2
IF( INCX.EQ.1 )THEN
DO 60, J = N, 1, -1
IF( X( J ).NE.ZERO )THEN
DO 50, I = N, J + 1, -1
X( I ) = X( I ) + X( J )*AP( K )
K = K - 1
50 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )*AP( K )
K = K - 1
ELSE
K = K - ( N - J + 1 )
END IF
60 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 80, J = N, 1, -1
IF( X( JX ).NE.ZERO )THEN
IX = KX
KK = K
DO 70, K = KK, KK - ( N - ( J + 1 ) ), -1
X( IX ) = X( IX ) + X( JX )*AP( K )
IX = IX - INCX
70 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )*AP( K )
K = K - 1
ELSE
K = K - ( N - J + 1 )
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
K = ( N*( N + 1 ) )/2
IF( INCX.EQ.1 )THEN
DO 110, J = N, 1, -1
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( J ) = X( J )*AP( K )
K = K - 1
DO 90, I = J - 1, 1, -1
X( J ) = X( J ) + AP( K )*X( I )
K = K - 1
90 CONTINUE
ELSE
IF( NOUNIT )
$ X( J ) = X( J )*CONJG( AP( K ) )
K = K - 1
DO 100, I = J - 1, 1, -1
X( J ) = X( J ) + CONJG( AP( K ) )*X( I )
K = K - 1
100 CONTINUE
END IF
110 CONTINUE
ELSE
JX = KX + ( N - 1 )*INCX
DO 140, J = N, 1, -1
IX = JX
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )*AP( K )
KK = K - 1
DO 120, K = KK, KK - J + 2, -1
IX = IX - INCX
X( JX ) = X( JX ) + AP( K )*X( IX )
120 CONTINUE
ELSE
IF( NOUNIT )
$ X( JX ) = X( JX )*CONJG( AP( K ) )
KK = K - 1
DO 130, K = KK, KK - J + 2, -1
IX = IX - INCX
X( JX ) = X( JX ) + CONJG( AP( K ) )*X( IX )
130 CONTINUE
END IF
JX = JX - INCX
140 CONTINUE
END IF
ELSE
K = 1
IF( INCX.EQ.1 )THEN
DO 170, J = 1, N
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( J ) = X( J )*AP( K )
K = K + 1
DO 150, I = J + 1, N
X( J ) = X( J ) + AP( K )*X( I )
K = K + 1
150 CONTINUE
ELSE
IF( NOUNIT )
$ X( J ) = X( J )*CONJG( AP( K ) )
K = K + 1
DO 160, I = J + 1, N
X( J ) = X( J ) + CONJG( AP( K ) )*X( I )
K = K + 1
160 CONTINUE
END IF
170 CONTINUE
ELSE
JX = KX
DO 200, J = 1, N
IX = JX
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )*AP( K )
KK = K + 1
DO 180, K = KK, KK + N - ( J + 1 )
IX = IX + INCX
X( JX ) = X( JX ) + AP( K )*X( IX )
180 CONTINUE
ELSE
IF( NOUNIT )
$ X( JX ) = X( JX )*CONJG( AP( K ) )
KK = K + 1
DO 190, K = KK, KK + N - ( J + 1 )
IX = IX + INCX
X( JX ) = X( JX ) + CONJG( AP( K ) )*X( IX )
190 CONTINUE
END IF
JX = JX + INCX
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ECTPMV.
*
END
*
************************************************************************
*
SUBROUTINE ECTPSV( UPLO, TRANS, DIAG, N, AP, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
COMPLEX*16 X( * )
COMPLEX AP( * )
* ..
*
* Purpose
* =======
*
* ECTPSV 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. Additional precision
* arithmetic is used in the computation.
*
* 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 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. At least double precision
* arithemtic is used in the computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
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 CONJG
* ..
* .. 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( 'ECTPSV', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = ( DIAG .EQ.'N' ).OR.( DIAG .EQ.'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
K = ( 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( K )
K = K - 1
DO 10, I = J - 1, 1, -1
X( I ) = X( I ) - X( J )*AP( K )
K = K - 1
10 CONTINUE
ELSE
K = K - J
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 )/AP( K )
IX = JX
KK = K - 1
DO 30, K = KK, KK - J + 2, -1
IX = IX - INCX
X( IX ) = X( IX ) - X( JX )*AP( K )
30 CONTINUE
ELSE
K = K - J
END IF
JX = JX - INCX
40 CONTINUE
END IF
ELSE
K = 1
IF( INCX.EQ.1 )THEN
DO 60, J = 1, N
IF( X( J ).NE.ZERO )THEN
IF( NOUNIT )
$ X( J ) = X( J )/AP( K )
K = K + 1
DO 50, I = J + 1, N
X( I ) = X( I ) - X( J )*AP( K )
K = K + 1
50 CONTINUE
ELSE
K = K + N - J + 1
END IF
60 CONTINUE
ELSE
JX = KX
DO 80, J = 1, N
IF( X( JX ).NE.ZERO )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )/AP( K )
IX = JX
KK = K + 1
DO 70, K = KK, KK + N - ( J + 1 )
IX = IX + INCX
X( IX ) = X( IX ) - X( JX )*AP( K )
70 CONTINUE
ELSE
K = K + N - J + 1
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
K = 1
IF( INCX.EQ.1 )THEN
DO 110, J = 1, N
IF( NOCONJ )THEN
DO 90, I = 1, J - 1
X( J ) = X( J ) - AP( K )*X( I )
K = K + 1
90 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/AP( K )
ELSE
DO 100, I = 1, J - 1
X( J ) = X( J ) - CONJG( AP( K ) )*X( I )
K = K + 1
100 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/CONJG( AP( K ) )
END IF
K = K + 1
110 CONTINUE
ELSE
JX = KX
DO 140, J = 1, N
IX = KX
KK = K
IF( NOCONJ )THEN
DO 120, K = KK, KK + J - 2
X( JX ) = X( JX ) - AP( K )*X( IX )
IX = IX + INCX
120 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/AP( K )
ELSE
DO 130, K = KK, KK + J - 2
X( JX ) = X( JX ) - CONJG( AP( K ) )*X( IX )
IX = IX + INCX
130 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/CONJG( AP( K ) )
END IF
K = K + 1
JX = JX + INCX
140 CONTINUE
END IF
ELSE
K = ( N*( N + 1 ) )/2
IF( INCX.EQ.1 )THEN
DO 170, J = N, 1, -1
IF( NOCONJ )THEN
DO 150, I = N, J + 1, -1
X( J ) = X( J ) - AP( K )*X( I )
K = K - 1
150 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/AP( K )
ELSE
DO 160, I = N, J + 1, -1
X( J ) = X( J ) - CONJG( AP( K ) )*X( I )
K = K - 1
160 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/CONJG( AP( K ) )
END IF
K = K - 1
170 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 200, J = N, 1, -1
IX = KX
KK = K
IF( NOCONJ )THEN
DO 180, K = KK, KK - ( N - ( J + 1 ) ), -1
X( JX ) = X( JX ) - AP( K )*X( IX )
IX = IX - INCX
180 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/AP( K )
ELSE
DO 190, K = KK, KK - ( N - ( J + 1 ) ), -1
X( JX ) = X( JX ) - CONJG( AP( K ) )*X( IX )
IX = IX - INCX
190 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/CONJG( AP( K ) )
END IF
K = K - 1
JX = JX - INCX
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ECTPSV.
*
END
*
************************************************************************
*
SUBROUTINE ECTBMV( 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 X( * )
COMPLEX A( LDA, * )
* ..
*
* Purpose
* =======
*
* ECTBMV performs one of the matrix-vector operations
*
* x := A*x, or x := A'*x, or x := conjg( A' )*x,
*
* where x is n element vector and A is an n by n unit, or non-unit,
* upper or lower triangular band matrix, with ( k + 1 ) diagonals.
* Additional precision arithmetic is used in the computation.
*
* 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 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. At least double precision arithmetic is
* used in the computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
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 MAX, MIN, CONJG
* ..
* .. 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( 'ECTBMV', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = ( DIAG .EQ.'N' ).OR.( DIAG .EQ.'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
L = KPLUS1 - J
DO 10, I = MAX( 1, J - K ), J - 1
X( I ) = X( I ) + X( J )*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
IX = KX
L = KPLUS1 - J
DO 30, I = MAX( 1, J - K ), J - 1
X( IX ) = X( IX ) + X( JX )*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
L = 1 - J
DO 50, I = MIN( N, J + K ), J + 1, -1
X( I ) = X( I ) + X( J )*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
IX = KX
L = 1 - J
DO 70, I = MIN( N, J + K ), J + 1, -1
X( IX ) = X( IX ) + X( JX )*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
L = KPLUS1 - J
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( J ) = X( J )*A( KPLUS1, J )
DO 90, I = J - 1, MAX( 1, J - K ), -1
X( J ) = X( J ) + A( L + I, J )*X( I )
90 CONTINUE
ELSE
IF( NOUNIT )
$ X( J ) = X( J )*CONJG( A( KPLUS1, J ) )
DO 100, I = J - 1, MAX( 1, J - K ), -1
X( J ) = X( J ) + CONJG( A( L + I, J ) )*X( I )
100 CONTINUE
END IF
110 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 140, J = N, 1, -1
KX = KX - INCX
IX = KX
L = KPLUS1 - J
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )*A( KPLUS1, J )
DO 120, I = J - 1, MAX( 1, J - K ), -1
X( JX ) = X( JX ) + A( L + I, J )*X( IX )
IX = IX - INCX
120 CONTINUE
ELSE
IF( NOUNIT )
$ X( JX ) = X( JX )*CONJG( A( KPLUS1, J ) )
DO 130, I = J - 1, MAX( 1, J - K ), -1
X( JX ) = X( JX ) +
$ CONJG( A( L + I, J ) )*X( IX )
IX = IX - INCX
130 CONTINUE
END IF
JX = JX - INCX
140 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 170, J = 1, N
L = 1 - J
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( J ) = X( J )*A( 1, J )
DO 150, I = J + 1, MIN( N, J + K )
X( J ) = X( J ) + A( L + I, J )*X( I )
150 CONTINUE
ELSE
IF( NOUNIT )
$ X( J ) = X( J )*CONJG( A( 1, J ) )
DO 160, I = J + 1, MIN( N, J + K )
X( J ) = X( J ) + CONJG( A( L + I, J ) )*X( I )
160 CONTINUE
END IF
170 CONTINUE
ELSE
JX = KX
DO 200, J = 1, N
KX = KX + INCX
IX = KX
L = 1 - J
IF( NOCONJ )THEN
IF( NOUNIT )
$ X( JX ) = X( JX )*A( 1, J )
DO 180, I = J + 1, MIN( N, J + K )
X( JX ) = X( JX ) + A( L + I, J )*X( IX )
IX = IX + INCX
180 CONTINUE
ELSE
IF( NOUNIT )
$ X( JX ) = X( JX )*CONJG( A( 1, J ) )
DO 190, I = J + 1, MIN( N, J + K )
X( JX ) = X( JX ) +
$ CONJG( A( L + I, J ) )*X( IX )
IX = IX + INCX
190 CONTINUE
END IF
JX = JX + INCX
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ECTBMV.
*
END
*
************************************************************************
*
SUBROUTINE ECTBSV( 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 X( * )
COMPLEX A( LDA, * )
* ..
*
* Purpose
* =======
*
* ECTBSV 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. Additional precision arithmetic is used in the
* computation.
*
* 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 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. At least double precision
* arithmetic is used in the computation of 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 20-July-1986.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
COMPLEX*16 ZERO
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
* .. Local Scalars ..
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 MAX, MIN, CONJG
* ..
* .. 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( 'ECTBSV', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOCONJ = LSAME( TRANS, 'T' )
NOUNIT = ( DIAG .EQ.'N' ).OR.( DIAG .EQ.'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 )
DO 10, I = J - 1, MAX( 1, J - K ), -1
X( I ) = X( I ) - X( J )*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 )
DO 30, I = J - 1, MAX( 1, J - K ), -1
X( IX ) = X( IX ) - X( JX )*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 )
DO 50, I = J + 1, MIN( N, J + K )
X( I ) = X( I ) - X( J )*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 )
DO 70, I = J + 1, MIN( N, J + K )
X( IX ) = X( IX ) - X( JX )*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
L = KPLUS1 - J
IF( NOCONJ )THEN
DO 90, I = MAX( 1, J - K ), J - 1
X( J ) = X( J ) - A( L + I, J )*X( I )
90 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/A( KPLUS1, J )
ELSE
DO 100, I = MAX( 1, J - K ), J - 1
X( J ) = X( J ) - CONJG( A( L + I, J ) )*X( I )
100 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/CONJG( A( KPLUS1, J ) )
END IF
110 CONTINUE
ELSE
JX = KX
DO 140, J = 1, N
IX = KX
L = KPLUS1 - J
IF( NOCONJ )THEN
DO 120, I = MAX( 1, J - K ), J - 1
X( JX ) = X( JX ) - A( L + I, J )*X( IX )
IX = IX + INCX
120 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/A( KPLUS1, J )
ELSE
DO 130, I = MAX( 1, J - K ), J - 1
X( JX ) = X( JX ) -
$ CONJG( A( L + I, J ) )*X( IX )
IX = IX + INCX
130 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/CONJG( A( KPLUS1, J ) )
END IF
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
L = 1 - J
IF( NOCONJ )THEN
DO 150, I = MIN( N, J + K ), J + 1, -1
X( J ) = X( J ) - A( L + I, J )*X( I )
150 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/A( 1, J )
ELSE
DO 160, I = MIN( N, J + K ), J + 1, -1
X( J ) = X( J ) - CONJG( A( L + I, J ) )*X( I )
160 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )/CONJG( A( 1, J ) )
END IF
170 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 200, J = N, 1, -1
IX = KX
L = 1 - J
IF( NOCONJ )THEN
DO 180, I = MIN( N, J + K ), J + 1, -1
X( JX ) = X( JX ) - A( L + I, J )*X( IX )
IX = IX - INCX
180 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/A( 1, J )
ELSE
DO 190, I = MIN( N, J + K ), J + 1, -1
X( JX ) = X( JX ) -
$ CONJG( A( L + I, J ) )*X( IX )
IX = IX - INCX
190 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )/CONJG( A( 1, J ) )
END IF
JX = JX - INCX
IF( ( N - J ).GE.K )
$ KX = KX - INCX
200 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of ECTBSV.
*
END