SUBROUTINE DLATM5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, \$ E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, \$ QBLCKB ) * * -- LAPACK test routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N, \$ PRTYPE, QBLCKA, QBLCKB DOUBLE PRECISION ALPHA * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ), \$ D( LDD, * ), E( LDE, * ), F( LDF, * ), \$ L( LDL, * ), R( LDR, * ) * .. * * Purpose * ======= * * DLATM5 generates matrices involved in the Generalized Sylvester * equation: * * A * R - L * B = C * D * R - L * E = F * * They also satisfy (the diagonalization condition) * * [ I -L ] ( [ A -C ], [ D -F ] ) [ I R ] = ( [ A ], [ D ] ) * [ I ] ( [ B ] [ E ] ) [ I ] ( [ B ] [ E ] ) * * * Arguments * ========= * * PRTYPE (input) INTEGER * "Points" to a certian type of the matrices to generate * (see futher details). * * M (input) INTEGER * Specifies the order of A and D and the number of rows in * C, F, R and L. * * N (input) INTEGER * Specifies the order of B and E and the number of columns in * C, F, R and L. * * A (output) DOUBLE PRECISION array, dimension (LDA, M). * On exit A M-by-M is initialized according to PRTYPE. * * LDA (input) INTEGER * The leading dimension of A. * * B (output) DOUBLE PRECISION array, dimension (LDB, N). * On exit B N-by-N is initialized according to PRTYPE. * * LDB (input) INTEGER * The leading dimension of B. * * C (output) DOUBLE PRECISION array, dimension (LDC, N). * On exit C M-by-N is initialized according to PRTYPE. * * LDC (input) INTEGER * The leading dimension of C. * * D (output) DOUBLE PRECISION array, dimension (LDD, M). * On exit D M-by-M is initialized according to PRTYPE. * * LDD (input) INTEGER * The leading dimension of D. * * E (output) DOUBLE PRECISION array, dimension (LDE, N). * On exit E N-by-N is initialized according to PRTYPE. * * LDE (input) INTEGER * The leading dimension of E. * * F (output) DOUBLE PRECISION array, dimension (LDF, N). * On exit F M-by-N is initialized according to PRTYPE. * * LDF (input) INTEGER * The leading dimension of F. * * R (output) DOUBLE PRECISION array, dimension (LDR, N). * On exit R M-by-N is initialized according to PRTYPE. * * LDR (input) INTEGER * The leading dimension of R. * * L (output) DOUBLE PRECISION array, dimension (LDL, N). * On exit L M-by-N is initialized according to PRTYPE. * * LDL (input) INTEGER * The leading dimension of L. * * ALPHA (input) DOUBLE PRECISION * Parameter used in generating PRTYPE = 1 and 5 matrices. * * QBLCKA (input) INTEGER * When PRTYPE = 3, specifies the distance between 2-by-2 * blocks on the diagonal in A. Otherwise, QBLCKA is not * referenced. QBLCKA > 1. * * QBLCKB (input) INTEGER * When PRTYPE = 3, specifies the distance between 2-by-2 * blocks on the diagonal in B. Otherwise, QBLCKB is not * referenced. QBLCKB > 1. * * * Further Details * =============== * * PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices * * A : if (i == j) then A(i, j) = 1.0 * if (j == i + 1) then A(i, j) = -1.0 * else A(i, j) = 0.0, i, j = 1...M * * B : if (i == j) then B(i, j) = 1.0 - ALPHA * if (j == i + 1) then B(i, j) = 1.0 * else B(i, j) = 0.0, i, j = 1...N * * D : if (i == j) then D(i, j) = 1.0 * else D(i, j) = 0.0, i, j = 1...M * * E : if (i == j) then E(i, j) = 1.0 * else E(i, j) = 0.0, i, j = 1...N * * L = R are chosen from [-10...10], * which specifies the right hand sides (C, F). * * PRTYPE = 2 or 3: Triangular and/or quasi- triangular. * * A : if (i <= j) then A(i, j) = [-1...1] * else A(i, j) = 0.0, i, j = 1...M * * if (PRTYPE = 3) then * A(k + 1, k + 1) = A(k, k) * A(k + 1, k) = [-1...1] * sign(A(k, k + 1) = -(sin(A(k + 1, k)) * k = 1, M - 1, QBLCKA * * B : if (i <= j) then B(i, j) = [-1...1] * else B(i, j) = 0.0, i, j = 1...N * * if (PRTYPE = 3) then * B(k + 1, k + 1) = B(k, k) * B(k + 1, k) = [-1...1] * sign(B(k, k + 1) = -(sign(B(k + 1, k)) * k = 1, N - 1, QBLCKB * * D : if (i <= j) then D(i, j) = [-1...1]. * else D(i, j) = 0.0, i, j = 1...M * * * E : if (i <= j) then D(i, j) = [-1...1] * else E(i, j) = 0.0, i, j = 1...N * * L, R are chosen from [-10...10], * which specifies the right hand sides (C, F). * * PRTYPE = 4 Full * A(i, j) = [-10...10] * D(i, j) = [-1...1] i,j = 1...M * B(i, j) = [-10...10] * E(i, j) = [-1...1] i,j = 1...N * R(i, j) = [-10...10] * L(i, j) = [-1...1] i = 1..M ,j = 1...N * * L, R specifies the right hand sides (C, F). * * PRTYPE = 5 special case common and/or close eigs. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO, TWENTY, HALF, TWO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0, TWENTY = 2.0D+1, \$ HALF = 0.5D+0, TWO = 2.0D+0 ) * .. * .. Local Scalars .. INTEGER I, J, K DOUBLE PRECISION IMEPS, REEPS * .. * .. Intrinsic Functions .. INTRINSIC DBLE, MOD, SIN * .. * .. External Subroutines .. EXTERNAL DGEMM * .. * .. Executable Statements .. * IF( PRTYPE.EQ.1 ) THEN DO 20 I = 1, M DO 10 J = 1, M IF( I.EQ.J ) THEN A( I, J ) = ONE D( I, J ) = ONE ELSE IF( I.EQ.J-1 ) THEN A( I, J ) = -ONE D( I, J ) = ZERO ELSE A( I, J ) = ZERO D( I, J ) = ZERO END IF 10 CONTINUE 20 CONTINUE * DO 40 I = 1, N DO 30 J = 1, N IF( I.EQ.J ) THEN B( I, J ) = ONE - ALPHA E( I, J ) = ONE ELSE IF( I.EQ.J-1 ) THEN B( I, J ) = ONE E( I, J ) = ZERO ELSE B( I, J ) = ZERO E( I, J ) = ZERO END IF 30 CONTINUE 40 CONTINUE * DO 60 I = 1, M DO 50 J = 1, N R( I, J ) = ( HALF-SIN( DBLE( I / J ) ) )*TWENTY L( I, J ) = R( I, J ) 50 CONTINUE 60 CONTINUE * ELSE IF( PRTYPE.EQ.2 .OR. PRTYPE.EQ.3 ) THEN DO 80 I = 1, M DO 70 J = 1, M IF( I.LE.J ) THEN A( I, J ) = ( HALF-SIN( DBLE( I ) ) )*TWO D( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*TWO ELSE A( I, J ) = ZERO D( I, J ) = ZERO END IF 70 CONTINUE 80 CONTINUE * DO 100 I = 1, N DO 90 J = 1, N IF( I.LE.J ) THEN B( I, J ) = ( HALF-SIN( DBLE( I+J ) ) )*TWO E( I, J ) = ( HALF-SIN( DBLE( J ) ) )*TWO ELSE B( I, J ) = ZERO E( I, J ) = ZERO END IF 90 CONTINUE 100 CONTINUE * DO 120 I = 1, M DO 110 J = 1, N R( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*TWENTY L( I, J ) = ( HALF-SIN( DBLE( I+J ) ) )*TWENTY 110 CONTINUE 120 CONTINUE * IF( PRTYPE.EQ.3 ) THEN IF( QBLCKA.LE.1 ) \$ QBLCKA = 2 DO 130 K = 1, M - 1, QBLCKA A( K+1, K+1 ) = A( K, K ) A( K+1, K ) = -SIN( A( K, K+1 ) ) 130 CONTINUE * IF( QBLCKB.LE.1 ) \$ QBLCKB = 2 DO 140 K = 1, N - 1, QBLCKB B( K+1, K+1 ) = B( K, K ) B( K+1, K ) = -SIN( B( K, K+1 ) ) 140 CONTINUE END IF * ELSE IF( PRTYPE.EQ.4 ) THEN DO 160 I = 1, M DO 150 J = 1, M A( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*TWENTY D( I, J ) = ( HALF-SIN( DBLE( I+J ) ) )*TWO 150 CONTINUE 160 CONTINUE * DO 180 I = 1, N DO 170 J = 1, N B( I, J ) = ( HALF-SIN( DBLE( I+J ) ) )*TWENTY E( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*TWO 170 CONTINUE 180 CONTINUE * DO 200 I = 1, M DO 190 J = 1, N R( I, J ) = ( HALF-SIN( DBLE( J / I ) ) )*TWENTY L( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*TWO 190 CONTINUE 200 CONTINUE * ELSE IF( PRTYPE.GE.5 ) THEN REEPS = HALF*TWO*TWENTY / ALPHA IMEPS = ( HALF-TWO ) / ALPHA DO 220 I = 1, M DO 210 J = 1, N R( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*ALPHA / TWENTY L( I, J ) = ( HALF-SIN( DBLE( I+J ) ) )*ALPHA / TWENTY 210 CONTINUE 220 CONTINUE * DO 230 I = 1, M D( I, I ) = ONE 230 CONTINUE * DO 240 I = 1, M IF( I.LE.4 ) THEN A( I, I ) = ONE IF( I.GT.2 ) \$ A( I, I ) = ONE + REEPS IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN A( I, I+1 ) = IMEPS ELSE IF( I.GT.1 ) THEN A( I, I-1 ) = -IMEPS END IF ELSE IF( I.LE.8 ) THEN IF( I.LE.6 ) THEN A( I, I ) = REEPS ELSE A( I, I ) = -REEPS END IF IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN A( I, I+1 ) = ONE ELSE IF( I.GT.1 ) THEN A( I, I-1 ) = -ONE END IF ELSE A( I, I ) = ONE IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN A( I, I+1 ) = IMEPS*2 ELSE IF( I.GT.1 ) THEN A( I, I-1 ) = -IMEPS*2 END IF END IF 240 CONTINUE * DO 250 I = 1, N E( I, I ) = ONE IF( I.LE.4 ) THEN B( I, I ) = -ONE IF( I.GT.2 ) \$ B( I, I ) = ONE - REEPS IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN B( I, I+1 ) = IMEPS ELSE IF( I.GT.1 ) THEN B( I, I-1 ) = -IMEPS END IF ELSE IF( I.LE.8 ) THEN IF( I.LE.6 ) THEN B( I, I ) = REEPS ELSE B( I, I ) = -REEPS END IF IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN B( I, I+1 ) = ONE + IMEPS ELSE IF( I.GT.1 ) THEN B( I, I-1 ) = -ONE - IMEPS END IF ELSE B( I, I ) = ONE - REEPS IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN B( I, I+1 ) = IMEPS*2 ELSE IF( I.GT.1 ) THEN B( I, I-1 ) = -IMEPS*2 END IF END IF 250 CONTINUE END IF * * Compute rhs (C, F) * CALL DGEMM( 'N', 'N', M, N, M, ONE, A, LDA, R, LDR, ZERO, C, LDC ) CALL DGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, B, LDB, ONE, C, LDC ) CALL DGEMM( 'N', 'N', M, N, M, ONE, D, LDD, R, LDR, ZERO, F, LDF ) CALL DGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, E, LDE, ONE, F, LDF ) * * End of DLATM5 * END