subroutine sgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, & ldc ) !*****************************************************************************80 ! !! SGEMM computes C = alpha * A * B and related operations. ! ! Discussion: ! ! SGEMM 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 N by N matrix. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 13 February 2014 ! ! Author: ! ! Original FORTRAN77 version by Jack Dongarra. ! This version by John Burkardt. ! ! Parameters: ! ! Input, character * 1 TRANSA, specifies the form of op( A ) to be used in ! the matrix multiplication as follows: ! 'N' or 'n', op ( A ) = A. ! 'T' or 't', op ( A ) = A'. ! 'C' or 'c', op ( A ) = A'. ! ! Input, character * 1 TRANSB, specifies the form of op ( B ) to be used in ! the matrix multiplication as follows: ! 'N' or 'n', op ( B ) = B. ! 'T' or 't', op ( B ) = B'. ! 'C' or 'c', op ( B ) = B'. ! ! Input, integer ( kind = 4 ) M, the number of rows of the matrix op ( A ) ! and of the matrix C. 0 <= M. ! ! Input, integer ( kind = 4 ) N, the number of columns of the matrix ! op ( B ) and the number of columns of the matrix C. 0 <= N. ! ! Input, integer ( kind = 4 ) K, the number of columns of the matrix ! op ( A ) and the number of rows of the matrix op ( B ). 0 <= K. ! ! Input, real ( kind = 4 ) ALPHA, the scalar multiplier ! for op ( A ) * op ( B ). ! ! Input, real ( kind = 4 ) A(LDA,KA), where: ! if TRANSA is 'N' or 'n', KA is equal to K, and the leading M by K ! part of the array contains A; ! if TRANSA is not 'N' or 'n', then KA is equal to M, and the leading ! K by M part of the array must contain the matrix A. ! ! Input, integer ( kind = 4 ) LDA, the first dimension of A as declared ! in the calling routine. When TRANSA = 'N' or 'n' then LDA must be at ! least max ( 1, M ), otherwise LDA must be at least max ( 1, K ). ! ! Input, real ( kind = 4 ) B(LDB,KB), where: ! if TRANSB is 'N' or 'n', kB is N, and the leading K by N ! part of the array contains B; ! if TRANSB is not 'N' or 'n', then KB is equal to K, and the leading ! n by k part of the array must contain the matrix B. ! ! Input, integer ( kind = 4 ) LDB, the first dimension of B as declared in ! the calling routine. When TRANSB = 'N' or 'n' then LDB must be at least ! max ( 1, K ), otherwise LDB must be at least max ( 1, N ). ! ! Input, real ( kind = 4 ) BETA, the scalar multiplier for C. ! ! Input/output, real ( kind = 4 ) C(LDC,N). ! On input, the leading M by N part of this array must contain the ! matrix C, except when BETA is zero, in which case C need not be set. ! On output, the array C is overwritten by the M by N matrix ! alpha * op ( A ) * op ( B ) + beta * C. ! ! Input, integer ( kind = 4 ) LDC, the first dimension of C as declared ! in the calling routine. max ( 1, M ) <= LDC. ! implicit none integer ( kind = 4 ) lda integer ( kind = 4 ) ldb integer ( kind = 4 ) ldc real ( kind = 4 ) a(lda,*) real ( kind = 4 ) alpha real ( kind = 4 ) b(ldb,*) real ( kind = 4 ) beta real ( kind = 4 ) c(ldc,*) integer ( kind = 4 ) i integer ( kind = 4 ) info integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) ncola integer ( kind = 4 ) nrowa integer ( kind = 4 ) nrowb logical nota logical notb real ( kind = 4 ) temp character * 1 transa character * 1 transb ! ! 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 = ( transa == 'N' .or. transa == 'n' ) if ( nota ) then nrowa = m ncola = k else nrowa = k ncola = m end if notb = ( transb == 'N' .or. transb == 'n' ) if ( notb ) then nrowb = k else nrowb = n end if ! ! Test the input parameters. ! info = 0 if ( transa /= 'N' .and. transa /= 'n' .and. & transa /= 'C' .and. transa /= 'c' .and. & transa /= 'T' .and. transa /= 't' ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'SGEMM - Fatal error!' write ( *, '(a)' ) ' Input TRANSA had illegal value.' stop 1 end if if ( transb /= 'N' .and. transb /= 'n' .and. & transb /= 'C' .and. transb /= 'c' .and. & transb /= 'T' .and. transb /= 't' ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'SGEMM - Fatal error!' write ( *, '(a)' ) ' Input TRANSB had illegal value.' stop 1 end if if ( m < 0 ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'SGEMM - Fatal error!' write ( *, '(a)' ) ' Input M had illegal value.' stop 1 end if if ( n < 0 ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'SGEMM - Fatal error!' write ( *, '(a)' ) ' Input N had illegal value.' stop 1 end if if ( k < 0 ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'SGEMM - Fatal error!' write ( *, '(a)' ) ' Input K had illegal value.' stop 1 end if if ( lda < max ( 1, nrowa ) ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'SGEMM - Fatal error!' write ( *, '(a)' ) ' Input LDA had illegal value.' stop 1 end if if ( ldb < max ( 1, nrowb ) ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'SGEMM - Fatal error!' write ( *, '(a)' ) ' Input LDB had illegal value.' stop 1 end if if ( ldc < max ( 1, m ) ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'SGEMM - Fatal error!' write ( *, '(a)' ) ' Input LDC had illegal value.' stop 1 end if ! ! Quick return if possible. ! if ( m == 0 ) then return end if if ( n == 0 ) then return end if if ( ( alpha == 0.0E+00 .or. k == 0 ) .and. beta == 1.0E+00 ) then return end if ! ! And if alpha is zero. ! if ( alpha == 0.0E+00 ) then if ( beta == 0.0E+00 ) then c(1:m,1:n) = 0.0E+00 else c(1:m,1:n) = beta * c(1:m,1:n) end if return end if ! ! Start the operations. ! if ( notb ) then ! ! Form C := alpha*A*B + beta*C. ! if ( nota ) then do j = 1, n if ( beta == 0.0E+00 ) then c(1:m,j) = 0.0E+00 else if ( beta /= 1.0E+00 ) then c(1:m,j) = beta * c(1:m,j) end if do l = 1, k if ( b(l,j) /= 0.0E+00 ) then c(1:m,j) = c(1:m,j) + alpha * b(l,j) * a(1:m,l) end if end do end do ! ! Form C := alpha*A'*B + beta*C ! else do j = 1, n do i = 1, m temp = dot_product ( a(1:k,i), b(1:k,j) ) if ( beta == 0.0E+00 ) then c(i,j) = alpha * temp else c(i,j) = alpha * temp + beta * c(i,j) end if end do end do end if ! ! Form C := alpha*A*B' + beta*C ! else if ( nota ) then do j = 1, n if ( beta == 0.0E+00 ) then c(1:m,j) = 0.0E+00 else if ( beta /= 1.0E+00 ) then c(1:m,j) = beta * c(1:m,j) end if do l = 1, k if ( b(j,l) /= 0.0E+00 ) then c(1:m,j) = c(1:m,j) + alpha * b(j,l) * a(1:m,l) end if end do end do ! ! Form C := alpha*A'*B' + beta*C ! else do j = 1, n do i = 1, m temp = dot_product ( a(1:k,i), b(j,1:k) ) if ( beta == 0.0E+00 ) then c(i,j) = alpha * temp else c(i,j) = alpha * temp + beta * c(i,j) end if end do end do end if end if return end SUBROUTINE SSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC) !! SSYMM performs C:=alpha*A*B+beta*C, A symmetric, B and C rectangular. ! ! .. Scalar Arguments .. REAL ALPHA,BETA INTEGER LDA,LDB,LDC,M,N CHARACTER SIDE,UPLO ! .. ! .. Array Arguments .. REAL A(LDA,*),B(LDB,*),C(LDC,*) ! .. ! ! Purpose ! ======= ! ! SSYMM 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 - REAL . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - REAL 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 - REAL 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 - REAL . ! 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 - REAL 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 .. REAL TEMP1,TEMP2 INTEGER I,INFO,J,K,NROWA LOGICAL UPPER ! .. ! .. Parameters .. REAL ONE,ZERO PARAMETER (ONE=1.0E+0,ZERO=0.0E+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('SSYMM ',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 SSYMM . ! END SUBROUTINE SSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) !! SSYR2K performs C:=alpha*A*TRANSPOSE(B)+alpha*B*TRANSPOSE(A)+beta*C, C symmetric. ! ! .. Scalar Arguments .. REAL ALPHA,BETA INTEGER K,LDA,LDB,LDC,N CHARACTER TRANS,UPLO ! .. ! .. Array Arguments .. REAL A(LDA,*),B(LDB,*),C(LDC,*) ! .. ! ! Purpose ! ======= ! ! SSYR2K 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. ! ! TRANS = 'C' or 'c' 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' 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 - REAL . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - REAL 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 - REAL 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 - REAL . ! On entry, BETA specifies the scalar beta. ! Unchanged on exit. ! ! C - REAL 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 .. REAL TEMP1,TEMP2 INTEGER I,INFO,J,L,NROWA LOGICAL UPPER ! .. ! .. Parameters .. REAL ONE,ZERO PARAMETER (ONE=1.0E+0,ZERO=0.0E+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')) .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('SSYR2K',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 SSYR2K. ! END SUBROUTINE SSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC) !! SSYRK performs C:=alpha*A*TRANSPOSE(A)+beta*C, C symmetric. ! ! .. Scalar Arguments .. REAL ALPHA,BETA INTEGER K,LDA,LDC,N CHARACTER TRANS,UPLO ! .. ! .. Array Arguments .. REAL A(LDA,*),C(LDC,*) ! .. ! ! Purpose ! ======= ! ! SSYRK 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. ! ! TRANS = 'C' or 'c' 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' or 'C' or 'c', K specifies the number ! of rows of the matrix A. K must be at least zero. ! Unchanged on exit. ! ! ALPHA - REAL . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - REAL 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 - REAL . ! On entry, BETA specifies the scalar beta. ! Unchanged on exit. ! ! C - REAL 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 .. REAL TEMP INTEGER I,INFO,J,L,NROWA LOGICAL UPPER ! .. ! .. Parameters .. REAL ONE,ZERO PARAMETER (ONE=1.0E+0,ZERO=0.0E+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')) .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('SSYRK ',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 SSYRK . ! END SUBROUTINE STRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) !! STRMM performs B:=A*B or B:=B*A, A triangular, B rectangular. ! ! .. Scalar Arguments .. REAL ALPHA INTEGER LDA,LDB,M,N CHARACTER DIAG,SIDE,TRANSA,UPLO ! .. ! .. Array Arguments .. REAL A(LDA,*),B(LDB,*) ! .. ! ! Purpose ! ======= ! ! STRMM 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. ! ! 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**T. ! ! 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 - REAL . ! 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 - REAL 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 - REAL 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 MAX ! .. ! .. Local Scalars .. REAL TEMP INTEGER I,INFO,J,K,NROWA LOGICAL LSIDE,NOUNIT,UPPER ! .. ! .. Parameters .. REAL ONE,ZERO PARAMETER (ONE=1.0E+0,ZERO=0.0E+0) ! .. ! ! 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.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('STRMM ',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. ! IF (UPPER) THEN DO 110 J = 1,N DO 100 I = M,1,-1 TEMP = B(I,J) IF (NOUNIT) TEMP = TEMP*A(I,I) 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) TEMP = TEMP*A(I,I) 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) TEMP = TEMP*A(J,J) 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) TEMP = TEMP*A(J,J) 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**T. ! 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) TEMP = TEMP*A(K,K) 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) TEMP = TEMP*A(K,K) 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 of STRMM . ! END SUBROUTINE STRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) !! STRSM performs B:=INVERSE(A)*C or B:=C*INVERSE(A), B and C rectangular, A triangular. ! ! .. Scalar Arguments .. REAL ALPHA INTEGER LDA,LDB,M,N CHARACTER DIAG,SIDE,TRANSA,UPLO ! .. ! .. Array Arguments .. REAL A(LDA,*),B(LDB,*) ! .. ! ! Purpose ! ======= ! ! STRSM 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. ! ! 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**T. ! ! 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 - REAL . ! 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 - REAL 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 - REAL 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 MAX ! .. ! .. Local Scalars .. REAL TEMP INTEGER I,INFO,J,K,NROWA LOGICAL LSIDE,NOUNIT,UPPER ! .. ! .. Parameters .. REAL ONE,ZERO PARAMETER (ONE=1.0E+0,ZERO=0.0E+0) ! .. ! ! 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.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('STRSM ',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. ! 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) TEMP = TEMP/A(I,I) 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) TEMP = TEMP/A(I,I) 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**T ). ! 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 of STRSM . ! END