Mercurial > octave-nkf
changeset 4347:024ef171aec3
[project @ 2003-02-20 23:31:46 by jwe]
| author | jwe |
|---|---|
| date | Thu, 20 Feb 2003 23:32:25 +0000 |
| parents | d39de791ef9c |
| children | 05415e529cef |
| files | libcruft/ChangeLog libcruft/blas/sgemm.f libcruft/blas/sgemv.f libcruft/blas/sscal.f libcruft/blas/ssyrk.f libcruft/blas/strsm.f |
| diffstat | 6 files changed, 1294 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Thu Feb 20 21:38:39 2003 +0000 +++ b/libcruft/ChangeLog Thu Feb 20 23:32:25 2003 +0000 @@ -1,3 +1,8 @@ +2003-02-20 John W. Eaton <jwe@bevo.che.wisc.edu> + + * blas/sgemm.f, blas/strsm.f, blas/ssyrk.f, blas/sscal.f, + blas/sgemv.f, blas/sdot.f: New files. + 2003-02-20 Paul Kienzle <pkienzle@users.sf.net> * dassl/ddaslv.f: Fortran doesn't use ; in statements.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/sgemm.f Thu Feb 20 23:32:25 2003 +0000 @@ -0,0 +1,313 @@ + SUBROUTINE SGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, + $ BETA, C, LDC ) +* .. Scalar Arguments .. + CHARACTER*1 TRANSA, TRANSB + INTEGER M, N, K, LDA, LDB, LDC + REAL ALPHA, BETA +* .. Array Arguments .. + REAL A( LDA, * ), B( LDB, * ), C( LDC, * ) +* .. +* +* Purpose +* ======= +* +* 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 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 - REAL . +* On entry, ALPHA specifies the scalar alpha. +* Unchanged on exit. +* +* A - REAL 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 - REAL 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 - 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 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 + REAL TEMP +* .. Parameters .. + REAL ONE , ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+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 .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( 'SGEMM ', 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 if 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 +* +* 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.EQ.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.EQ.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.EQ.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 of SGEMM . +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/sgemv.f Thu Feb 20 23:32:25 2003 +0000 @@ -0,0 +1,261 @@ + SUBROUTINE SGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, + $ BETA, Y, INCY ) +* .. Scalar Arguments .. + REAL ALPHA, BETA + INTEGER INCX, INCY, LDA, M, N + CHARACTER*1 TRANS +* .. Array Arguments .. + REAL A( LDA, * ), X( * ), Y( * ) +* .. +* +* Purpose +* ======= +* +* SGEMV performs one of the matrix-vector operations +* +* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, +* +* where alpha and beta are scalars, x and y are vectors and A is an +* m by n matrix. +* +* Parameters +* ========== +* +* TRANS - CHARACTER*1. +* On entry, TRANS specifies the operation to be performed as +* follows: +* +* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. +* +* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. +* +* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. +* +* Unchanged on exit. +* +* M - INTEGER. +* On entry, M specifies the number of rows of the matrix A. +* M must be at least zero. +* Unchanged on exit. +* +* N - INTEGER. +* On entry, N specifies the number of columns of the matrix A. +* N must be at least zero. +* Unchanged on exit. +* +* ALPHA - REAL . +* On entry, ALPHA specifies the scalar alpha. +* Unchanged on exit. +* +* A - REAL array of DIMENSION ( LDA, n ). +* Before entry, the leading m by n part of the array A must +* contain the matrix of coefficients. +* Unchanged on exit. +* +* LDA - INTEGER. +* On entry, LDA specifies the first dimension of A as declared +* in the calling (sub) program. LDA must be at least +* max( 1, m ). +* Unchanged on exit. +* +* X - REAL array of DIMENSION at least +* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' +* and at least +* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. +* Before entry, the incremented array X must contain the +* vector x. +* Unchanged on exit. +* +* INCX - INTEGER. +* On entry, INCX specifies the increment for the elements of +* X. INCX must not be zero. +* Unchanged on exit. +* +* BETA - REAL . +* On entry, BETA specifies the scalar beta. When BETA is +* supplied as zero then Y need not be set on input. +* Unchanged on exit. +* +* Y - REAL array of DIMENSION at least +* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' +* and at least +* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. +* Before entry with BETA non-zero, the incremented array Y +* must contain the vector y. On exit, Y is overwritten by the +* updated vector y. +* +* INCY - INTEGER. +* On entry, INCY specifies the increment for the elements of +* Y. INCY must not be zero. +* Unchanged on exit. +* +* +* Level 2 Blas routine. +* +* -- Written on 22-October-1986. +* Jack Dongarra, Argonne National Lab. +* Jeremy Du Croz, Nag Central Office. +* Sven Hammarling, Nag Central Office. +* Richard Hanson, Sandia National Labs. +* +* +* .. Parameters .. + REAL ONE , ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) +* .. Local Scalars .. + REAL TEMP + INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. External Subroutines .. + EXTERNAL XERBLA +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + IF ( .NOT.LSAME( TRANS, 'N' ).AND. + $ .NOT.LSAME( TRANS, 'T' ).AND. + $ .NOT.LSAME( TRANS, 'C' ) )THEN + INFO = 1 + ELSE IF( M.LT.0 )THEN + INFO = 2 + ELSE IF( N.LT.0 )THEN + INFO = 3 + ELSE IF( LDA.LT.MAX( 1, M ) )THEN + INFO = 6 + ELSE IF( INCX.EQ.0 )THEN + INFO = 8 + ELSE IF( INCY.EQ.0 )THEN + INFO = 11 + END IF + IF( INFO.NE.0 )THEN + CALL XERBLA( 'SGEMV ', 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 +* +* Set LENX and LENY, the lengths of the vectors x and y, and set +* up the start points in X and Y. +* + IF( LSAME( TRANS, 'N' ) )THEN + LENX = N + LENY = M + ELSE + LENX = M + LENY = N + END IF + IF( INCX.GT.0 )THEN + KX = 1 + ELSE + KX = 1 - ( LENX - 1 )*INCX + END IF + IF( INCY.GT.0 )THEN + KY = 1 + ELSE + KY = 1 - ( LENY - 1 )*INCY + END IF +* +* Start the operations. In this version the elements of A are +* accessed sequentially with one pass through A. +* +* First form y := beta*y. +* + IF( BETA.NE.ONE )THEN + IF( INCY.EQ.1 )THEN + IF( BETA.EQ.ZERO )THEN + DO 10, I = 1, LENY + Y( I ) = ZERO + 10 CONTINUE + ELSE + DO 20, I = 1, LENY + Y( I ) = BETA*Y( I ) + 20 CONTINUE + END IF + ELSE + IY = KY + IF( BETA.EQ.ZERO )THEN + DO 30, I = 1, LENY + Y( IY ) = ZERO + IY = IY + INCY + 30 CONTINUE + ELSE + DO 40, I = 1, LENY + Y( IY ) = BETA*Y( IY ) + IY = IY + INCY + 40 CONTINUE + END IF + END IF + END IF + IF( ALPHA.EQ.ZERO ) + $ RETURN + IF( LSAME( TRANS, 'N' ) )THEN +* +* Form y := alpha*A*x + y. +* + JX = KX + IF( INCY.EQ.1 )THEN + DO 60, J = 1, N + IF( X( JX ).NE.ZERO )THEN + TEMP = ALPHA*X( JX ) + DO 50, I = 1, M + Y( I ) = Y( I ) + TEMP*A( I, J ) + 50 CONTINUE + END IF + JX = JX + INCX + 60 CONTINUE + ELSE + DO 80, J = 1, N + IF( X( JX ).NE.ZERO )THEN + TEMP = ALPHA*X( JX ) + IY = KY + DO 70, I = 1, M + Y( IY ) = Y( IY ) + TEMP*A( I, J ) + IY = IY + INCY + 70 CONTINUE + END IF + JX = JX + INCX + 80 CONTINUE + END IF + ELSE +* +* Form y := alpha*A'*x + y. +* + JY = KY + IF( INCX.EQ.1 )THEN + DO 100, J = 1, N + TEMP = ZERO + DO 90, I = 1, M + TEMP = TEMP + A( I, J )*X( I ) + 90 CONTINUE + Y( JY ) = Y( JY ) + ALPHA*TEMP + JY = JY + INCY + 100 CONTINUE + ELSE + DO 120, J = 1, N + TEMP = ZERO + IX = KX + DO 110, I = 1, M + TEMP = TEMP + A( I, J )*X( IX ) + IX = IX + INCX + 110 CONTINUE + Y( JY ) = Y( JY ) + ALPHA*TEMP + JY = JY + INCY + 120 CONTINUE + END IF + END IF +* + RETURN +* +* End of SGEMV . +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/sscal.f Thu Feb 20 23:32:25 2003 +0000 @@ -0,0 +1,43 @@ + subroutine sscal(n,sa,sx,incx) +c +c scales a vector by a constant. +c uses unrolled loops for increment equal to 1. +c jack dongarra, linpack, 3/11/78. +c modified 3/93 to return if incx .le. 0. +c modified 12/3/93, array(1) declarations changed to array(*) +c + real sa,sx(*) + integer i,incx,m,mp1,n,nincx +c + if( n.le.0 .or. incx.le.0 )return + if(incx.eq.1)go to 20 +c +c code for increment not equal to 1 +c + nincx = n*incx + do 10 i = 1,nincx,incx + sx(i) = sa*sx(i) + 10 continue + return +c +c code for increment equal to 1 +c +c +c clean-up loop +c + 20 m = mod(n,5) + if( m .eq. 0 ) go to 40 + do 30 i = 1,m + sx(i) = sa*sx(i) + 30 continue + if( n .lt. 5 ) return + 40 mp1 = m + 1 + do 50 i = mp1,n,5 + sx(i) = sa*sx(i) + sx(i + 1) = sa*sx(i + 1) + sx(i + 2) = sa*sx(i + 2) + sx(i + 3) = sa*sx(i + 3) + sx(i + 4) = sa*sx(i + 4) + 50 continue + return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/ssyrk.f Thu Feb 20 23:32:25 2003 +0000 @@ -0,0 +1,294 @@ + SUBROUTINE SSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, + $ BETA, C, LDC ) +* .. Scalar Arguments .. + CHARACTER*1 UPLO, TRANS + INTEGER N, K, LDA, LDC + REAL ALPHA, BETA +* .. Array Arguments .. + REAL A( LDA, * ), C( LDC, * ) +* .. +* +* Purpose +* ======= +* +* SSYRK 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 - 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. +* +* +* 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 + REAL TEMP +* .. Parameters .. + REAL ONE , ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+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 .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' + 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'*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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/strsm.f Thu Feb 20 23:32:25 2003 +0000 @@ -0,0 +1,378 @@ + SUBROUTINE STRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, + $ B, LDB ) +* .. Scalar Arguments .. + CHARACTER*1 SIDE, UPLO, TRANSA, DIAG + INTEGER M, N, LDA, LDB + REAL ALPHA +* .. 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'. +* +* 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 - 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. +* +* +* 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 + REAL TEMP +* .. Parameters .. + REAL ONE , ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+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 .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( 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' )*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' ). +* + 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