SUBROUTINE ZGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) !! ZGEMM performs C:=alpha*A*B+beta*C, A, B, C rectangular. ! ! .. Scalar Arguments .. DOUBLE COMPLEX ALPHA,BETA INTEGER K,LDA,LDB,LDC,M,N CHARACTER TRANSA,TRANSB ! .. ! .. Array Arguments .. DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) ! .. ! ! Purpose ! ======= ! ! ZGEMM performs one of the matrix-matrix operations ! ! C := alpha*op( A )*op( B ) + beta*C, ! ! where op( X ) is one of ! ! op( X ) = X or op( X ) = X**T or op( X ) = X**H, ! ! alpha and beta are scalars, and A, B and C are matrices, with op( A ) ! an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. ! ! Arguments ! ========== ! ! TRANSA - CHARACTER*1. ! On entry, TRANSA specifies the form of op( A ) to be used in ! the matrix multiplication as follows: ! ! TRANSA = 'N' or 'n', op( A ) = A. ! ! TRANSA = 'T' or 't', op( A ) = A**T. ! ! TRANSA = 'C' or 'c', op( A ) = A**H. ! ! Unchanged on exit. ! ! TRANSB - CHARACTER*1. ! On entry, TRANSB specifies the form of op( B ) to be used in ! the matrix multiplication as follows: ! ! TRANSB = 'N' or 'n', op( B ) = B. ! ! TRANSB = 'T' or 't', op( B ) = B**T. ! ! TRANSB = 'C' or 'c', op( B ) = B**H. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix ! op( A ) and of the matrix C. M must be at least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix ! op( B ) and the number of columns of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry, K specifies the number of columns of the matrix ! op( A ) and the number of rows of the matrix op( B ). K must ! be at least zero. ! Unchanged on exit. ! ! ALPHA - COMPLEX*16 . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is ! k when TRANSA = 'N' or 'n', and is m otherwise. ! Before entry with TRANSA = 'N' or 'n', the leading m by k ! part of the array A must contain the matrix A, otherwise ! the leading k by m part of the array A must contain the ! matrix A. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When TRANSA = 'N' or 'n' then ! LDA must be at least max( 1, m ), otherwise LDA must be at ! least max( 1, k ). ! Unchanged on exit. ! ! B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is ! n when TRANSB = 'N' or 'n', and is k otherwise. ! Before entry with TRANSB = 'N' or 'n', the leading k by n ! part of the array B must contain the matrix B, otherwise ! the leading n by k part of the array B must contain the ! matrix B. ! Unchanged on exit. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. When TRANSB = 'N' or 'n' then ! LDB must be at least max( 1, k ), otherwise LDB must be at ! least max( 1, n ). ! Unchanged on exit. ! ! BETA - COMPLEX*16 . ! On entry, BETA specifies the scalar beta. When BETA is ! supplied as zero then C need not be set on input. ! Unchanged on exit. ! ! Input/output, COMPLEX*16 C(LDC,N). ! Before entry, the leading m by n part of the array C must ! contain the matrix C, except when beta is zero, in which ! case C need not be set. ! On exit, the array C is overwritten by the m by n matrix ! ( alpha*op( A )*op( B ) + beta*C ). ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, m ). ! Unchanged on exit. ! ! Further Details ! =============== ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ===================================================================== ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC DCONJG,MAX ! .. ! .. Local Scalars .. DOUBLE COMPLEX TEMP INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB LOGICAL CONJA,CONJB,NOTA,NOTB ! .. ! .. Parameters .. DOUBLE COMPLEX ONE PARAMETER (ONE= (1.0D+0,0.0D+0)) DOUBLE COMPLEX ZERO PARAMETER (ZERO= (0.0D+0,0.0D+0)) ! .. ! ! Set NOTA and NOTB as true if A and B respectively are not ! conjugated or transposed, set CONJA and CONJB as true if A and ! B respectively are to be transposed but not conjugated and set ! NROWA, NCOLA and NROWB as the number of rows and columns of A ! and the number of rows of B respectively. ! NOTA = LSAME(TRANSA,'N') NOTB = LSAME(TRANSB,'N') CONJA = LSAME(TRANSA,'C') CONJB = LSAME(TRANSB,'C') IF (NOTA) THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF (NOTB) THEN NROWB = K ELSE NROWB = N END IF ! ! Test the input parameters. ! INFO = 0 IF ((.NOT.NOTA) .AND. (.NOT.CONJA) .AND. & (.NOT.LSAME(TRANSA,'T'))) THEN INFO = 1 ELSE IF ((.NOT.NOTB) .AND. (.NOT.CONJB) .AND. & (.NOT.LSAME(TRANSB,'T'))) THEN INFO = 2 ELSE IF (M.LT.0) THEN INFO = 3 ELSE IF (N.LT.0) THEN INFO = 4 ELSE IF (K.LT.0) THEN INFO = 5 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN INFO = 8 ELSE IF (LDB.LT.MAX(1,NROWB)) THEN INFO = 10 ELSE IF (LDC.LT.MAX(1,M)) THEN INFO = 13 END IF IF (INFO.NE.0) THEN CALL XERBLA('ZGEMM ',INFO) RETURN END IF ! ! Quick return if possible. ! IF ((M.EQ.0) .OR. (N.EQ.0) .OR. & (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN ! ! And when alpha.eq.zero. ! IF (ALPHA.EQ.ZERO) THEN IF (BETA.EQ.ZERO) THEN DO 20 J = 1,N DO 10 I = 1,M C(I,J) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40 J = 1,N DO 30 I = 1,M C(I,J) = BETA*C(I,J) 30 CONTINUE 40 CONTINUE END IF RETURN END IF ! ! Start the operations. ! IF (NOTB) THEN IF (NOTA) THEN ! ! Form C := alpha*A*B + beta*C. ! DO 90 J = 1,N IF (BETA.EQ.ZERO) THEN DO 50 I = 1,M C(I,J) = ZERO 50 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 60 I = 1,M C(I,J) = BETA*C(I,J) 60 CONTINUE END IF DO 80 L = 1,K IF (B(L,J).NE.ZERO) THEN TEMP = ALPHA*B(L,J) DO 70 I = 1,M C(I,J) = C(I,J) + TEMP*A(I,L) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE IF (CONJA) THEN ! ! Form C := alpha*A**H*B + beta*C. ! DO 120 J = 1,N DO 110 I = 1,M TEMP = ZERO DO 100 L = 1,K TEMP = TEMP + DCONJG(A(L,I))*B(L,J) 100 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP ELSE C(I,J) = ALPHA*TEMP + BETA*C(I,J) END IF 110 CONTINUE 120 CONTINUE ELSE ! ! Form C := alpha*A**T*B + beta*C ! DO 150 J = 1,N DO 140 I = 1,M TEMP = ZERO DO 130 L = 1,K TEMP = TEMP + A(L,I)*B(L,J) 130 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP ELSE C(I,J) = ALPHA*TEMP + BETA*C(I,J) END IF 140 CONTINUE 150 CONTINUE END IF ELSE IF (NOTA) THEN IF (CONJB) THEN ! ! Form C := alpha*A*B**H + beta*C. ! DO 200 J = 1,N IF (BETA.EQ.ZERO) THEN DO 160 I = 1,M C(I,J) = ZERO 160 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 170 I = 1,M C(I,J) = BETA*C(I,J) 170 CONTINUE END IF DO 190 L = 1,K IF (B(J,L).NE.ZERO) THEN TEMP = ALPHA*DCONJG(B(J,L)) DO 180 I = 1,M C(I,J) = C(I,J) + TEMP*A(I,L) 180 CONTINUE END IF 190 CONTINUE 200 CONTINUE ELSE ! ! Form C := alpha*A*B**T + beta*C ! DO 250 J = 1,N IF (BETA.EQ.ZERO) THEN DO 210 I = 1,M C(I,J) = ZERO 210 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 220 I = 1,M C(I,J) = BETA*C(I,J) 220 CONTINUE END IF DO 240 L = 1,K IF (B(J,L).NE.ZERO) THEN TEMP = ALPHA*B(J,L) DO 230 I = 1,M C(I,J) = C(I,J) + TEMP*A(I,L) 230 CONTINUE END IF 240 CONTINUE 250 CONTINUE END IF ELSE IF (CONJA) THEN IF (CONJB) THEN ! ! Form C := alpha*A**H*B**H + beta*C. ! DO 280 J = 1,N DO 270 I = 1,M TEMP = ZERO DO 260 L = 1,K TEMP = TEMP + DCONJG(A(L,I))*DCONJG(B(J,L)) 260 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP ELSE C(I,J) = ALPHA*TEMP + BETA*C(I,J) END IF 270 CONTINUE 280 CONTINUE ELSE ! ! Form C := alpha*A**H*B**T + beta*C ! DO 310 J = 1,N DO 300 I = 1,M TEMP = ZERO DO 290 L = 1,K TEMP = TEMP + DCONJG(A(L,I))*B(J,L) 290 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP ELSE C(I,J) = ALPHA*TEMP + BETA*C(I,J) END IF 300 CONTINUE 310 CONTINUE END IF ELSE IF (CONJB) THEN ! ! Form C := alpha*A**T*B**H + beta*C ! DO 340 J = 1,N DO 330 I = 1,M TEMP = ZERO DO 320 L = 1,K TEMP = TEMP + A(L,I)*DCONJG(B(J,L)) 320 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP ELSE C(I,J) = ALPHA*TEMP + BETA*C(I,J) END IF 330 CONTINUE 340 CONTINUE ELSE ! ! Form C := alpha*A**T*B**T + beta*C ! DO 370 J = 1,N DO 360 I = 1,M TEMP = ZERO DO 350 L = 1,K TEMP = TEMP + A(L,I)*B(J,L) 350 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP ELSE C(I,J) = ALPHA*TEMP + BETA*C(I,J) END IF 360 CONTINUE 370 CONTINUE END IF END IF ! RETURN ! ! End of ZGEMM . ! END SUBROUTINE ZHEMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC) !! ZHEMM performs C:= alpha*A*B+beta*C, for A hermitian. ! ! .. Scalar Arguments .. DOUBLE COMPLEX ALPHA,BETA INTEGER LDA,LDB,LDC,M,N CHARACTER SIDE,UPLO ! .. ! .. Array Arguments .. DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) ! .. ! ! Purpose ! ======= ! ! ZHEMM performs one of the matrix-matrix operations ! ! C := alpha*A*B + beta*C, ! ! or ! ! C := alpha*B*A + beta*C, ! ! where alpha and beta are scalars, A is an hermitian matrix and B and ! C are m by n matrices. ! ! Arguments ! ========== ! ! SIDE - CHARACTER*1. ! On entry, SIDE specifies whether the hermitian matrix A ! appears on the left or right in the operation as follows: ! ! SIDE = 'L' or 'l' C := alpha*A*B + beta*C, ! ! SIDE = 'R' or 'r' C := alpha*B*A + beta*C, ! ! Unchanged on exit. ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the hermitian matrix A is to be ! referenced as follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of the ! hermitian matrix is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of the ! hermitian matrix is to be referenced. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix C. ! M must be at least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix C. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - COMPLEX*16 . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is ! m when SIDE = 'L' or 'l' and is n otherwise. ! Before entry with SIDE = 'L' or 'l', the m by m part of ! the array A must contain the hermitian matrix, such that ! when UPLO = 'U' or 'u', the leading m by m upper triangular ! part of the array A must contain the upper triangular part ! of the hermitian matrix and the strictly lower triangular ! part of A is not referenced, and when UPLO = 'L' or 'l', ! the leading m by m lower triangular part of the array A ! must contain the lower triangular part of the hermitian ! matrix and the strictly upper triangular part of A is not ! referenced. ! Before entry with SIDE = 'R' or 'r', the n by n part of ! the array A must contain the hermitian matrix, such that ! when UPLO = 'U' or 'u', the leading n by n upper triangular ! part of the array A must contain the upper triangular part ! of the hermitian matrix and the strictly lower triangular ! part of A is not referenced, and when UPLO = 'L' or 'l', ! the leading n by n lower triangular part of the array A ! must contain the lower triangular part of the hermitian ! matrix and the strictly upper triangular part of A is not ! referenced. ! Note that the imaginary parts of the diagonal elements need ! not be set, they are assumed to be zero. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When SIDE = 'L' or 'l' then ! LDA must be at least max( 1, m ), otherwise LDA must be at ! least max( 1, n ). ! Unchanged on exit. ! ! B - COMPLEX*16 array of DIMENSION ( LDB, n ). ! Before entry, the leading m by n part of the array B must ! contain the matrix B. ! Unchanged on exit. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. LDB must be at least ! max( 1, m ). ! Unchanged on exit. ! ! BETA - COMPLEX*16 . ! On entry, BETA specifies the scalar beta. When BETA is ! supplied as zero then C need not be set on input. ! Unchanged on exit. ! ! C - COMPLEX*16 array of DIMENSION ( LDC, n ). ! Before entry, the leading m by n part of the array C must ! contain the matrix C, except when beta is zero, in which ! case C need not be set on entry. ! On exit, the array C is overwritten by the m by n updated ! matrix. ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, m ). ! Unchanged on exit. ! ! Further Details ! =============== ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ===================================================================== ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC DBLE,DCONJG,MAX ! .. ! .. Local Scalars .. DOUBLE COMPLEX TEMP1,TEMP2 INTEGER I,INFO,J,K,NROWA LOGICAL UPPER ! .. ! .. Parameters .. DOUBLE COMPLEX ONE PARAMETER (ONE= (1.0D+0,0.0D+0)) DOUBLE COMPLEX ZERO PARAMETER (ZERO= (0.0D+0,0.0D+0)) ! .. ! ! Set NROWA as the number of rows of A. ! IF (LSAME(SIDE,'L')) THEN NROWA = M ELSE NROWA = N END IF UPPER = LSAME(UPLO,'U') ! ! Test the input parameters. ! INFO = 0 IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN INFO = 1 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN INFO = 2 ELSE IF (M.LT.0) THEN INFO = 3 ELSE IF (N.LT.0) THEN INFO = 4 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN INFO = 7 ELSE IF (LDB.LT.MAX(1,M)) THEN INFO = 9 ELSE IF (LDC.LT.MAX(1,M)) THEN INFO = 12 END IF IF (INFO.NE.0) THEN CALL XERBLA('ZHEMM ',INFO) RETURN END IF ! ! Quick return if possible. ! IF ((M.EQ.0) .OR. (N.EQ.0) .OR. & ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN ! ! And when alpha.eq.zero. ! IF (ALPHA.EQ.ZERO) THEN IF (BETA.EQ.ZERO) THEN DO 20 J = 1,N DO 10 I = 1,M C(I,J) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40 J = 1,N DO 30 I = 1,M C(I,J) = BETA*C(I,J) 30 CONTINUE 40 CONTINUE END IF RETURN END IF ! ! Start the operations. ! IF (LSAME(SIDE,'L')) THEN ! ! Form C := alpha*A*B + beta*C. ! IF (UPPER) THEN DO 70 J = 1,N DO 60 I = 1,M TEMP1 = ALPHA*B(I,J) TEMP2 = ZERO DO 50 K = 1,I - 1 C(K,J) = C(K,J) + TEMP1*A(K,I) TEMP2 = TEMP2 + B(K,J)*DCONJG(A(K,I)) 50 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = TEMP1*DBLE(A(I,I)) + ALPHA*TEMP2 ELSE C(I,J) = BETA*C(I,J) + TEMP1*DBLE(A(I,I)) + & ALPHA*TEMP2 END IF 60 CONTINUE 70 CONTINUE ELSE DO 100 J = 1,N DO 90 I = M,1,-1 TEMP1 = ALPHA*B(I,J) TEMP2 = ZERO DO 80 K = I + 1,M C(K,J) = C(K,J) + TEMP1*A(K,I) TEMP2 = TEMP2 + B(K,J)*DCONJG(A(K,I)) 80 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = TEMP1*DBLE(A(I,I)) + ALPHA*TEMP2 ELSE C(I,J) = BETA*C(I,J) + TEMP1*DBLE(A(I,I)) + & ALPHA*TEMP2 END IF 90 CONTINUE 100 CONTINUE END IF ELSE ! ! Form C := alpha*B*A + beta*C. ! DO 170 J = 1,N TEMP1 = ALPHA*DBLE(A(J,J)) IF (BETA.EQ.ZERO) THEN DO 110 I = 1,M C(I,J) = TEMP1*B(I,J) 110 CONTINUE ELSE DO 120 I = 1,M C(I,J) = BETA*C(I,J) + TEMP1*B(I,J) 120 CONTINUE END IF DO 140 K = 1,J - 1 IF (UPPER) THEN TEMP1 = ALPHA*A(K,J) ELSE TEMP1 = ALPHA*DCONJG(A(J,K)) END IF DO 130 I = 1,M C(I,J) = C(I,J) + TEMP1*B(I,K) 130 CONTINUE 140 CONTINUE DO 160 K = J + 1,N IF (UPPER) THEN TEMP1 = ALPHA*DCONJG(A(J,K)) ELSE TEMP1 = ALPHA*A(K,J) END IF DO 150 I = 1,M C(I,J) = C(I,J) + TEMP1*B(I,K) 150 CONTINUE 160 CONTINUE 170 CONTINUE END IF ! RETURN ! ! End of ZHEMM . ! END SUBROUTINE ZHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) !! ZHER2K performs C := alpha*A*conjg(B')+conjg(alpha)*b*conjg(A')+beta*C, for C hermitian. ! ! .. Scalar Arguments .. DOUBLE COMPLEX ALPHA DOUBLE PRECISION BETA INTEGER K,LDA,LDB,LDC,N CHARACTER TRANS,UPLO ! .. ! .. Array Arguments .. DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) ! .. ! ! Purpose ! ======= ! ! ZHER2K performs one of the hermitian rank 2k operations ! ! C := alpha*A*B**H + conjg( alpha )*B*A**H + beta*C, ! ! or ! ! C := alpha*A**H*B + conjg( alpha )*B**H*A + beta*C, ! ! where alpha and beta are scalars with beta real, C is an n by n ! hermitian matrix and A and B are n by k matrices in the first case ! and k by n matrices in the second case. ! ! Arguments ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array C is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of C ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of C ! is to be referenced. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' C := alpha*A*B**H + ! conjg( alpha )*B*A**H + ! beta*C. ! ! TRANS = 'C' or 'c' C := alpha*A**H*B + ! conjg( alpha )*B**H*A + ! beta*C. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry with TRANS = 'N' or 'n', K specifies the number ! of columns of the matrices A and B, and on entry with ! TRANS = 'C' or 'c', K specifies the number of rows of the ! matrices A and B. K must be at least zero. ! Unchanged on exit. ! ! ALPHA - COMPLEX*16 . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is ! k when TRANS = 'N' or 'n', and is n otherwise. ! Before entry with TRANS = 'N' or 'n', the leading n by k ! part of the array A must contain the matrix A, otherwise ! the leading k by n part of the array A must contain the ! matrix A. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When TRANS = 'N' or 'n' ! then LDA must be at least max( 1, n ), otherwise LDA must ! be at least max( 1, k ). ! Unchanged on exit. ! ! B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is ! k when TRANS = 'N' or 'n', and is n otherwise. ! Before entry with TRANS = 'N' or 'n', the leading n by k ! part of the array B must contain the matrix B, otherwise ! the leading k by n part of the array B must contain the ! matrix B. ! Unchanged on exit. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. When TRANS = 'N' or 'n' ! then LDB must be at least max( 1, n ), otherwise LDB must ! be at least max( 1, k ). ! Unchanged on exit. ! ! BETA - DOUBLE PRECISION . ! On entry, BETA specifies the scalar beta. ! Unchanged on exit. ! ! C - COMPLEX*16 array of DIMENSION ( LDC, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array C must contain the upper ! triangular part of the hermitian matrix and the strictly ! lower triangular part of C is not referenced. On exit, the ! upper triangular part of the array C is overwritten by the ! upper triangular part of the updated matrix. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array C must contain the lower ! triangular part of the hermitian matrix and the strictly ! upper triangular part of C is not referenced. On exit, the ! lower triangular part of the array C is overwritten by the ! lower triangular part of the updated matrix. ! Note that the imaginary parts of the diagonal elements need ! not be set, they are assumed to be zero, and on exit they ! are set to zero. ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, n ). ! Unchanged on exit. ! ! Further Details ! =============== ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1. ! Ed Anderson, Cray Research Inc. ! ! ===================================================================== ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC DBLE,DCONJG,MAX ! .. ! .. Local Scalars .. DOUBLE COMPLEX TEMP1,TEMP2 INTEGER I,INFO,J,L,NROWA LOGICAL UPPER ! .. ! .. Parameters .. DOUBLE PRECISION ONE PARAMETER (ONE=1.0D+0) DOUBLE COMPLEX ZERO PARAMETER (ZERO= (0.0D+0,0.0D+0)) ! .. ! ! Test the input parameters. ! IF (LSAME(TRANS,'N')) THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME(UPLO,'U') ! INFO = 0 IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN INFO = 1 ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND. & (.NOT.LSAME(TRANS,'C'))) THEN INFO = 2 ELSE IF (N.LT.0) THEN INFO = 3 ELSE IF (K.LT.0) THEN INFO = 4 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN INFO = 7 ELSE IF (LDB.LT.MAX(1,NROWA)) THEN INFO = 9 ELSE IF (LDC.LT.MAX(1,N)) THEN INFO = 12 END IF IF (INFO.NE.0) THEN CALL XERBLA('ZHER2K',INFO) RETURN END IF ! ! Quick return if possible. ! IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR. & (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN ! ! And when alpha.eq.zero. ! IF (ALPHA.EQ.ZERO) THEN IF (UPPER) THEN IF (BETA.EQ.DBLE(ZERO)) THEN DO 20 J = 1,N DO 10 I = 1,J C(I,J) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40 J = 1,N DO 30 I = 1,J - 1 C(I,J) = BETA*C(I,J) 30 CONTINUE C(J,J) = BETA*DBLE(C(J,J)) 40 CONTINUE END IF ELSE IF (BETA.EQ.DBLE(ZERO)) THEN DO 60 J = 1,N DO 50 I = J,N C(I,J) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80 J = 1,N C(J,J) = BETA*DBLE(C(J,J)) DO 70 I = J + 1,N C(I,J) = BETA*C(I,J) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF ! ! Start the operations. ! IF (LSAME(TRANS,'N')) THEN ! ! Form C := alpha*A*B**H + conjg( alpha )*B*A**H + ! C. ! IF (UPPER) THEN DO 130 J = 1,N IF (BETA.EQ.DBLE(ZERO)) THEN DO 90 I = 1,J C(I,J) = ZERO 90 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 100 I = 1,J - 1 C(I,J) = BETA*C(I,J) 100 CONTINUE C(J,J) = BETA*DBLE(C(J,J)) ELSE C(J,J) = DBLE(C(J,J)) END IF DO 120 L = 1,K IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN TEMP1 = ALPHA*DCONJG(B(J,L)) TEMP2 = DCONJG(ALPHA*A(J,L)) DO 110 I = 1,J - 1 C(I,J) = C(I,J) + A(I,L)*TEMP1 + & B(I,L)*TEMP2 110 CONTINUE C(J,J) = DBLE(C(J,J)) + & DBLE(A(J,L)*TEMP1+B(J,L)*TEMP2) END IF 120 CONTINUE 130 CONTINUE ELSE DO 180 J = 1,N IF (BETA.EQ.DBLE(ZERO)) THEN DO 140 I = J,N C(I,J) = ZERO 140 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 150 I = J + 1,N C(I,J) = BETA*C(I,J) 150 CONTINUE C(J,J) = BETA*DBLE(C(J,J)) ELSE C(J,J) = DBLE(C(J,J)) END IF DO 170 L = 1,K IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN TEMP1 = ALPHA*DCONJG(B(J,L)) TEMP2 = DCONJG(ALPHA*A(J,L)) DO 160 I = J + 1,N C(I,J) = C(I,J) + A(I,L)*TEMP1 + & B(I,L)*TEMP2 160 CONTINUE C(J,J) = DBLE(C(J,J)) + & DBLE(A(J,L)*TEMP1+B(J,L)*TEMP2) END IF 170 CONTINUE 180 CONTINUE END IF ELSE ! ! Form C := alpha*A**H*B + conjg( alpha )*B**H*A + ! C. ! IF (UPPER) THEN DO 210 J = 1,N DO 200 I = 1,J TEMP1 = ZERO TEMP2 = ZERO DO 190 L = 1,K TEMP1 = TEMP1 + DCONJG(A(L,I))*B(L,J) TEMP2 = TEMP2 + DCONJG(B(L,I))*A(L,J) 190 CONTINUE IF (I.EQ.J) THEN IF (BETA.EQ.DBLE(ZERO)) THEN C(J,J) = DBLE(ALPHA*TEMP1+ & DCONJG(ALPHA)*TEMP2) ELSE C(J,J) = BETA*DBLE(C(J,J)) + & DBLE(ALPHA*TEMP1+ & DCONJG(ALPHA)*TEMP2) END IF ELSE IF (BETA.EQ.DBLE(ZERO)) THEN C(I,J) = ALPHA*TEMP1 + DCONJG(ALPHA)*TEMP2 ELSE C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 + & DCONJG(ALPHA)*TEMP2 END IF END IF 200 CONTINUE 210 CONTINUE ELSE DO 240 J = 1,N DO 230 I = J,N TEMP1 = ZERO TEMP2 = ZERO DO 220 L = 1,K TEMP1 = TEMP1 + DCONJG(A(L,I))*B(L,J) TEMP2 = TEMP2 + DCONJG(B(L,I))*A(L,J) 220 CONTINUE IF (I.EQ.J) THEN IF (BETA.EQ.DBLE(ZERO)) THEN C(J,J) = DBLE(ALPHA*TEMP1+ & DCONJG(ALPHA)*TEMP2) ELSE C(J,J) = BETA*DBLE(C(J,J)) + & DBLE(ALPHA*TEMP1+ & DCONJG(ALPHA)*TEMP2) END IF ELSE IF (BETA.EQ.DBLE(ZERO)) THEN C(I,J) = ALPHA*TEMP1 + DCONJG(ALPHA)*TEMP2 ELSE C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 + & DCONJG(ALPHA)*TEMP2 END IF END IF 230 CONTINUE 240 CONTINUE END IF END IF ! RETURN ! ! End of ZHER2K. ! END SUBROUTINE ZHERK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC) !! ZHERK performs C := alpha*A*conjg( A' ) + beta*C, C hermitian. ! ! .. Scalar Arguments .. DOUBLE PRECISION ALPHA,BETA INTEGER K,LDA,LDC,N CHARACTER TRANS,UPLO ! .. ! .. Array Arguments .. DOUBLE COMPLEX A(LDA,*),C(LDC,*) ! .. ! ! Purpose ! ======= ! ! ZHERK performs one of the hermitian rank k operations ! ! C := alpha*A*A**H + beta*C, ! ! or ! ! C := alpha*A**H*A + beta*C, ! ! where alpha and beta are real scalars, C is an n by n hermitian ! matrix and A is an n by k matrix in the first case and a k by n ! matrix in the second case. ! ! Arguments ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array C is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of C ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of C ! is to be referenced. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' C := alpha*A*A**H + beta*C. ! ! TRANS = 'C' or 'c' C := alpha*A**H*A + beta*C. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry with TRANS = 'N' or 'n', K specifies the number ! of columns of the matrix A, and on entry with ! TRANS = 'C' or 'c', K specifies the number of rows of the ! matrix A. K must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is ! k when TRANS = 'N' or 'n', and is n otherwise. ! Before entry with TRANS = 'N' or 'n', the leading n by k ! part of the array A must contain the matrix A, otherwise ! the leading k by n part of the array A must contain the ! matrix A. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When TRANS = 'N' or 'n' ! then LDA must be at least max( 1, n ), otherwise LDA must ! be at least max( 1, k ). ! Unchanged on exit. ! ! BETA - DOUBLE PRECISION. ! On entry, BETA specifies the scalar beta. ! Unchanged on exit. ! ! C - COMPLEX*16 array of DIMENSION ( LDC, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array C must contain the upper ! triangular part of the hermitian matrix and the strictly ! lower triangular part of C is not referenced. On exit, the ! upper triangular part of the array C is overwritten by the ! upper triangular part of the updated matrix. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array C must contain the lower ! triangular part of the hermitian matrix and the strictly ! upper triangular part of C is not referenced. On exit, the ! lower triangular part of the array C is overwritten by the ! lower triangular part of the updated matrix. ! Note that the imaginary parts of the diagonal elements need ! not be set, they are assumed to be zero, and on exit they ! are set to zero. ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, n ). ! Unchanged on exit. ! ! Further Details ! =============== ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1. ! Ed Anderson, Cray Research Inc. ! ! ===================================================================== ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC DBLE,DCMPLX,DCONJG,MAX ! .. ! .. Local Scalars .. DOUBLE COMPLEX TEMP DOUBLE PRECISION RTEMP INTEGER I,INFO,J,L,NROWA LOGICAL UPPER ! .. ! .. Parameters .. DOUBLE PRECISION ONE,ZERO PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) ! .. ! ! Test the input parameters. ! IF (LSAME(TRANS,'N')) THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME(UPLO,'U') ! INFO = 0 IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN INFO = 1 ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND. & (.NOT.LSAME(TRANS,'C'))) THEN INFO = 2 ELSE IF (N.LT.0) THEN INFO = 3 ELSE IF (K.LT.0) THEN INFO = 4 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN INFO = 7 ELSE IF (LDC.LT.MAX(1,N)) THEN INFO = 10 END IF IF (INFO.NE.0) THEN CALL XERBLA('ZHERK ',INFO) RETURN END IF ! ! Quick return if possible. ! IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR. & (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN ! ! And when alpha.eq.zero. ! IF (ALPHA.EQ.ZERO) THEN IF (UPPER) THEN IF (BETA.EQ.ZERO) THEN DO 20 J = 1,N DO 10 I = 1,J C(I,J) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40 J = 1,N DO 30 I = 1,J - 1 C(I,J) = BETA*C(I,J) 30 CONTINUE C(J,J) = BETA*DBLE(C(J,J)) 40 CONTINUE END IF ELSE IF (BETA.EQ.ZERO) THEN DO 60 J = 1,N DO 50 I = J,N C(I,J) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80 J = 1,N C(J,J) = BETA*DBLE(C(J,J)) DO 70 I = J + 1,N C(I,J) = BETA*C(I,J) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF ! ! Start the operations. ! IF (LSAME(TRANS,'N')) THEN ! ! Form C := alpha*A*A**H + beta*C. ! IF (UPPER) THEN DO 130 J = 1,N IF (BETA.EQ.ZERO) THEN DO 90 I = 1,J C(I,J) = ZERO 90 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 100 I = 1,J - 1 C(I,J) = BETA*C(I,J) 100 CONTINUE C(J,J) = BETA*DBLE(C(J,J)) ELSE C(J,J) = DBLE(C(J,J)) END IF DO 120 L = 1,K IF (A(J,L).NE.DCMPLX(ZERO)) THEN TEMP = ALPHA*DCONJG(A(J,L)) DO 110 I = 1,J - 1 C(I,J) = C(I,J) + TEMP*A(I,L) 110 CONTINUE C(J,J) = DBLE(C(J,J)) + DBLE(TEMP*A(I,L)) END IF 120 CONTINUE 130 CONTINUE ELSE DO 180 J = 1,N IF (BETA.EQ.ZERO) THEN DO 140 I = J,N C(I,J) = ZERO 140 CONTINUE ELSE IF (BETA.NE.ONE) THEN C(J,J) = BETA*DBLE(C(J,J)) DO 150 I = J + 1,N C(I,J) = BETA*C(I,J) 150 CONTINUE ELSE C(J,J) = DBLE(C(J,J)) END IF DO 170 L = 1,K IF (A(J,L).NE.DCMPLX(ZERO)) THEN TEMP = ALPHA*DCONJG(A(J,L)) C(J,J) = DBLE(C(J,J)) + DBLE(TEMP*A(J,L)) DO 160 I = J + 1,N C(I,J) = C(I,J) + TEMP*A(I,L) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE END IF ELSE ! ! Form C := alpha*A**H*A + beta*C. ! IF (UPPER) THEN DO 220 J = 1,N DO 200 I = 1,J - 1 TEMP = ZERO DO 190 L = 1,K TEMP = TEMP + DCONJG(A(L,I))*A(L,J) 190 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP ELSE C(I,J) = ALPHA*TEMP + BETA*C(I,J) END IF 200 CONTINUE RTEMP = ZERO DO 210 L = 1,K RTEMP = RTEMP + DCONJG(A(L,J))*A(L,J) 210 CONTINUE IF (BETA.EQ.ZERO) THEN C(J,J) = ALPHA*RTEMP ELSE C(J,J) = ALPHA*RTEMP + BETA*DBLE(C(J,J)) END IF 220 CONTINUE ELSE DO 260 J = 1,N RTEMP = ZERO DO 230 L = 1,K RTEMP = RTEMP + DCONJG(A(L,J))*A(L,J) 230 CONTINUE IF (BETA.EQ.ZERO) THEN C(J,J) = ALPHA*RTEMP ELSE C(J,J) = ALPHA*RTEMP + BETA*DBLE(C(J,J)) END IF DO 250 I = J + 1,N TEMP = ZERO DO 240 L = 1,K TEMP = TEMP + DCONJG(A(L,I))*A(L,J) 240 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP ELSE C(I,J) = ALPHA*TEMP + BETA*C(I,J) END IF 250 CONTINUE 260 CONTINUE END IF END IF ! RETURN ! ! End of ZHERK . ! END SUBROUTINE ZSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC) !! ZSYMM performs C:=alpha*A*B+beta*C, A symmetric, B and C rectangular. ! ! .. Scalar Arguments .. DOUBLE COMPLEX ALPHA,BETA INTEGER LDA,LDB,LDC,M,N CHARACTER SIDE,UPLO ! .. ! .. Array Arguments .. DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) ! .. ! ! Purpose ! ======= ! ! ZSYMM performs one of the matrix-matrix operations ! ! C := alpha*A*B + beta*C, ! ! or ! ! C := alpha*B*A + beta*C, ! ! where alpha and beta are scalars, A is a symmetric matrix and B and ! C are m by n matrices. ! ! Arguments ! ========== ! ! SIDE - CHARACTER*1. ! On entry, SIDE specifies whether the symmetric matrix A ! appears on the left or right in the operation as follows: ! ! SIDE = 'L' or 'l' C := alpha*A*B + beta*C, ! ! SIDE = 'R' or 'r' C := alpha*B*A + beta*C, ! ! Unchanged on exit. ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the symmetric matrix A is to be ! referenced as follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of the ! symmetric matrix is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of the ! symmetric matrix is to be referenced. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix C. ! M must be at least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix C. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - COMPLEX*16 . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is ! m when SIDE = 'L' or 'l' and is n otherwise. ! Before entry with SIDE = 'L' or 'l', the m by m part of ! the array A must contain the symmetric matrix, such that ! when UPLO = 'U' or 'u', the leading m by m upper triangular ! part of the array A must contain the upper triangular part ! of the symmetric matrix and the strictly lower triangular ! part of A is not referenced, and when UPLO = 'L' or 'l', ! the leading m by m lower triangular part of the array A ! must contain the lower triangular part of the symmetric ! matrix and the strictly upper triangular part of A is not ! referenced. ! Before entry with SIDE = 'R' or 'r', the n by n part of ! the array A must contain the symmetric matrix, such that ! when UPLO = 'U' or 'u', the leading n by n upper triangular ! part of the array A must contain the upper triangular part ! of the symmetric matrix and the strictly lower triangular ! part of A is not referenced, and when UPLO = 'L' or 'l', ! the leading n by n lower triangular part of the array A ! must contain the lower triangular part of the symmetric ! matrix and the strictly upper triangular part of A is not ! referenced. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When SIDE = 'L' or 'l' then ! LDA must be at least max( 1, m ), otherwise LDA must be at ! least max( 1, n ). ! Unchanged on exit. ! ! B - COMPLEX*16 array of DIMENSION ( LDB, n ). ! Before entry, the leading m by n part of the array B must ! contain the matrix B. ! Unchanged on exit. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. LDB must be at least ! max( 1, m ). ! Unchanged on exit. ! ! BETA - COMPLEX*16 . ! On entry, BETA specifies the scalar beta. When BETA is ! supplied as zero then C need not be set on input. ! Unchanged on exit. ! ! C - COMPLEX*16 array of DIMENSION ( LDC, n ). ! Before entry, the leading m by n part of the array C must ! contain the matrix C, except when beta is zero, in which ! case C need not be set on entry. ! On exit, the array C is overwritten by the m by n updated ! matrix. ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, m ). ! Unchanged on exit. ! ! Further Details ! =============== ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ===================================================================== ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Local Scalars .. DOUBLE COMPLEX TEMP1,TEMP2 INTEGER I,INFO,J,K,NROWA LOGICAL UPPER ! .. ! .. Parameters .. DOUBLE COMPLEX ONE PARAMETER (ONE= (1.0D+0,0.0D+0)) DOUBLE COMPLEX ZERO PARAMETER (ZERO= (0.0D+0,0.0D+0)) ! .. ! ! Set NROWA as the number of rows of A. ! IF (LSAME(SIDE,'L')) THEN NROWA = M ELSE NROWA = N END IF UPPER = LSAME(UPLO,'U') ! ! Test the input parameters. ! INFO = 0 IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN INFO = 1 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN INFO = 2 ELSE IF (M.LT.0) THEN INFO = 3 ELSE IF (N.LT.0) THEN INFO = 4 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN INFO = 7 ELSE IF (LDB.LT.MAX(1,M)) THEN INFO = 9 ELSE IF (LDC.LT.MAX(1,M)) THEN INFO = 12 END IF IF (INFO.NE.0) THEN CALL XERBLA('ZSYMM ',INFO) RETURN END IF ! ! Quick return if possible. ! IF ((M.EQ.0) .OR. (N.EQ.0) .OR. & ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN ! ! And when alpha.eq.zero. ! IF (ALPHA.EQ.ZERO) THEN IF (BETA.EQ.ZERO) THEN DO 20 J = 1,N DO 10 I = 1,M C(I,J) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40 J = 1,N DO 30 I = 1,M C(I,J) = BETA*C(I,J) 30 CONTINUE 40 CONTINUE END IF RETURN END IF ! ! Start the operations. ! IF (LSAME(SIDE,'L')) THEN ! ! Form C := alpha*A*B + beta*C. ! IF (UPPER) THEN DO 70 J = 1,N DO 60 I = 1,M TEMP1 = ALPHA*B(I,J) TEMP2 = ZERO DO 50 K = 1,I - 1 C(K,J) = C(K,J) + TEMP1*A(K,I) TEMP2 = TEMP2 + B(K,J)*A(K,I) 50 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2 ELSE C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) + & ALPHA*TEMP2 END IF 60 CONTINUE 70 CONTINUE ELSE DO 100 J = 1,N DO 90 I = M,1,-1 TEMP1 = ALPHA*B(I,J) TEMP2 = ZERO DO 80 K = I + 1,M C(K,J) = C(K,J) + TEMP1*A(K,I) TEMP2 = TEMP2 + B(K,J)*A(K,I) 80 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2 ELSE C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) + & ALPHA*TEMP2 END IF 90 CONTINUE 100 CONTINUE END IF ELSE ! ! Form C := alpha*B*A + beta*C. ! DO 170 J = 1,N TEMP1 = ALPHA*A(J,J) IF (BETA.EQ.ZERO) THEN DO 110 I = 1,M C(I,J) = TEMP1*B(I,J) 110 CONTINUE ELSE DO 120 I = 1,M C(I,J) = BETA*C(I,J) + TEMP1*B(I,J) 120 CONTINUE END IF DO 140 K = 1,J - 1 IF (UPPER) THEN TEMP1 = ALPHA*A(K,J) ELSE TEMP1 = ALPHA*A(J,K) END IF DO 130 I = 1,M C(I,J) = C(I,J) + TEMP1*B(I,K) 130 CONTINUE 140 CONTINUE DO 160 K = J + 1,N IF (UPPER) THEN TEMP1 = ALPHA*A(J,K) ELSE TEMP1 = ALPHA*A(K,J) END IF DO 150 I = 1,M C(I,J) = C(I,J) + TEMP1*B(I,K) 150 CONTINUE 160 CONTINUE 170 CONTINUE END IF ! RETURN ! ! End of ZSYMM . ! END SUBROUTINE ZSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) !! ZSYR2K performs C := alpha*A*B' + alpha*B*A' + beta*C, C symmetric. ! ! .. Scalar Arguments .. DOUBLE COMPLEX ALPHA,BETA INTEGER K,LDA,LDB,LDC,N CHARACTER TRANS,UPLO ! .. ! .. Array Arguments .. DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) ! .. ! ! Purpose ! ======= ! ! ZSYR2K performs one of the symmetric rank 2k operations ! ! C := alpha*A*B**T + alpha*B*A**T + beta*C, ! ! or ! ! C := alpha*A**T*B + alpha*B**T*A + beta*C, ! ! where alpha and beta are scalars, C is an n by n symmetric matrix ! and A and B are n by k matrices in the first case and k by n ! matrices in the second case. ! ! Arguments ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array C is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of C ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of C ! is to be referenced. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' C := alpha*A*B**T + alpha*B*A**T + ! beta*C. ! ! TRANS = 'T' or 't' C := alpha*A**T*B + alpha*B**T*A + ! beta*C. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry with TRANS = 'N' or 'n', K specifies the number ! of columns of the matrices A and B, and on entry with ! TRANS = 'T' or 't', K specifies the number of rows of the ! matrices A and B. K must be at least zero. ! Unchanged on exit. ! ! ALPHA - COMPLEX*16 . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is ! k when TRANS = 'N' or 'n', and is n otherwise. ! Before entry with TRANS = 'N' or 'n', the leading n by k ! part of the array A must contain the matrix A, otherwise ! the leading k by n part of the array A must contain the ! matrix A. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When TRANS = 'N' or 'n' ! then LDA must be at least max( 1, n ), otherwise LDA must ! be at least max( 1, k ). ! Unchanged on exit. ! ! B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is ! k when TRANS = 'N' or 'n', and is n otherwise. ! Before entry with TRANS = 'N' or 'n', the leading n by k ! part of the array B must contain the matrix B, otherwise ! the leading k by n part of the array B must contain the ! matrix B. ! Unchanged on exit. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. When TRANS = 'N' or 'n' ! then LDB must be at least max( 1, n ), otherwise LDB must ! be at least max( 1, k ). ! Unchanged on exit. ! ! BETA - COMPLEX*16 . ! On entry, BETA specifies the scalar beta. ! Unchanged on exit. ! ! C - COMPLEX*16 array of DIMENSION ( LDC, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array C must contain the upper ! triangular part of the symmetric matrix and the strictly ! lower triangular part of C is not referenced. On exit, the ! upper triangular part of the array C is overwritten by the ! upper triangular part of the updated matrix. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array C must contain the lower ! triangular part of the symmetric matrix and the strictly ! upper triangular part of C is not referenced. On exit, the ! lower triangular part of the array C is overwritten by the ! lower triangular part of the updated matrix. ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, n ). ! Unchanged on exit. ! ! Further Details ! =============== ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ===================================================================== ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Local Scalars .. DOUBLE COMPLEX TEMP1,TEMP2 INTEGER I,INFO,J,L,NROWA LOGICAL UPPER ! .. ! .. Parameters .. DOUBLE COMPLEX ONE PARAMETER (ONE= (1.0D+0,0.0D+0)) DOUBLE COMPLEX ZERO PARAMETER (ZERO= (0.0D+0,0.0D+0)) ! .. ! ! Test the input parameters. ! IF (LSAME(TRANS,'N')) THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME(UPLO,'U') ! INFO = 0 IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN INFO = 1 ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND. & (.NOT.LSAME(TRANS,'T'))) THEN INFO = 2 ELSE IF (N.LT.0) THEN INFO = 3 ELSE IF (K.LT.0) THEN INFO = 4 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN INFO = 7 ELSE IF (LDB.LT.MAX(1,NROWA)) THEN INFO = 9 ELSE IF (LDC.LT.MAX(1,N)) THEN INFO = 12 END IF IF (INFO.NE.0) THEN CALL XERBLA('ZSYR2K',INFO) RETURN END IF ! ! Quick return if possible. ! IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR. & (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN ! ! And when alpha.eq.zero. ! IF (ALPHA.EQ.ZERO) THEN IF (UPPER) THEN IF (BETA.EQ.ZERO) THEN DO 20 J = 1,N DO 10 I = 1,J C(I,J) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40 J = 1,N DO 30 I = 1,J C(I,J) = BETA*C(I,J) 30 CONTINUE 40 CONTINUE END IF ELSE IF (BETA.EQ.ZERO) THEN DO 60 J = 1,N DO 50 I = J,N C(I,J) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80 J = 1,N DO 70 I = J,N C(I,J) = BETA*C(I,J) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF ! ! Start the operations. ! IF (LSAME(TRANS,'N')) THEN ! ! Form C := alpha*A*B**T + alpha*B*A**T + C. ! IF (UPPER) THEN DO 130 J = 1,N IF (BETA.EQ.ZERO) THEN DO 90 I = 1,J C(I,J) = ZERO 90 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 100 I = 1,J C(I,J) = BETA*C(I,J) 100 CONTINUE END IF DO 120 L = 1,K IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN TEMP1 = ALPHA*B(J,L) TEMP2 = ALPHA*A(J,L) DO 110 I = 1,J C(I,J) = C(I,J) + A(I,L)*TEMP1 + & B(I,L)*TEMP2 110 CONTINUE END IF 120 CONTINUE 130 CONTINUE ELSE DO 180 J = 1,N IF (BETA.EQ.ZERO) THEN DO 140 I = J,N C(I,J) = ZERO 140 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 150 I = J,N C(I,J) = BETA*C(I,J) 150 CONTINUE END IF DO 170 L = 1,K IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN TEMP1 = ALPHA*B(J,L) TEMP2 = ALPHA*A(J,L) DO 160 I = J,N C(I,J) = C(I,J) + A(I,L)*TEMP1 + & B(I,L)*TEMP2 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE END IF ELSE ! ! Form C := alpha*A**T*B + alpha*B**T*A + C. ! IF (UPPER) THEN DO 210 J = 1,N DO 200 I = 1,J TEMP1 = ZERO TEMP2 = ZERO DO 190 L = 1,K TEMP1 = TEMP1 + A(L,I)*B(L,J) TEMP2 = TEMP2 + B(L,I)*A(L,J) 190 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2 ELSE C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 + & ALPHA*TEMP2 END IF 200 CONTINUE 210 CONTINUE ELSE DO 240 J = 1,N DO 230 I = J,N TEMP1 = ZERO TEMP2 = ZERO DO 220 L = 1,K TEMP1 = TEMP1 + A(L,I)*B(L,J) TEMP2 = TEMP2 + B(L,I)*A(L,J) 220 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2 ELSE C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 + & ALPHA*TEMP2 END IF 230 CONTINUE 240 CONTINUE END IF END IF ! RETURN ! ! End of ZSYR2K. ! END SUBROUTINE ZSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC) !! ZSYRK performs C := alpha*A*A' + beta*C, C is symmetric. ! ! .. Scalar Arguments .. DOUBLE COMPLEX ALPHA,BETA INTEGER K,LDA,LDC,N CHARACTER TRANS,UPLO ! .. ! .. Array Arguments .. DOUBLE COMPLEX A(LDA,*),C(LDC,*) ! .. ! ! Purpose ! ======= ! ! ZSYRK performs one of the symmetric rank k operations ! ! C := alpha*A*A**T + beta*C, ! ! or ! ! C := alpha*A**T*A + beta*C, ! ! where alpha and beta are scalars, C is an n by n symmetric matrix ! and A is an n by k matrix in the first case and a k by n matrix ! in the second case. ! ! Arguments ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array C is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of C ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of C ! is to be referenced. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' C := alpha*A*A**T + beta*C. ! ! TRANS = 'T' or 't' C := alpha*A**T*A + beta*C. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry with TRANS = 'N' or 'n', K specifies the number ! of columns of the matrix A, and on entry with ! TRANS = 'T' or 't', K specifies the number of rows of the ! matrix A. K must be at least zero. ! Unchanged on exit. ! ! ALPHA - COMPLEX*16 . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is ! k when TRANS = 'N' or 'n', and is n otherwise. ! Before entry with TRANS = 'N' or 'n', the leading n by k ! part of the array A must contain the matrix A, otherwise ! the leading k by n part of the array A must contain the ! matrix A. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When TRANS = 'N' or 'n' ! then LDA must be at least max( 1, n ), otherwise LDA must ! be at least max( 1, k ). ! Unchanged on exit. ! ! BETA - COMPLEX*16 . ! On entry, BETA specifies the scalar beta. ! Unchanged on exit. ! ! C - COMPLEX*16 array of DIMENSION ( LDC, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array C must contain the upper ! triangular part of the symmetric matrix and the strictly ! lower triangular part of C is not referenced. On exit, the ! upper triangular part of the array C is overwritten by the ! upper triangular part of the updated matrix. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array C must contain the lower ! triangular part of the symmetric matrix and the strictly ! upper triangular part of C is not referenced. On exit, the ! lower triangular part of the array C is overwritten by the ! lower triangular part of the updated matrix. ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, n ). ! Unchanged on exit. ! ! Further Details ! =============== ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ===================================================================== ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Local Scalars .. DOUBLE COMPLEX TEMP INTEGER I,INFO,J,L,NROWA LOGICAL UPPER ! .. ! .. Parameters .. DOUBLE COMPLEX ONE PARAMETER (ONE= (1.0D+0,0.0D+0)) DOUBLE COMPLEX ZERO PARAMETER (ZERO= (0.0D+0,0.0D+0)) ! .. ! ! Test the input parameters. ! IF (LSAME(TRANS,'N')) THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME(UPLO,'U') ! INFO = 0 IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN INFO = 1 ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND. & (.NOT.LSAME(TRANS,'T'))) THEN INFO = 2 ELSE IF (N.LT.0) THEN INFO = 3 ELSE IF (K.LT.0) THEN INFO = 4 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN INFO = 7 ELSE IF (LDC.LT.MAX(1,N)) THEN INFO = 10 END IF IF (INFO.NE.0) THEN CALL XERBLA('ZSYRK ',INFO) RETURN END IF ! ! Quick return if possible. ! IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR. & (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN ! ! And when alpha.eq.zero. ! IF (ALPHA.EQ.ZERO) THEN IF (UPPER) THEN IF (BETA.EQ.ZERO) THEN DO 20 J = 1,N DO 10 I = 1,J C(I,J) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40 J = 1,N DO 30 I = 1,J C(I,J) = BETA*C(I,J) 30 CONTINUE 40 CONTINUE END IF ELSE IF (BETA.EQ.ZERO) THEN DO 60 J = 1,N DO 50 I = J,N C(I,J) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80 J = 1,N DO 70 I = J,N C(I,J) = BETA*C(I,J) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF ! ! Start the operations. ! IF (LSAME(TRANS,'N')) THEN ! ! Form C := alpha*A*A**T + beta*C. ! IF (UPPER) THEN DO 130 J = 1,N IF (BETA.EQ.ZERO) THEN DO 90 I = 1,J C(I,J) = ZERO 90 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 100 I = 1,J C(I,J) = BETA*C(I,J) 100 CONTINUE END IF DO 120 L = 1,K IF (A(J,L).NE.ZERO) THEN TEMP = ALPHA*A(J,L) DO 110 I = 1,J C(I,J) = C(I,J) + TEMP*A(I,L) 110 CONTINUE END IF 120 CONTINUE 130 CONTINUE ELSE DO 180 J = 1,N IF (BETA.EQ.ZERO) THEN DO 140 I = J,N C(I,J) = ZERO 140 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 150 I = J,N C(I,J) = BETA*C(I,J) 150 CONTINUE END IF DO 170 L = 1,K IF (A(J,L).NE.ZERO) THEN TEMP = ALPHA*A(J,L) DO 160 I = J,N C(I,J) = C(I,J) + TEMP*A(I,L) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE END IF ELSE ! ! Form C := alpha*A**T*A + beta*C. ! IF (UPPER) THEN DO 210 J = 1,N DO 200 I = 1,J TEMP = ZERO DO 190 L = 1,K TEMP = TEMP + A(L,I)*A(L,J) 190 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP ELSE C(I,J) = ALPHA*TEMP + BETA*C(I,J) END IF 200 CONTINUE 210 CONTINUE ELSE DO 240 J = 1,N DO 230 I = J,N TEMP = ZERO DO 220 L = 1,K TEMP = TEMP + A(L,I)*A(L,J) 220 CONTINUE IF (BETA.EQ.ZERO) THEN C(I,J) = ALPHA*TEMP ELSE C(I,J) = ALPHA*TEMP + BETA*C(I,J) END IF 230 CONTINUE 240 CONTINUE END IF END IF ! RETURN ! ! End of ZSYRK . ! END SUBROUTINE ZTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) !! ZTRMM performs B := alpha*op( A ) * B. where A is triangular. ! ! .. Scalar Arguments .. DOUBLE COMPLEX ALPHA INTEGER LDA,LDB,M,N CHARACTER DIAG,SIDE,TRANSA,UPLO ! .. ! .. Array Arguments .. DOUBLE COMPLEX A(LDA,*),B(LDB,*) ! .. ! ! Purpose ! ======= ! ! ZTRMM performs one of the matrix-matrix operations ! ! B := alpha*op( A )*B, or B := alpha*B*op( A ) ! ! where alpha is a scalar, B is an m by n matrix, A is a unit, or ! non-unit, upper or lower triangular matrix and op( A ) is one of ! ! op( A ) = A or op( A ) = A**T or op( A ) = A**H. ! ! Arguments ! ========== ! ! SIDE - CHARACTER*1. ! On entry, SIDE specifies whether op( A ) multiplies B from ! the left or right as follows: ! ! SIDE = 'L' or 'l' B := alpha*op( A )*B. ! ! SIDE = 'R' or 'r' B := alpha*B*op( A ). ! ! Unchanged on exit. ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix A is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANSA - CHARACTER*1. ! On entry, TRANSA specifies the form of op( A ) to be used in ! the matrix multiplication as follows: ! ! TRANSA = 'N' or 'n' op( A ) = A. ! ! TRANSA = 'T' or 't' op( A ) = A**T. ! ! TRANSA = 'C' or 'c' op( A ) = A**H. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit triangular ! as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of B. M must be at ! least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of B. N must be ! at least zero. ! Unchanged on exit. ! ! ALPHA - COMPLEX*16 . ! On entry, ALPHA specifies the scalar alpha. When alpha is ! zero then A is not referenced and B need not be set before ! entry. ! Unchanged on exit. ! ! A - COMPLEX*16 array of DIMENSION ( LDA, k ), where k is m ! when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. ! Before entry with UPLO = 'U' or 'u', the leading k by k ! upper triangular part of the array A must contain the upper ! triangular matrix and the strictly lower triangular part of ! A is not referenced. ! Before entry with UPLO = 'L' or 'l', the leading k by k ! lower triangular part of the array A must contain the lower ! triangular matrix and the strictly upper triangular part of ! A is not referenced. ! Note that when DIAG = 'U' or 'u', the diagonal elements of ! A are not referenced either, but are assumed to be unity. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When SIDE = 'L' or 'l' then ! LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' ! then LDA must be at least max( 1, n ). ! Unchanged on exit. ! ! B - COMPLEX*16 array of DIMENSION ( LDB, n ). ! Before entry, the leading m by n part of the array B must ! contain the matrix B, and on exit is overwritten by the ! transformed matrix. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. LDB must be at least ! max( 1, m ). ! Unchanged on exit. ! ! Further Details ! =============== ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ===================================================================== ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC DCONJG,MAX ! .. ! .. Local Scalars .. DOUBLE COMPLEX TEMP INTEGER I,INFO,J,K,NROWA LOGICAL LSIDE,NOCONJ,NOUNIT,UPPER ! .. ! .. Parameters .. DOUBLE COMPLEX ONE PARAMETER (ONE= (1.0D+0,0.0D+0)) DOUBLE COMPLEX ZERO PARAMETER (ZERO= (0.0D+0,0.0D+0)) ! .. ! ! Test the input parameters. ! LSIDE = LSAME(SIDE,'L') IF (LSIDE) THEN NROWA = M ELSE NROWA = N END IF NOCONJ = LSAME(TRANSA,'T') NOUNIT = LSAME(DIAG,'N') UPPER = LSAME(UPLO,'U') ! INFO = 0 IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN INFO = 1 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN INFO = 2 ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND. & (.NOT.LSAME(TRANSA,'T')) .AND. & (.NOT.LSAME(TRANSA,'C'))) THEN INFO = 3 ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN INFO = 4 ELSE IF (M.LT.0) THEN INFO = 5 ELSE IF (N.LT.0) THEN INFO = 6 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN INFO = 9 ELSE IF (LDB.LT.MAX(1,M)) THEN INFO = 11 END IF IF (INFO.NE.0) THEN CALL XERBLA('ZTRMM ',INFO) RETURN END IF ! ! Quick return if possible. ! IF (M.EQ.0 .OR. N.EQ.0) RETURN ! ! And when alpha.eq.zero. ! IF (ALPHA.EQ.ZERO) THEN DO 20 J = 1,N DO 10 I = 1,M B(I,J) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF ! ! Start the operations. ! IF (LSIDE) THEN IF (LSAME(TRANSA,'N')) THEN ! ! Form B := alpha*A*B. ! IF (UPPER) THEN DO 50 J = 1,N DO 40 K = 1,M IF (B(K,J).NE.ZERO) THEN TEMP = ALPHA*B(K,J) DO 30 I = 1,K - 1 B(I,J) = B(I,J) + TEMP*A(I,K) 30 CONTINUE IF (NOUNIT) TEMP = TEMP*A(K,K) B(K,J) = TEMP END IF 40 CONTINUE 50 CONTINUE ELSE DO 80 J = 1,N DO 70 K = M,1,-1 IF (B(K,J).NE.ZERO) THEN TEMP = ALPHA*B(K,J) B(K,J) = TEMP IF (NOUNIT) B(K,J) = B(K,J)*A(K,K) DO 60 I = K + 1,M B(I,J) = B(I,J) + TEMP*A(I,K) 60 CONTINUE END IF 70 CONTINUE 80 CONTINUE END IF ELSE ! ! Form B := alpha*A**T*B or B := alpha*A**H*B. ! IF (UPPER) THEN DO 120 J = 1,N DO 110 I = M,1,-1 TEMP = B(I,J) IF (NOCONJ) THEN IF (NOUNIT) TEMP = TEMP*A(I,I) DO 90 K = 1,I - 1 TEMP = TEMP + A(K,I)*B(K,J) 90 CONTINUE ELSE IF (NOUNIT) TEMP = TEMP*DCONJG(A(I,I)) DO 100 K = 1,I - 1 TEMP = TEMP + DCONJG(A(K,I))*B(K,J) 100 CONTINUE END IF B(I,J) = ALPHA*TEMP 110 CONTINUE 120 CONTINUE ELSE DO 160 J = 1,N DO 150 I = 1,M TEMP = B(I,J) IF (NOCONJ) THEN IF (NOUNIT) TEMP = TEMP*A(I,I) DO 130 K = I + 1,M TEMP = TEMP + A(K,I)*B(K,J) 130 CONTINUE ELSE IF (NOUNIT) TEMP = TEMP*DCONJG(A(I,I)) DO 140 K = I + 1,M TEMP = TEMP + DCONJG(A(K,I))*B(K,J) 140 CONTINUE END IF B(I,J) = ALPHA*TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF (LSAME(TRANSA,'N')) THEN ! ! Form B := alpha*B*A. ! IF (UPPER) THEN DO 200 J = N,1,-1 TEMP = ALPHA IF (NOUNIT) TEMP = TEMP*A(J,J) DO 170 I = 1,M B(I,J) = TEMP*B(I,J) 170 CONTINUE DO 190 K = 1,J - 1 IF (A(K,J).NE.ZERO) THEN TEMP = ALPHA*A(K,J) DO 180 I = 1,M B(I,J) = B(I,J) + TEMP*B(I,K) 180 CONTINUE END IF 190 CONTINUE 200 CONTINUE ELSE DO 240 J = 1,N TEMP = ALPHA IF (NOUNIT) TEMP = TEMP*A(J,J) DO 210 I = 1,M B(I,J) = TEMP*B(I,J) 210 CONTINUE DO 230 K = J + 1,N IF (A(K,J).NE.ZERO) THEN TEMP = ALPHA*A(K,J) DO 220 I = 1,M B(I,J) = B(I,J) + TEMP*B(I,K) 220 CONTINUE END IF 230 CONTINUE 240 CONTINUE END IF ELSE ! ! Form B := alpha*B*A**T or B := alpha*B*A**H. ! IF (UPPER) THEN DO 280 K = 1,N DO 260 J = 1,K - 1 IF (A(J,K).NE.ZERO) THEN IF (NOCONJ) THEN TEMP = ALPHA*A(J,K) ELSE TEMP = ALPHA*DCONJG(A(J,K)) END IF DO 250 I = 1,M B(I,J) = B(I,J) + TEMP*B(I,K) 250 CONTINUE END IF 260 CONTINUE TEMP = ALPHA IF (NOUNIT) THEN IF (NOCONJ) THEN TEMP = TEMP*A(K,K) ELSE TEMP = TEMP*DCONJG(A(K,K)) END IF END IF IF (TEMP.NE.ONE) THEN DO 270 I = 1,M B(I,K) = TEMP*B(I,K) 270 CONTINUE END IF 280 CONTINUE ELSE DO 320 K = N,1,-1 DO 300 J = K + 1,N IF (A(J,K).NE.ZERO) THEN IF (NOCONJ) THEN TEMP = ALPHA*A(J,K) ELSE TEMP = ALPHA*DCONJG(A(J,K)) END IF DO 290 I = 1,M B(I,J) = B(I,J) + TEMP*B(I,K) 290 CONTINUE END IF 300 CONTINUE TEMP = ALPHA IF (NOUNIT) THEN IF (NOCONJ) THEN TEMP = TEMP*A(K,K) ELSE TEMP = TEMP*DCONJG(A(K,K)) END IF END IF IF (TEMP.NE.ONE) THEN DO 310 I = 1,M B(I,K) = TEMP*B(I,K) 310 CONTINUE END IF 320 CONTINUE END IF END IF END IF ! RETURN ! ! End of ZTRMM . ! END SUBROUTINE ZTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) !! ZTRSM solves op( A ) * x = alpha * B, where A is triangular. ! ! .. Scalar Arguments .. DOUBLE COMPLEX ALPHA INTEGER LDA,LDB,M,N CHARACTER DIAG,SIDE,TRANSA,UPLO ! .. ! .. Array Arguments .. DOUBLE COMPLEX A(LDA,*),B(LDB,*) ! .. ! ! Purpose ! ======= ! ! ZTRSM solves one of the matrix equations ! ! op( A )*X = alpha*B, or X*op( A ) = alpha*B, ! ! where alpha is a scalar, X and B are m by n matrices, A is a unit, or ! non-unit, upper or lower triangular matrix and op( A ) is one of ! ! op( A ) = A or op( A ) = A**T or op( A ) = A**H. ! ! The matrix X is overwritten on B. ! ! Arguments ! ========== ! ! SIDE - CHARACTER*1. ! On entry, SIDE specifies whether op( A ) appears on the left ! or right of X as follows: ! ! SIDE = 'L' or 'l' op( A )*X = alpha*B. ! ! SIDE = 'R' or 'r' X*op( A ) = alpha*B. ! ! Unchanged on exit. ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix A is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANSA - CHARACTER*1. ! On entry, TRANSA specifies the form of op( A ) to be used in ! the matrix multiplication as follows: ! ! TRANSA = 'N' or 'n' op( A ) = A. ! ! TRANSA = 'T' or 't' op( A ) = A**T. ! ! TRANSA = 'C' or 'c' op( A ) = A**H. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit triangular ! as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of B. M must be at ! least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of B. N must be ! at least zero. ! Unchanged on exit. ! ! ALPHA - COMPLEX*16 . ! On entry, ALPHA specifies the scalar alpha. When alpha is ! zero then A is not referenced and B need not be set before ! entry. ! Unchanged on exit. ! ! A - COMPLEX*16 array of DIMENSION ( LDA, k ), where k is m ! when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. ! Before entry with UPLO = 'U' or 'u', the leading k by k ! upper triangular part of the array A must contain the upper ! triangular matrix and the strictly lower triangular part of ! A is not referenced. ! Before entry with UPLO = 'L' or 'l', the leading k by k ! lower triangular part of the array A must contain the lower ! triangular matrix and the strictly upper triangular part of ! A is not referenced. ! Note that when DIAG = 'U' or 'u', the diagonal elements of ! A are not referenced either, but are assumed to be unity. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When SIDE = 'L' or 'l' then ! LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' ! then LDA must be at least max( 1, n ). ! Unchanged on exit. ! ! B - COMPLEX*16 array of DIMENSION ( LDB, n ). ! Before entry, the leading m by n part of the array B must ! contain the right-hand side matrix B, and on exit is ! overwritten by the solution matrix X. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. LDB must be at least ! max( 1, m ). ! Unchanged on exit. ! ! Further Details ! =============== ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ===================================================================== ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC DCONJG,MAX ! .. ! .. Local Scalars .. DOUBLE COMPLEX TEMP INTEGER I,INFO,J,K,NROWA LOGICAL LSIDE,NOCONJ,NOUNIT,UPPER ! .. ! .. Parameters .. DOUBLE COMPLEX ONE PARAMETER (ONE= (1.0D+0,0.0D+0)) DOUBLE COMPLEX ZERO PARAMETER (ZERO= (0.0D+0,0.0D+0)) ! .. ! ! Test the input parameters. ! LSIDE = LSAME(SIDE,'L') IF (LSIDE) THEN NROWA = M ELSE NROWA = N END IF NOCONJ = LSAME(TRANSA,'T') NOUNIT = LSAME(DIAG,'N') UPPER = LSAME(UPLO,'U') ! INFO = 0 IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN INFO = 1 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN INFO = 2 ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND. & (.NOT.LSAME(TRANSA,'T')) .AND. & (.NOT.LSAME(TRANSA,'C'))) THEN INFO = 3 ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN INFO = 4 ELSE IF (M.LT.0) THEN INFO = 5 ELSE IF (N.LT.0) THEN INFO = 6 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN INFO = 9 ELSE IF (LDB.LT.MAX(1,M)) THEN INFO = 11 END IF IF (INFO.NE.0) THEN CALL XERBLA('ZTRSM ',INFO) RETURN END IF ! ! Quick return if possible. ! IF (M.EQ.0 .OR. N.EQ.0) RETURN ! ! And when alpha.eq.zero. ! IF (ALPHA.EQ.ZERO) THEN DO 20 J = 1,N DO 10 I = 1,M B(I,J) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF ! ! Start the operations. ! IF (LSIDE) THEN IF (LSAME(TRANSA,'N')) THEN ! ! Form B := alpha*inv( A )*B. ! IF (UPPER) THEN DO 60 J = 1,N IF (ALPHA.NE.ONE) THEN DO 30 I = 1,M B(I,J) = ALPHA*B(I,J) 30 CONTINUE END IF DO 50 K = M,1,-1 IF (B(K,J).NE.ZERO) THEN IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) DO 40 I = 1,K - 1 B(I,J) = B(I,J) - B(K,J)*A(I,K) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100 J = 1,N IF (ALPHA.NE.ONE) THEN DO 70 I = 1,M B(I,J) = ALPHA*B(I,J) 70 CONTINUE END IF DO 90 K = 1,M IF (B(K,J).NE.ZERO) THEN IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) DO 80 I = K + 1,M B(I,J) = B(I,J) - B(K,J)*A(I,K) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE ! ! Form B := alpha*inv( A**T )*B ! or B := alpha*inv( A**H )*B. ! IF (UPPER) THEN DO 140 J = 1,N DO 130 I = 1,M TEMP = ALPHA*B(I,J) IF (NOCONJ) THEN DO 110 K = 1,I - 1 TEMP = TEMP - A(K,I)*B(K,J) 110 CONTINUE IF (NOUNIT) TEMP = TEMP/A(I,I) ELSE DO 120 K = 1,I - 1 TEMP = TEMP - DCONJG(A(K,I))*B(K,J) 120 CONTINUE IF (NOUNIT) TEMP = TEMP/DCONJG(A(I,I)) END IF B(I,J) = TEMP 130 CONTINUE 140 CONTINUE ELSE DO 180 J = 1,N DO 170 I = M,1,-1 TEMP = ALPHA*B(I,J) IF (NOCONJ) THEN DO 150 K = I + 1,M TEMP = TEMP - A(K,I)*B(K,J) 150 CONTINUE IF (NOUNIT) TEMP = TEMP/A(I,I) ELSE DO 160 K = I + 1,M TEMP = TEMP - DCONJG(A(K,I))*B(K,J) 160 CONTINUE IF (NOUNIT) TEMP = TEMP/DCONJG(A(I,I)) END IF B(I,J) = TEMP 170 CONTINUE 180 CONTINUE END IF END IF ELSE IF (LSAME(TRANSA,'N')) THEN ! ! Form B := alpha*B*inv( A ). ! IF (UPPER) THEN DO 230 J = 1,N IF (ALPHA.NE.ONE) THEN DO 190 I = 1,M B(I,J) = ALPHA*B(I,J) 190 CONTINUE END IF DO 210 K = 1,J - 1 IF (A(K,J).NE.ZERO) THEN DO 200 I = 1,M B(I,J) = B(I,J) - A(K,J)*B(I,K) 200 CONTINUE END IF 210 CONTINUE IF (NOUNIT) THEN TEMP = ONE/A(J,J) DO 220 I = 1,M B(I,J) = TEMP*B(I,J) 220 CONTINUE END IF 230 CONTINUE ELSE DO 280 J = N,1,-1 IF (ALPHA.NE.ONE) THEN DO 240 I = 1,M B(I,J) = ALPHA*B(I,J) 240 CONTINUE END IF DO 260 K = J + 1,N IF (A(K,J).NE.ZERO) THEN DO 250 I = 1,M B(I,J) = B(I,J) - A(K,J)*B(I,K) 250 CONTINUE END IF 260 CONTINUE IF (NOUNIT) THEN TEMP = ONE/A(J,J) DO 270 I = 1,M B(I,J) = TEMP*B(I,J) 270 CONTINUE END IF 280 CONTINUE END IF ELSE ! ! Form B := alpha*B*inv( A**T ) ! or B := alpha*B*inv( A**H ). ! IF (UPPER) THEN DO 330 K = N,1,-1 IF (NOUNIT) THEN IF (NOCONJ) THEN TEMP = ONE/A(K,K) ELSE TEMP = ONE/DCONJG(A(K,K)) END IF DO 290 I = 1,M B(I,K) = TEMP*B(I,K) 290 CONTINUE END IF DO 310 J = 1,K - 1 IF (A(J,K).NE.ZERO) THEN IF (NOCONJ) THEN TEMP = A(J,K) ELSE TEMP = DCONJG(A(J,K)) END IF DO 300 I = 1,M B(I,J) = B(I,J) - TEMP*B(I,K) 300 CONTINUE END IF 310 CONTINUE IF (ALPHA.NE.ONE) THEN DO 320 I = 1,M B(I,K) = ALPHA*B(I,K) 320 CONTINUE END IF 330 CONTINUE ELSE DO 380 K = 1,N IF (NOUNIT) THEN IF (NOCONJ) THEN TEMP = ONE/A(K,K) ELSE TEMP = ONE/DCONJG(A(K,K)) END IF DO 340 I = 1,M B(I,K) = TEMP*B(I,K) 340 CONTINUE END IF DO 360 J = K + 1,N IF (A(J,K).NE.ZERO) THEN IF (NOCONJ) THEN TEMP = A(J,K) ELSE TEMP = DCONJG(A(J,K)) END IF DO 350 I = 1,M B(I,J) = B(I,J) - TEMP*B(I,K) 350 CONTINUE END IF 360 CONTINUE IF (ALPHA.NE.ONE) THEN DO 370 I = 1,M B(I,K) = ALPHA*B(I,K) 370 CONTINUE END IF 380 CONTINUE END IF END IF END IF ! RETURN ! ! End of ZTRSM . ! END