* ************************************************************************ * * 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