SUBROUTINE MVMTLS( TRANS, N, M, X, LDX, Y, LDY ) * .. * .. Scalar Arugments .. INTEGER N, M, LDX, LDY, TRANS DOUBLE PRECISION X( LDX, * ), Y( LDY, * ) * .. * .. Common blocks .. DOUBLE PRECISION BLK( 18, 90 ) COMMON BLK * * Purpose * ======= * * Compute * * Y(:,1:M) = op(A)*X(:,1:M) * * where op(A) is A or A' (the transpose of A). The matrix A is the * Tolosa matrix. * * The common block BLK store the matrix data. The read in file * readb.f must be called to read in the data file block. * * Arguments * ========= * * TRANS (input) INTEGER * If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M) * If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M) * * N (input) * The order of the matrix A, N has to be 90+5*K for a negative * integer K. * * M (input) * The number of columns of X to multiply. * * X (input) DOUBLE PRECISION array, dimension (LDX,M) * contains the matrix X. * * LDX (input) INTEGER * The leading dimension of the array X, LDX >= max(1,N) * * Y (output) DOUBLE PRECISION array, dimension (LDY,M) * Y contains the product of the matrix A with X. * * LDY (input) INTEGER * The leading dimension of the array Y, LDY >= max(1,N) * * ================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO, THREE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, $ THREE = 3.0D+0 ) DOUBLE PRECISION LAMBD1 PARAMETER ( LAMBD1 = -0.12098D+02 ) DOUBLE PRECISION PARAM1, PARAM2 PARAMETER ( PARAM1 = 0.299D+0, PARAM2 = 1.0D-3 ) * .. * .. Local Scalars .. INTEGER I, K, TB DOUBLE PRECISION ALPHA, C1, C2, LAMBD2, LAMBD3, OMEGA * .. * .. Intrinsic functions .. INTRINSIC DBLE * .. * .. External Subroutines .. EXTERNAL DGEMM * .. * .. Executable Statements .. * * .. Initialize the parameters * LAMBD2 = TWO*LAMBD1 LAMBD3 = THREE*LAMBD1 * IF( TRANS.EQ.0 )THEN * * Y(:,1:M) = A*X(:,1:M) * IF ( N.EQ.90 ) THEN * DO 20 K = 1, M DO 10 I = 1,18 Y( I,K ) = X( 18+I,K ) Y( 2*18+I,K ) = X( 18+I,K ) + LAMBD1*X( 2*18+I,K ) Y( 3*18+I,K ) = X( 18+I,K ) + LAMBD2*X( 3*18+I,K ) Y( 4*18+I,K ) = X( 18+I,K ) + LAMBD3*X( 4*18+I,K ) 10 CONTINUE 20 CONTINUE * CALL DGEMM( 'N', 'N', 18, M, 90, ONE, BLK, 18, $ X, LDX, ZERO, Y( 18+1,1 ), LDY ) * ELSE IF ( N.GT.90 ) THEN * * N > 90 * TB = N/5 DO 40 K = 1, M DO 30 I = 1,TB Y( I,K ) = X( TB+I,K ) Y( 2*TB+I,K ) = X( TB+I,K ) + LAMBD1*X( 2*TB+I,K ) Y( 3*TB+I,K ) = X( TB+I,K ) + LAMBD2*X( 3*TB+I,K ) Y( 4*TB+I,K ) = X( TB+I,K ) + LAMBD3*X( 4*TB+I,K ) 30 CONTINUE 40 CONTINUE * CALL DGEMM( 'N', 'N', 18, M, 18, ONE, BLK, 18, $ X, LDX, ZERO, Y( TB+1,1 ), LDY ) * CALL DGEMM( 'N', 'N', 18, M, 18, ONE, BLK( 1,18+1 ), 18, $ X( TB+1,1 ), LDX, ONE, Y( TB+1,1 ), LDY ) * CALL DGEMM( 'N', 'N', 18, M, 18, ONE, BLK( 1,2*18+1 ), $ 18, X( 2*TB+1,1 ), LDX, ONE, Y( TB+1,1 ), $ LDY ) * CALL DGEMM( 'N', 'N', 18, M, 18, ONE, BLK( 1,3*18+1 ), $ 18, X( 3*TB+1,1 ), LDX, ONE, Y( TB+1,1 ), $ LDY ) * CALL DGEMM( 'N', 'N', 18, M, 18, ONE, BLK( 1,4*18+1 ), $ 18, X( 4*TB+1,1 ), LDX, ONE, Y( TB+1,1 ), $ LDY ) * C1 = PARAM1 / DBLE( TB-18 ) C2 = PARAM2 - C1 DO 60 K = 1, M DO 50 I = 1, TB - 18 OMEGA = DBLE( 150+6*I ) ALPHA = C1*DBLE( I ) + C2 Y( TB+18+I, K ) = -OMEGA*OMEGA*X( 18+I, K ) $ -TWO*ALPHA*OMEGA*X( TB+18+I, K ) 50 CONTINUE 60 CONTINUE * END IF * ELSE IF( TRANS.EQ.1 )THEN * * Compute Y(:,1:M) = A'*X(:,1:M) * IF ( N.EQ.90 )THEN * CALL DGEMM( 'T', 'N', 90, M, 18, ONE, BLK, 18, $ X( 18+1, 1 ), LDX, ZERO, Y, LDY ) * DO 80 K = 1, M DO 70 I = 1, 18 Y( 18+I, K ) = X( I, K ) + Y( 18+I, K ) + $ X( 2*18+I, K ) + X( 3*18+I, K ) + $ X( 4*18+I, K ) Y( 2*18+I, K ) = Y( 2*18+I, K ) + $ LAMBD1*X( 2*18+I, K ) Y( 3*18+I, K ) = Y( 3*18+I, K ) + $ LAMBD2*X( 3*18+I, K ) Y( 4*18+I, K ) = Y( 4*18+I, K ) + $ LAMBD3*X( 4*18+I, K ) 70 CONTINUE 80 CONTINUE * ELSE IF ( N.GT.90 ) THEN * * N > 90 * TB = N / 5 CALL DGEMM( 'T', 'N', 18, M, 18, ONE, BLK, 18, $ X( TB+1, 1 ), LDX, ZERO, Y, LDY ) * CALL DGEMM( 'T', 'N', 18, M, 18, ONE, BLK( 1,18+1 ), $ 18, X( TB+1, 1 ), LDX, ZERO, Y( TB+1, 1 ), $ LDY ) * CALL DGEMM( 'T', 'N', 18, M, 18, ONE, BLK( 1,2*18+1 ), $ 18, X( TB+1, 1 ), LDX, ZERO, Y( 2*TB+1, 1 ), $ LDY ) * CALL DGEMM( 'T', 'N', 18, M, 18, ONE, BLK( 1,3*18+1 ), $ 18, X( TB+1, 1 ), LDX, ZERO, Y( 3*TB+1, 1 ), $ LDY ) * CALL DGEMM( 'T', 'N', 18, M, 18, ONE, BLK( 1,4*18+1 ), $ 18, X( TB+1, 1 ), LDX, ZERO, Y( 4*TB+1, 1 ), $ LDY ) * C1 = PARAM1 / DBLE( TB - 18 ) C2 = PARAM2 - C1 DO 100 K = 1, M DO 90 I = 1, TB - 18 OMEGA = DBLE( 150 + 6*I ) ALPHA = C1*DBLE( I ) + C2 Y( 18+I, K ) = -OMEGA*OMEGA*X( TB+18+I, K ) Y( TB+18+I, K ) = -TWO*ALPHA*OMEGA*X( TB+18+I, K ) 90 CONTINUE 100 CONTINUE * DO 140 K = 1, M * DO 110 I = 1, TB Y( TB+I, K ) = X( I, K ) + Y( TB+I, K ) + $ X( 2*TB+I, K ) + X( 3*TB+I, K ) + $ X( 4*TB+I, K ) 110 CONTINUE * DO 120 I = 1, 18 Y( 2*TB+I, K ) = Y( 2*TB+I, K ) + $ LAMBD1*X( 2*TB+I, K ) Y( 3*TB+I, K ) = Y( 3*TB+I, K ) + $ LAMBD2*X( 3*TB+I, K ) Y( 4*TB+I, K ) = Y( 4*TB+I, K ) + $ LAMBD3*X( 4*TB+I, K ) 120 CONTINUE * DO 130 I = 1, TB - 18 Y( 2*TB+18+I, K ) = LAMBD1*X( 2*TB+18+I, K ) Y( 3*TB+18+I, K ) = LAMBD2*X( 3*TB+18+I, K ) Y( 4*TB+18+I, K ) = LAMBD3*X( 4*TB+18+I, K ) 130 CONTINUE * 140 CONTINUE * END IF * END IF * RETURN * * END IF MVMTLS * END * * the following routines are from BLAS * SUBROUTINE DGEMM ( 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 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 .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( 'DGEMM ', 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 DGEMM . * END SUBROUTINE XERBLA( SRNAME, INFO ) * * -- LAPACK auxiliary routine (preliminary version) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER*6 SRNAME INTEGER INFO * .. * * Purpose * ======= * * XERBLA is an error handler for the LAPACK routines. * It is called by an LAPACK routine if an input parameter has an * invalid value. A message is printed and execution stops. * * Installers may consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Arguments * ========= * * SRNAME (input) CHARACTER*6 * The name of the routine which called XERBLA. * * INFO (input) INTEGER * The position of the invalid parameter in the parameter list * of the calling routine. * * WRITE( *, FMT = 9999 )SRNAME, INFO * STOP * 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', $ 'an illegal value' ) * * End of XERBLA * END LOGICAL FUNCTION LSAME( CA, CB ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER CA, CB * .. * * Purpose * ======= * * LSAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * ===================================================================== * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA.EQ.CB IF( LSAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB * * RETURN * * End of LSAME * END