subroutine dtrmm ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB ) ! !******************************************************************************* ! !! DTRMM ! character SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB DOUBLE PRECISION ALPHA ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) ! .. ! ! Purpose ! ======= ! ! DTRMM 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'. ! ! Parameters ! ========== ! ! SIDE - CHARACTER*1. ! On entry, SIDE specifies whether op( A ) multiplies B from ! the left or right as follows: ! ! SIDE = 'L' or 'l' B := alpha*op( A )*B. ! ! SIDE = 'R' or 'r' B := alpha*B*op( A ). ! ! Unchanged on exit. ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix A is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANSA - CHARACTER*1. ! On entry, TRANSA specifies the form of op( A ) to be used in ! the matrix multiplication as follows: ! ! TRANSA = 'N' or 'n' op( A ) = A. ! ! TRANSA = 'T' or 't' op( A ) = A'. ! ! TRANSA = 'C' or 'c' op( A ) = A'. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit triangular ! as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of B. M must be at ! least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of B. N must be ! at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION array of DIMENSION ( LDB, n ). ! Before entry, the leading m by n part of the array B must ! contain the matrix B, and on exit is overwritten by the ! transformed matrix. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. LDB must be at least ! max( 1, m ). ! Unchanged on exit. ! ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! LSIDE = LSAME( SIDE , 'L' ) if ( LSIDE ) then NROWA = M else NROWA = N end if 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 < 0 ) then INFO = 5 else if ( N < 0 ) then INFO = 6 else if ( LDA < MAX( 1, NROWA ) ) then INFO = 9 else if ( LDB < MAX( 1, M ) ) then INFO = 11 end if if ( INFO.NE.0 ) then CALL XERBLA( 'DTRMM ', INFO ) return end if ! ! Quick return if possible. ! if ( N == 0 ) then return end if ! ! And when alpha.eq.zero. ! if ( ALPHA == 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 ) then TEMP = TEMP*A( K, K ) end if 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 ) then B( K, J ) = B( K, J )*A( K, K ) end if DO 60, I = K + 1, M B( I, J ) = B( I, J ) + TEMP*A( I, K ) 60 CONTINUE end if 70 CONTINUE 80 CONTINUE end if else ! ! Form B := alpha*A'*B. ! if ( UPPER ) then DO 110, J = 1, N DO 100, I = M, 1, -1 TEMP = B( I, J ) if ( NOUNIT ) then TEMP = TEMP*A( I, I ) end if DO 90, K = 1, I - 1 TEMP = TEMP + A( K, I )*B( K, J ) 90 CONTINUE B( I, J ) = ALPHA*TEMP 100 CONTINUE 110 CONTINUE else DO 140, J = 1, N DO 130, I = 1, M TEMP = B( I, J ) if ( NOUNIT ) then TEMP = TEMP*A( I, I ) end if DO 120, K = I + 1, M TEMP = TEMP + A( K, I )*B( K, J ) 120 CONTINUE B( I, J ) = ALPHA*TEMP 130 CONTINUE 140 CONTINUE end if end if else if ( LSAME( TRANSA, 'N' ) ) then ! ! Form B := alpha*B*A. ! if ( UPPER ) then DO 180, J = N, 1, -1 TEMP = ALPHA if ( NOUNIT ) then TEMP = TEMP*A( J, J ) end if DO 150, I = 1, M B( I, J ) = TEMP*B( I, J ) 150 CONTINUE DO 170, K = 1, J - 1 if ( A( K, J ).NE.ZERO ) then TEMP = ALPHA*A( K, J ) DO 160, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 160 CONTINUE end if 170 CONTINUE 180 CONTINUE else DO 220, J = 1, N TEMP = ALPHA if ( NOUNIT ) then TEMP = TEMP*A( J, J ) end if DO 190, I = 1, M B( I, J ) = TEMP*B( I, J ) 190 CONTINUE DO 210, K = J + 1, N if ( A( K, J ).NE.ZERO ) then TEMP = ALPHA*A( K, J ) DO 200, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 200 CONTINUE end if 210 CONTINUE 220 CONTINUE end if else ! ! Form B := alpha*B*A'. ! if ( UPPER ) then DO 260, K = 1, N DO 240, J = 1, K - 1 if ( A( J, K ).NE.ZERO ) then TEMP = ALPHA*A( J, K ) DO 230, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 230 CONTINUE end if 240 CONTINUE TEMP = ALPHA if ( NOUNIT ) then TEMP = TEMP*A( K, K ) end if if ( TEMP.NE.ONE ) then DO 250, I = 1, M B( I, K ) = TEMP*B( I, K ) 250 CONTINUE end if 260 CONTINUE else DO 300, K = N, 1, -1 DO 280, J = K + 1, N if ( A( J, K ).NE.ZERO ) then TEMP = ALPHA*A( J, K ) DO 270, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 270 CONTINUE end if 280 CONTINUE TEMP = ALPHA if ( NOUNIT ) then TEMP = TEMP*A( K, K ) end if if ( TEMP.NE.ONE ) then DO 290, I = 1, M B( I, K ) = TEMP*B( I, K ) 290 CONTINUE end if 300 CONTINUE end if end if end if return end subroutine dtrsm ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB ) ! !******************************************************************************* ! !! DTRSM ! character SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB DOUBLE PRECISION ALPHA ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) ! .. ! ! Purpose ! ======= ! ! DTRSM 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'. ! ! The matrix X is overwritten on B. ! ! Parameters ! ========== ! ! SIDE - CHARACTER*1. ! On entry, SIDE specifies whether op( A ) appears on the left ! or right of X as follows: ! ! SIDE = 'L' or 'l' op( A )*X = alpha*B. ! ! SIDE = 'R' or 'r' X*op( A ) = alpha*B. ! ! Unchanged on exit. ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix A is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANSA - CHARACTER*1. ! On entry, TRANSA specifies the form of op( A ) to be used in ! the matrix multiplication as follows: ! ! TRANSA = 'N' or 'n' op( A ) = A. ! ! TRANSA = 'T' or 't' op( A ) = A'. ! ! TRANSA = 'C' or 'c' op( A ) = A'. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit triangular ! as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of B. M must be at ! least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of B. N must be ! at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION array of DIMENSION ( LDB, n ). ! Before entry, the leading m by n part of the array B must ! contain the right-hand side matrix B, and on exit is ! overwritten by the solution matrix X. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. LDB must be at least ! max( 1, m ). ! Unchanged on exit. ! ! ! Level 3 Blas routine. ! ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! LSIDE = LSAME( SIDE , 'L' ) if ( LSIDE ) then NROWA = M else NROWA = N end if 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 < 0 ) then INFO = 5 else if ( N < 0 ) then INFO = 6 else if ( LDA < MAX( 1, NROWA ) ) then INFO = 9 else if ( LDB < MAX( 1, M ) ) then INFO = 11 end if if ( INFO.NE.0 ) then CALL XERBLA( 'DTRSM ', INFO ) return end if ! ! Quick return if possible. ! if ( N == 0 ) then return end if ! ! And when alpha.eq.zero. ! if ( ALPHA == 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 ) then B( K, J ) = B( K, J )/A( K, K ) end if 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 ) then B( K, J ) = B( K, J )/A( K, K ) end if DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE end if 90 CONTINUE 100 CONTINUE end if else ! ! Form B := alpha*inv( A' )*B. ! if ( UPPER ) then DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE if ( NOUNIT ) then TEMP = TEMP/A( I, I ) end if B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE else DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE if ( NOUNIT ) then TEMP = TEMP/A( I, I ) end if B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE end if end if else if ( LSAME( TRANSA, 'N' ) ) then ! ! Form B := alpha*B*inv( A ). ! if ( UPPER ) then DO 210, J = 1, N if ( ALPHA.NE.ONE ) then DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE end if DO 190, K = 1, J - 1 if ( A( K, J ).NE.ZERO ) then DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE end if 190 CONTINUE if ( NOUNIT ) then TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE end if 210 CONTINUE else DO 260, J = N, 1, -1 if ( ALPHA.NE.ONE ) then DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE end if DO 240, K = J + 1, N if ( A( K, J ).NE.ZERO ) then DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE end if 240 CONTINUE if ( NOUNIT ) then TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE end if 260 CONTINUE end if else ! ! Form B := alpha*B*inv( A' ). ! if ( UPPER ) then DO 310, K = N, 1, -1 if ( NOUNIT ) then TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE end if DO 290, J = 1, K - 1 if ( A( J, K ).NE.ZERO ) then TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE end if 290 CONTINUE if ( ALPHA.NE.ONE ) then DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE end if 310 CONTINUE else DO 360, K = 1, N if ( NOUNIT ) then TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE end if DO 340, J = K + 1, N if ( A( J, K ).NE.ZERO ) then TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE end if 340 CONTINUE if ( ALPHA.NE.ONE ) then DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE end if 360 CONTINUE end if end if end if return end subroutine dsymm ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, BETA, C, LDC ) ! !******************************************************************************* ! !! DSYMM ! character SIDE, UPLO INTEGER M, N, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) ! .. ! ! Purpose ! ======= ! ! DSYMM performs one of the matrix-matrix operations ! ! C := alpha*A*B + beta*C, ! ! or ! ! C := alpha*B*A + beta*C, ! ! where alpha and beta are scalars, A is a symmetric matrix and B and ! C are m by n matrices. ! ! Parameters ! ========== ! ! SIDE - CHARACTER*1. ! On entry, SIDE specifies whether the symmetric matrix A ! appears on the left or right in the operation as follows: ! ! SIDE = 'L' or 'l' C := alpha*A*B + beta*C, ! ! SIDE = 'R' or 'r' C := alpha*B*A + beta*C, ! ! Unchanged on exit. ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the symmetric matrix A is to be ! referenced as follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of the ! symmetric matrix is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of the ! symmetric matrix is to be referenced. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix C. ! M must be at least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix C. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION. ! 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 - DOUBLE PRECISION array of DIMENSION ( LDC, n ). ! Before entry, the leading m by n part of the array C must ! contain the matrix C, except when beta is zero, in which ! case C need not be set on entry. ! On exit, the array C is overwritten by the m by n updated ! matrix. ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, m ). ! Unchanged on exit. ! ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP1, TEMP2 ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Executable Statements .. ! ! Set NROWA as the number of rows of A. ! if ( LSAME( SIDE, 'L' ) ) then NROWA = M else NROWA = N end if UPPER = LSAME( UPLO, 'U' ) ! ! Test the input parameters. ! INFO = 0 if ( ( .NOT.LSAME( SIDE, 'L' ) ) .and. & ( .NOT.LSAME( SIDE, 'R' ) ) ) then INFO = 1 else if ( ( .NOT.UPPER ) .and. & ( .NOT.LSAME( UPLO, 'L' ) ) ) then INFO = 2 else if ( M < 0 ) then INFO = 3 else if ( N < 0 ) then INFO = 4 else if ( LDA < MAX( 1, NROWA ) ) then INFO = 7 else if ( LDB < MAX( 1, M ) ) then INFO = 9 else if ( LDC < MAX( 1, M ) ) then INFO = 12 end if if ( INFO.NE.0 ) then CALL XERBLA( 'DSYMM ', INFO ) return end if ! ! Quick return if possible. ! if ( ( M == 0 ) .or. ( N == 0 ).OR. & ( ( ALPHA == ZERO ) .and. ( BETA == ONE ) ) ) then return end if ! ! And when alpha.eq.zero. ! if ( ALPHA == ZERO ) then if ( BETA == 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 == 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 == 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 == 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 subroutine dsyrk ( UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC ) ! !******************************************************************************* ! !! DSYRK ! character UPLO, TRANS INTEGER N, K, LDA, LDC DOUBLE PRECISION ALPHA, BETA ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ) ! .. ! ! Purpose ! ======= ! ! DSYRK performs one of the symmetric rank k operations ! ! C := alpha*A*A' + beta*C, ! ! or ! ! C := alpha*A'*A + beta*C, ! ! where alpha and beta are scalars, C is an n by n symmetric matrix ! and A is an n by k matrix in the first case and a k by n matrix ! in the second case. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array C is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of C ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of C ! is to be referenced. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. ! ! TRANS = 'T' or 't' C := alpha*A'*A + beta*C. ! ! TRANS = 'C' or 'c' C := alpha*A'*A + beta*C. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry with TRANS = 'N' or 'n', K specifies the number ! of columns of the matrix A, and on entry with ! TRANS = 'T' or 't' or '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 - DOUBLE PRECISION 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 - DOUBLE PRECISION array of DIMENSION ( LDC, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array C must contain the upper ! triangular part of the symmetric matrix and the strictly ! lower triangular part of C is not referenced. On exit, the ! upper triangular part of the array C is overwritten by the ! upper triangular part of the updated matrix. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array C must contain the lower ! triangular part of the symmetric matrix and the strictly ! upper triangular part of C is not referenced. On exit, the ! lower triangular part of the array C is overwritten by the ! lower triangular part of the updated matrix. ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, n ). ! Unchanged on exit. ! ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA DOUBLE PRECISION TEMP ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! if ( LSAME( TRANS, 'N' ) ) then NROWA = N else NROWA = K end if UPPER = LSAME( UPLO, 'U' ) ! INFO = 0 if ( ( .NOT.UPPER ) .and. & ( .NOT.LSAME( UPLO , 'L' ) ) ) then INFO = 1 else if ( ( .NOT.LSAME( TRANS, 'N' ) ) .and. & ( .NOT.LSAME( TRANS, 'T' ) ) .and. & ( .NOT.LSAME( TRANS, 'C' ) ) ) then INFO = 2 else if ( N < 0 ) then INFO = 3 else if ( K < 0 ) then INFO = 4 else if ( LDA < MAX( 1, NROWA ) ) then INFO = 7 else if ( LDC < MAX( 1, N ) ) then INFO = 10 end if if ( INFO.NE.0 ) then CALL XERBLA( 'DSYRK ', INFO ) return end if ! ! Quick return if possible. ! if ( ( N == 0 ) .or. & ( ( ( ALPHA == ZERO ) .or. ( K == 0 ) ) .and. ( BETA == ONE ) ) ) then return end if ! ! And when alpha.eq.zero. ! if ( ALPHA == ZERO ) then if ( UPPER ) then if ( BETA == 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 == ZERO ) then DO 60, J = 1, N DO 50, I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE else DO 80, J = 1, N DO 70, I = J, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE end if end if return end if ! ! Start the operations. ! if ( LSAME( TRANS, 'N' ) ) then ! ! Form C := alpha*A*A' + beta*C. ! if ( UPPER ) then DO 130, J = 1, N if ( BETA == 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 == ZERO ) then DO 140, I = J, N C( I, J ) = ZERO 140 CONTINUE else if ( BETA.NE.ONE ) then DO 150, I = J, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE end if DO 170, L = 1, K if ( A( J, L ).NE.ZERO ) then TEMP = ALPHA*A( J, L ) DO 160, I = J, N C( I, J ) = C( I, J ) + TEMP*A( I, L ) 160 CONTINUE end if 170 CONTINUE 180 CONTINUE end if else ! ! Form C := alpha*A'*A + beta*C. ! if ( UPPER ) then DO 210, J = 1, N DO 200, I = 1, J TEMP = ZERO DO 190, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 190 CONTINUE if ( BETA == 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 == 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 subroutine dsyr2k ( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC ) ! !******************************************************************************* ! !! DSYR2K ! character UPLO, TRANS INTEGER N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) ! .. ! ! Purpose ! ======= ! ! DSYR2K performs one of the symmetric rank 2k operations ! ! C := alpha*A*B' + alpha*B*A' + beta*C, ! ! or ! ! C := alpha*A'*B + alpha*B'*A + beta*C, ! ! where alpha and beta are scalars, C is an n by n symmetric matrix ! and A and B are n by k matrices in the first case and k by n ! matrices in the second case. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array C is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of C ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of C ! is to be referenced. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' C := alpha*A*B' + alpha*B*A' + ! beta*C. ! ! TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A + ! beta*C. ! ! TRANS = 'C' or 'c' C := alpha*A'*B + alpha*B'*A + ! beta*C. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry with TRANS = 'N' or 'n', K specifies the number ! of columns of the matrices A and B, and on entry with ! TRANS = 'T' or 't' or '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 - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION array of DIMENSION ( LDC, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array C must contain the upper ! triangular part of the symmetric matrix and the strictly ! lower triangular part of C is not referenced. On exit, the ! upper triangular part of the array C is overwritten by the ! upper triangular part of the updated matrix. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array C must contain the lower ! triangular part of the symmetric matrix and the strictly ! upper triangular part of C is not referenced. On exit, the ! lower triangular part of the array C is overwritten by the ! lower triangular part of the updated matrix. ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, n ). ! Unchanged on exit. ! ! ! Level 3 Blas routine. ! ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA DOUBLE PRECISION TEMP1, TEMP2 ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! if ( LSAME( TRANS, 'N' ) ) then NROWA = N else NROWA = K end if UPPER = LSAME( UPLO, 'U' ) ! INFO = 0 if ( ( .NOT.UPPER ) .and. & ( .NOT.LSAME( UPLO , 'L' ) ) ) then INFO = 1 else if ( ( .NOT.LSAME( TRANS, 'N' ) ) .and. & ( .NOT.LSAME( TRANS, 'T' ) ) .and. & ( .NOT.LSAME( TRANS, 'C' ) ) ) then INFO = 2 else if ( N < 0 ) then INFO = 3 else if ( K < 0 ) then INFO = 4 else if ( LDA < MAX( 1, NROWA ) ) then INFO = 7 else if ( LDB < MAX( 1, NROWA ) ) then INFO = 9 else if ( LDC < MAX( 1, N ) ) then INFO = 12 end if if ( INFO.NE.0 ) then CALL XERBLA( 'DSYR2K', INFO ) return end if ! ! Quick return if possible. ! if ( ( N == 0 ) .or. & ( ( ( ALPHA == ZERO ) .or. ( K == 0 ) ) .and. ( BETA == ONE ) ) ) then return end if ! ! And when alpha.eq.zero. ! if ( ALPHA == ZERO ) then if ( UPPER ) then if ( BETA == 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 == ZERO ) then DO 60, J = 1, N DO 50, I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE else DO 80, J = 1, N DO 70, I = J, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE end if end if return end if ! ! Start the operations. ! if ( LSAME( TRANS, 'N' ) ) then ! ! Form C := alpha*A*B' + alpha*B*A' + C. ! if ( UPPER ) then DO 130, J = 1, N if ( BETA == 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 == ZERO ) then DO 140, I = J, N C( I, J ) = ZERO 140 CONTINUE else if ( BETA.NE.ONE ) then DO 150, I = J, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE end if DO 170, L = 1, K if ( ( A( J, L ).NE.ZERO ) .or. ( B( J, L ).NE.ZERO ) ) then TEMP1 = ALPHA*B( J, L ) TEMP2 = ALPHA*A( J, L ) DO 160, I = J, N C( I, J ) = C( I, J ) + A( I, L )*TEMP1 + B( I, L )*TEMP2 160 CONTINUE end if 170 CONTINUE 180 CONTINUE end if else ! ! Form C := alpha*A'*B + alpha*B'*A + C. ! if ( UPPER ) then DO 210, J = 1, N DO 200, I = 1, J TEMP1 = ZERO TEMP2 = ZERO DO 190, L = 1, K TEMP1 = TEMP1 + A( L, I )*B( L, J ) TEMP2 = TEMP2 + B( L, I )*A( L, J ) 190 CONTINUE if ( BETA == 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 == 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 subroutine dgemm ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC ) ! !******************************************************************************* ! !! DGEMM ! character TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) ! .. ! ! Purpose ! ======= ! ! DGEMM 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', ! ! alpha and beta are scalars, and A, B and C are matrices, with op( A ) ! an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. ! ! Parameters ! ========== ! ! TRANSA - CHARACTER*1. ! On entry, TRANSA specifies the form of op( A ) to be used in ! the matrix multiplication as follows: ! ! TRANSA = 'N' or 'n', op( A ) = A. ! ! TRANSA = 'T' or 't', op( A ) = A'. ! ! TRANSA = 'C' or 'c', op( A ) = A'. ! ! Unchanged on exit. ! ! TRANSB - CHARACTER*1. ! On entry, TRANSB specifies the form of op( B ) to be used in ! the matrix multiplication as follows: ! ! TRANSB = 'N' or 'n', op( B ) = B. ! ! TRANSB = 'T' or 't', op( B ) = B'. ! ! TRANSB = 'C' or 'c', op( B ) = B'. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix ! op( A ) and of the matrix C. M must be at least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix ! op( B ) and the number of columns of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry, K specifies the number of columns of the matrix ! op( A ) and the number of rows of the matrix op( B ). K must ! be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION. ! 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 - DOUBLE PRECISION array of DIMENSION ( LDC, n ). ! Before entry, the leading m by n part of the array C must ! contain the matrix C, except when beta is zero, in which ! case C need not be set on entry. ! On exit, the array C is overwritten by the m by n matrix ! ( alpha*op( A )*op( B ) + beta*C ). ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, m ). ! Unchanged on exit. ! ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB DOUBLE PRECISION TEMP ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Executable Statements .. ! ! Set NOTA and NOTB as true if A and B respectively are not ! transposed 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' ) 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.LSAME( TRANSA, 'C' ) ) .and. & ( .NOT.LSAME( TRANSA, 'T' ) ) ) then INFO = 1 else if ( ( .NOT.NOTB ) .and. & ( .NOT.LSAME( TRANSB, 'C' ) ) .and. & ( .NOT.LSAME( TRANSB, 'T' ) ) ) then INFO = 2 else if ( M < 0 ) then INFO = 3 else if ( N < 0 ) then INFO = 4 else if ( K < 0 ) then INFO = 5 else if ( LDA < MAX( 1, NROWA ) ) then INFO = 8 else if ( LDB < MAX( 1, NROWB ) ) then INFO = 10 else if ( LDC < MAX( 1, M ) ) then INFO = 13 end if if ( INFO.NE.0 ) then CALL XERBLA( 'DGEMM ', INFO ) return end if ! ! Quick return if possible. ! if ( ( M == 0 ) .or. ( N == 0 ).OR. & ( ( ( ALPHA == ZERO ) .or. ( K == 0 ) ) .and. ( BETA == ONE ) ) ) then return end if ! ! And if alpha.eq.zero. ! if ( ALPHA == ZERO ) then if ( BETA == 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 == 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 ! ! Form C := alpha*A'*B + beta*C ! DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE if ( BETA == ZERO ) then C( I, J ) = ALPHA*TEMP else C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) end if 110 CONTINUE 120 CONTINUE end if else if ( NOTA ) then ! ! Form C := alpha*A*B' + beta*C ! DO 170, J = 1, N if ( BETA == ZERO ) then DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE else if ( BETA.NE.ONE ) then DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE end if DO 160, L = 1, K if ( B( J, L ).NE.ZERO ) then TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE end if 160 CONTINUE 170 CONTINUE else ! ! Form C := alpha*A'*B' + beta*C ! DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE if ( BETA == ZERO ) then C( I, J ) = ALPHA*TEMP else C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) end if 190 CONTINUE 200 CONTINUE end if end if return end