LAPACK: TESTING/MATGEN/clatm5.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLATM5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD,
00002      $                   E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA,
00003      $                   QBLCKB )
00004 *
00005 *  -- LAPACK test routine (version 3.1) --
00006 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
00011      $                   PRTYPE, QBLCKA, QBLCKB
00012       REAL               ALPHA
00013 *     ..
00014 *     .. Array Arguments ..
00015       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * ),
00016      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
00017      $                   L( LDL, * ), R( LDR, * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  CLATM5 generates matrices involved in the Generalized Sylvester
00024 *  equation:
00025 *
00026 *      A * R - L * B = C
00027 *      D * R - L * E = F
00028 *
00029 *  They also satisfy (the diagonalization condition)
00030 *
00031 *   [ I -L ] ( [ A  -C ], [ D -F ] ) [ I  R ] = ( [ A    ], [ D    ] )
00032 *   [    I ] ( [     B ]  [    E ] ) [    I ]   ( [    B ]  [    E ] )
00033 *
00034 *
00035 *  Arguments
00036 *  =========
00037 *
00038 *  PRTYPE  (input) INTEGER
00039 *          "Points" to a certian type of the matrices to generate
00040 *          (see futher details).
00041 *
00042 *  M       (input) INTEGER
00043 *          Specifies the order of A and D and the number of rows in
00044 *          C, F,  R and L.
00045 *
00046 *  N       (input) INTEGER
00047 *          Specifies the order of B and E and the number of columns in
00048 *          C, F, R and L.
00049 *
00050 *  A       (output) COMPLEX array, dimension (LDA, M).
00051 *          On exit A M-by-M is initialized according to PRTYPE.
00052 *
00053 *  LDA     (input) INTEGER
00054 *          The leading dimension of A.
00055 *
00056 *  B       (output) COMPLEX array, dimension (LDB, N).
00057 *          On exit B N-by-N is initialized according to PRTYPE.
00058 *
00059 *  LDB     (input) INTEGER
00060 *          The leading dimension of B.
00061 *
00062 *  C       (output) COMPLEX array, dimension (LDC, N).
00063 *          On exit C M-by-N is initialized according to PRTYPE.
00064 *
00065 *  LDC     (input) INTEGER
00066 *          The leading dimension of C.
00067 *
00068 *  D       (output) COMPLEX array, dimension (LDD, M).
00069 *          On exit D M-by-M is initialized according to PRTYPE.
00070 *
00071 *  LDD     (input) INTEGER
00072 *          The leading dimension of D.
00073 *
00074 *  E       (output) COMPLEX array, dimension (LDE, N).
00075 *          On exit E N-by-N is initialized according to PRTYPE.
00076 *
00077 *  LDE     (input) INTEGER
00078 *          The leading dimension of E.
00079 *
00080 *  F       (output) COMPLEX array, dimension (LDF, N).
00081 *          On exit F M-by-N is initialized according to PRTYPE.
00082 *
00083 *  LDF     (input) INTEGER
00084 *          The leading dimension of F.
00085 *
00086 *  R       (output) COMPLEX array, dimension (LDR, N).
00087 *          On exit R M-by-N is initialized according to PRTYPE.
00088 *
00089 *  LDR     (input) INTEGER
00090 *          The leading dimension of R.
00091 *
00092 *  L       (output) COMPLEX array, dimension (LDL, N).
00093 *          On exit L M-by-N is initialized according to PRTYPE.
00094 *
00095 *  LDL     (input) INTEGER
00096 *          The leading dimension of L.
00097 *
00098 *  ALPHA   (input) REAL
00099 *          Parameter used in generating PRTYPE = 1 and 5 matrices.
00100 *
00101 *  QBLCKA  (input) INTEGER
00102 *          When PRTYPE = 3, specifies the distance between 2-by-2
00103 *          blocks on the diagonal in A. Otherwise, QBLCKA is not
00104 *          referenced. QBLCKA > 1.
00105 *
00106 *  QBLCKB  (input) INTEGER
00107 *          When PRTYPE = 3, specifies the distance between 2-by-2
00108 *          blocks on the diagonal in B. Otherwise, QBLCKB is not
00109 *          referenced. QBLCKB > 1.
00110 *
00111 *
00112 *  Further Details
00113 *  ===============
00114 *
00115 *  PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices
00116 *
00117 *             A : if (i == j) then A(i, j) = 1.0
00118 *                 if (j == i + 1) then A(i, j) = -1.0
00119 *                 else A(i, j) = 0.0,            i, j = 1...M
00120 *
00121 *             B : if (i == j) then B(i, j) = 1.0 - ALPHA
00122 *                 if (j == i + 1) then B(i, j) = 1.0
00123 *                 else B(i, j) = 0.0,            i, j = 1...N
00124 *
00125 *             D : if (i == j) then D(i, j) = 1.0
00126 *                 else D(i, j) = 0.0,            i, j = 1...M
00127 *
00128 *             E : if (i == j) then E(i, j) = 1.0
00129 *                 else E(i, j) = 0.0,            i, j = 1...N
00130 *
00131 *             L =  R are chosen from [-10...10],
00132 *                  which specifies the right hand sides (C, F).
00133 *
00134 *  PRTYPE = 2 or 3: Triangular and/or quasi- triangular.
00135 *
00136 *             A : if (i <= j) then A(i, j) = [-1...1]
00137 *                 else A(i, j) = 0.0,             i, j = 1...M
00138 *
00139 *                 if (PRTYPE = 3) then
00140 *                    A(k + 1, k + 1) = A(k, k)
00141 *                    A(k + 1, k) = [-1...1]
00142 *                    sign(A(k, k + 1) = -(sin(A(k + 1, k))
00143 *                        k = 1, M - 1, QBLCKA
00144 *
00145 *             B : if (i <= j) then B(i, j) = [-1...1]
00146 *                 else B(i, j) = 0.0,            i, j = 1...N
00147 *
00148 *                 if (PRTYPE = 3) then
00149 *                    B(k + 1, k + 1) = B(k, k)
00150 *                    B(k + 1, k) = [-1...1]
00151 *                    sign(B(k, k + 1) = -(sign(B(k + 1, k))
00152 *                        k = 1, N - 1, QBLCKB
00153 *
00154 *             D : if (i <= j) then D(i, j) = [-1...1].
00155 *                 else D(i, j) = 0.0,            i, j = 1...M
00156 *
00157 *
00158 *             E : if (i <= j) then D(i, j) = [-1...1]
00159 *                 else E(i, j) = 0.0,            i, j = 1...N
00160 *
00161 *                 L, R are chosen from [-10...10],
00162 *                 which specifies the right hand sides (C, F).
00163 *
00164 *  PRTYPE = 4 Full
00165 *             A(i, j) = [-10...10]
00166 *             D(i, j) = [-1...1]    i,j = 1...M
00167 *             B(i, j) = [-10...10]
00168 *             E(i, j) = [-1...1]    i,j = 1...N
00169 *             R(i, j) = [-10...10]
00170 *             L(i, j) = [-1...1]    i = 1..M ,j = 1...N
00171 *
00172 *             L, R specifies the right hand sides (C, F).
00173 *
00174 *  PRTYPE = 5 special case common and/or close eigs.
00175 *
00176 *  =====================================================================
00177 *
00178 *     .. Parameters ..
00179       COMPLEX            ONE, TWO, ZERO, HALF, TWENTY
00180       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
00181      $                   TWO = ( 2.0E+0, 0.0E+0 ),
00182      $                   ZERO = ( 0.0E+0, 0.0E+0 ),
00183      $                   HALF = ( 0.5E+0, 0.0E+0 ),
00184      $                   TWENTY = ( 2.0E+1, 0.0E+0 ) )
00185 *     ..
00186 *     .. Local Scalars ..
00187       INTEGER            I, J, K
00188       COMPLEX            IMEPS, REEPS
00189 *     ..
00190 *     .. Intrinsic Functions ..
00191       INTRINSIC          CMPLX, MOD, SIN
00192 *     ..
00193 *     .. External Subroutines ..
00194       EXTERNAL           CGEMM
00195 *     ..
00196 *     .. Executable Statements ..
00197 *
00198       IF( PRTYPE.EQ.1 ) THEN
00199          DO 20 I = 1, M
00200             DO 10 J = 1, M
00201                IF( I.EQ.J ) THEN
00202                   A( I, J ) = ONE
00203                   D( I, J ) = ONE
00204                ELSE IF( I.EQ.J-1 ) THEN
00205                   A( I, J ) = -ONE
00206                   D( I, J ) = ZERO
00207                ELSE
00208                   A( I, J ) = ZERO
00209                   D( I, J ) = ZERO
00210                END IF
00211    10       CONTINUE
00212    20    CONTINUE
00213 *
00214          DO 40 I = 1, N
00215             DO 30 J = 1, N
00216                IF( I.EQ.J ) THEN
00217                   B( I, J ) = ONE - ALPHA
00218                   E( I, J ) = ONE
00219                ELSE IF( I.EQ.J-1 ) THEN
00220                   B( I, J ) = ONE
00221                   E( I, J ) = ZERO
00222                ELSE
00223                   B( I, J ) = ZERO
00224                   E( I, J ) = ZERO
00225                END IF
00226    30       CONTINUE
00227    40    CONTINUE
00228 *
00229          DO 60 I = 1, M
00230             DO 50 J = 1, N
00231                R( I, J ) = ( HALF-SIN( CMPLX( I / J ) ) )*TWENTY
00232                L( I, J ) = R( I, J )
00233    50       CONTINUE
00234    60    CONTINUE
00235 *
00236       ELSE IF( PRTYPE.EQ.2 .OR. PRTYPE.EQ.3 ) THEN
00237          DO 80 I = 1, M
00238             DO 70 J = 1, M
00239                IF( I.LE.J ) THEN
00240                   A( I, J ) = ( HALF-SIN( CMPLX( I ) ) )*TWO
00241                   D( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*TWO
00242                ELSE
00243                   A( I, J ) = ZERO
00244                   D( I, J ) = ZERO
00245                END IF
00246    70       CONTINUE
00247    80    CONTINUE
00248 *
00249          DO 100 I = 1, N
00250             DO 90 J = 1, N
00251                IF( I.LE.J ) THEN
00252                   B( I, J ) = ( HALF-SIN( CMPLX( I+J ) ) )*TWO
00253                   E( I, J ) = ( HALF-SIN( CMPLX( J ) ) )*TWO
00254                ELSE
00255                   B( I, J ) = ZERO
00256                   E( I, J ) = ZERO
00257                END IF
00258    90       CONTINUE
00259   100    CONTINUE
00260 *
00261          DO 120 I = 1, M
00262             DO 110 J = 1, N
00263                R( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*TWENTY
00264                L( I, J ) = ( HALF-SIN( CMPLX( I+J ) ) )*TWENTY
00265   110       CONTINUE
00266   120    CONTINUE
00267 *
00268          IF( PRTYPE.EQ.3 ) THEN
00269             IF( QBLCKA.LE.1 )
00270      $         QBLCKA = 2
00271             DO 130 K = 1, M - 1, QBLCKA
00272                A( K+1, K+1 ) = A( K, K )
00273                A( K+1, K ) = -SIN( A( K, K+1 ) )
00274   130       CONTINUE
00275 *
00276             IF( QBLCKB.LE.1 )
00277      $         QBLCKB = 2
00278             DO 140 K = 1, N - 1, QBLCKB
00279                B( K+1, K+1 ) = B( K, K )
00280                B( K+1, K ) = -SIN( B( K, K+1 ) )
00281   140       CONTINUE
00282          END IF
00283 *
00284       ELSE IF( PRTYPE.EQ.4 ) THEN
00285          DO 160 I = 1, M
00286             DO 150 J = 1, M
00287                A( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*TWENTY
00288                D( I, J ) = ( HALF-SIN( CMPLX( I+J ) ) )*TWO
00289   150       CONTINUE
00290   160    CONTINUE
00291 *
00292          DO 180 I = 1, N
00293             DO 170 J = 1, N
00294                B( I, J ) = ( HALF-SIN( CMPLX( I+J ) ) )*TWENTY
00295                E( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*TWO
00296   170       CONTINUE
00297   180    CONTINUE
00298 *
00299          DO 200 I = 1, M
00300             DO 190 J = 1, N
00301                R( I, J ) = ( HALF-SIN( CMPLX( J / I ) ) )*TWENTY
00302                L( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*TWO
00303   190       CONTINUE
00304   200    CONTINUE
00305 *
00306       ELSE IF( PRTYPE.GE.5 ) THEN
00307          REEPS = HALF*TWO*TWENTY / ALPHA
00308          IMEPS = ( HALF-TWO ) / ALPHA
00309          DO 220 I = 1, M
00310             DO 210 J = 1, N
00311                R( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*ALPHA / TWENTY
00312                L( I, J ) = ( HALF-SIN( CMPLX( I+J ) ) )*ALPHA / TWENTY
00313   210       CONTINUE
00314   220    CONTINUE
00315 *
00316          DO 230 I = 1, M
00317             D( I, I ) = ONE
00318   230    CONTINUE
00319 *
00320          DO 240 I = 1, M
00321             IF( I.LE.4 ) THEN
00322                A( I, I ) = ONE
00323                IF( I.GT.2 )
00324      $            A( I, I ) = ONE + REEPS
00325                IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
00326                   A( I, I+1 ) = IMEPS
00327                ELSE IF( I.GT.1 ) THEN
00328                   A( I, I-1 ) = -IMEPS
00329                END IF
00330             ELSE IF( I.LE.8 ) THEN
00331                IF( I.LE.6 ) THEN
00332                   A( I, I ) = REEPS
00333                ELSE
00334                   A( I, I ) = -REEPS
00335                END IF
00336                IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
00337                   A( I, I+1 ) = ONE
00338                ELSE IF( I.GT.1 ) THEN
00339                   A( I, I-1 ) = -ONE
00340                END IF
00341             ELSE
00342                A( I, I ) = ONE
00343                IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
00344                   A( I, I+1 ) = IMEPS*2
00345                ELSE IF( I.GT.1 ) THEN
00346                   A( I, I-1 ) = -IMEPS*2
00347                END IF
00348             END IF
00349   240    CONTINUE
00350 *
00351          DO 250 I = 1, N
00352             E( I, I ) = ONE
00353             IF( I.LE.4 ) THEN
00354                B( I, I ) = -ONE
00355                IF( I.GT.2 )
00356      $            B( I, I ) = ONE - REEPS
00357                IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
00358                   B( I, I+1 ) = IMEPS
00359                ELSE IF( I.GT.1 ) THEN
00360                   B( I, I-1 ) = -IMEPS
00361                END IF
00362             ELSE IF( I.LE.8 ) THEN
00363                IF( I.LE.6 ) THEN
00364                   B( I, I ) = REEPS
00365                ELSE
00366                   B( I, I ) = -REEPS
00367                END IF
00368                IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
00369                   B( I, I+1 ) = ONE + IMEPS
00370                ELSE IF( I.GT.1 ) THEN
00371                   B( I, I-1 ) = -ONE - IMEPS
00372                END IF
00373             ELSE
00374                B( I, I ) = ONE - REEPS
00375                IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
00376                   B( I, I+1 ) = IMEPS*2
00377                ELSE IF( I.GT.1 ) THEN
00378                   B( I, I-1 ) = -IMEPS*2
00379                END IF
00380             END IF
00381   250    CONTINUE
00382       END IF
00383 *
00384 *     Compute rhs (C, F)
00385 *
00386       CALL CGEMM( 'N', 'N', M, N, M, ONE, A, LDA, R, LDR, ZERO, C, LDC )
00387       CALL CGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, B, LDB, ONE, C, LDC )
00388       CALL CGEMM( 'N', 'N', M, N, M, ONE, D, LDD, R, LDR, ZERO, F, LDF )
00389       CALL CGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, E, LDE, ONE, F, LDF )
00390 *
00391 *     End of CLATM5
00392 *
00393       END