LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
clatm5.f
Go to the documentation of this file.
00001 *> \brief \b CLATM5
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE CLATM5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD,
00012 *                          E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA,
00013 *                          QBLCKB )
00014 * 
00015 *       .. Scalar Arguments ..
00016 *       INTEGER            LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
00017 *      $                   PRTYPE, QBLCKA, QBLCKB
00018 *       REAL               ALPHA
00019 *       ..
00020 *       .. Array Arguments ..
00021 *       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * ),
00022 *      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
00023 *      $                   L( LDL, * ), R( LDR, * )
00024 *       ..
00025 *  
00026 *
00027 *> \par Purpose:
00028 *  =============
00029 *>
00030 *> \verbatim
00031 *>
00032 *> CLATM5 generates matrices involved in the Generalized Sylvester
00033 *> equation:
00034 *>
00035 *>     A * R - L * B = C
00036 *>     D * R - L * E = F
00037 *>
00038 *> They also satisfy (the diagonalization condition)
00039 *>
00040 *>  [ I -L ] ( [ A  -C ], [ D -F ] ) [ I  R ] = ( [ A    ], [ D    ] )
00041 *>  [    I ] ( [     B ]  [    E ] ) [    I ]   ( [    B ]  [    E ] )
00042 *>
00043 *> \endverbatim
00044 *
00045 *  Arguments:
00046 *  ==========
00047 *
00048 *> \param[in] PRTYPE
00049 *> \verbatim
00050 *>          PRTYPE is INTEGER
00051 *>          "Points" to a certian type of the matrices to generate
00052 *>          (see futher details).
00053 *> \endverbatim
00054 *>
00055 *> \param[in] M
00056 *> \verbatim
00057 *>          M is INTEGER
00058 *>          Specifies the order of A and D and the number of rows in
00059 *>          C, F,  R and L.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] N
00063 *> \verbatim
00064 *>          N is INTEGER
00065 *>          Specifies the order of B and E and the number of columns in
00066 *>          C, F, R and L.
00067 *> \endverbatim
00068 *>
00069 *> \param[out] A
00070 *> \verbatim
00071 *>          A is COMPLEX array, dimension (LDA, M).
00072 *>          On exit A M-by-M is initialized according to PRTYPE.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] LDA
00076 *> \verbatim
00077 *>          LDA is INTEGER
00078 *>          The leading dimension of A.
00079 *> \endverbatim
00080 *>
00081 *> \param[out] B
00082 *> \verbatim
00083 *>          B is COMPLEX array, dimension (LDB, N).
00084 *>          On exit B N-by-N is initialized according to PRTYPE.
00085 *> \endverbatim
00086 *>
00087 *> \param[in] LDB
00088 *> \verbatim
00089 *>          LDB is INTEGER
00090 *>          The leading dimension of B.
00091 *> \endverbatim
00092 *>
00093 *> \param[out] C
00094 *> \verbatim
00095 *>          C is COMPLEX array, dimension (LDC, N).
00096 *>          On exit C M-by-N is initialized according to PRTYPE.
00097 *> \endverbatim
00098 *>
00099 *> \param[in] LDC
00100 *> \verbatim
00101 *>          LDC is INTEGER
00102 *>          The leading dimension of C.
00103 *> \endverbatim
00104 *>
00105 *> \param[out] D
00106 *> \verbatim
00107 *>          D is COMPLEX array, dimension (LDD, M).
00108 *>          On exit D M-by-M is initialized according to PRTYPE.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] LDD
00112 *> \verbatim
00113 *>          LDD is INTEGER
00114 *>          The leading dimension of D.
00115 *> \endverbatim
00116 *>
00117 *> \param[out] E
00118 *> \verbatim
00119 *>          E is COMPLEX array, dimension (LDE, N).
00120 *>          On exit E N-by-N is initialized according to PRTYPE.
00121 *> \endverbatim
00122 *>
00123 *> \param[in] LDE
00124 *> \verbatim
00125 *>          LDE is INTEGER
00126 *>          The leading dimension of E.
00127 *> \endverbatim
00128 *>
00129 *> \param[out] F
00130 *> \verbatim
00131 *>          F is COMPLEX array, dimension (LDF, N).
00132 *>          On exit F M-by-N is initialized according to PRTYPE.
00133 *> \endverbatim
00134 *>
00135 *> \param[in] LDF
00136 *> \verbatim
00137 *>          LDF is INTEGER
00138 *>          The leading dimension of F.
00139 *> \endverbatim
00140 *>
00141 *> \param[out] R
00142 *> \verbatim
00143 *>          R is COMPLEX array, dimension (LDR, N).
00144 *>          On exit R M-by-N is initialized according to PRTYPE.
00145 *> \endverbatim
00146 *>
00147 *> \param[in] LDR
00148 *> \verbatim
00149 *>          LDR is INTEGER
00150 *>          The leading dimension of R.
00151 *> \endverbatim
00152 *>
00153 *> \param[out] L
00154 *> \verbatim
00155 *>          L is COMPLEX array, dimension (LDL, N).
00156 *>          On exit L M-by-N is initialized according to PRTYPE.
00157 *> \endverbatim
00158 *>
00159 *> \param[in] LDL
00160 *> \verbatim
00161 *>          LDL is INTEGER
00162 *>          The leading dimension of L.
00163 *> \endverbatim
00164 *>
00165 *> \param[in] ALPHA
00166 *> \verbatim
00167 *>          ALPHA is REAL
00168 *>          Parameter used in generating PRTYPE = 1 and 5 matrices.
00169 *> \endverbatim
00170 *>
00171 *> \param[in] QBLCKA
00172 *> \verbatim
00173 *>          QBLCKA is INTEGER
00174 *>          When PRTYPE = 3, specifies the distance between 2-by-2
00175 *>          blocks on the diagonal in A. Otherwise, QBLCKA is not
00176 *>          referenced. QBLCKA > 1.
00177 *> \endverbatim
00178 *>
00179 *> \param[in] QBLCKB
00180 *> \verbatim
00181 *>          QBLCKB is INTEGER
00182 *>          When PRTYPE = 3, specifies the distance between 2-by-2
00183 *>          blocks on the diagonal in B. Otherwise, QBLCKB is not
00184 *>          referenced. QBLCKB > 1.
00185 *> \endverbatim
00186 *
00187 *  Authors:
00188 *  ========
00189 *
00190 *> \author Univ. of Tennessee 
00191 *> \author Univ. of California Berkeley 
00192 *> \author Univ. of Colorado Denver 
00193 *> \author NAG Ltd. 
00194 *
00195 *> \date November 2011
00196 *
00197 *> \ingroup complex_matgen
00198 *
00199 *> \par Further Details:
00200 *  =====================
00201 *>
00202 *> \verbatim
00203 *>
00204 *>  PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices
00205 *>
00206 *>             A : if (i == j) then A(i, j) = 1.0
00207 *>                 if (j == i + 1) then A(i, j) = -1.0
00208 *>                 else A(i, j) = 0.0,            i, j = 1...M
00209 *>
00210 *>             B : if (i == j) then B(i, j) = 1.0 - ALPHA
00211 *>                 if (j == i + 1) then B(i, j) = 1.0
00212 *>                 else B(i, j) = 0.0,            i, j = 1...N
00213 *>
00214 *>             D : if (i == j) then D(i, j) = 1.0
00215 *>                 else D(i, j) = 0.0,            i, j = 1...M
00216 *>
00217 *>             E : if (i == j) then E(i, j) = 1.0
00218 *>                 else E(i, j) = 0.0,            i, j = 1...N
00219 *>
00220 *>             L =  R are chosen from [-10...10],
00221 *>                  which specifies the right hand sides (C, F).
00222 *>
00223 *>  PRTYPE = 2 or 3: Triangular and/or quasi- triangular.
00224 *>
00225 *>             A : if (i <= j) then A(i, j) = [-1...1]
00226 *>                 else A(i, j) = 0.0,             i, j = 1...M
00227 *>
00228 *>                 if (PRTYPE = 3) then
00229 *>                    A(k + 1, k + 1) = A(k, k)
00230 *>                    A(k + 1, k) = [-1...1]
00231 *>                    sign(A(k, k + 1) = -(sin(A(k + 1, k))
00232 *>                        k = 1, M - 1, QBLCKA
00233 *>
00234 *>             B : if (i <= j) then B(i, j) = [-1...1]
00235 *>                 else B(i, j) = 0.0,            i, j = 1...N
00236 *>
00237 *>                 if (PRTYPE = 3) then
00238 *>                    B(k + 1, k + 1) = B(k, k)
00239 *>                    B(k + 1, k) = [-1...1]
00240 *>                    sign(B(k, k + 1) = -(sign(B(k + 1, k))
00241 *>                        k = 1, N - 1, QBLCKB
00242 *>
00243 *>             D : if (i <= j) then D(i, j) = [-1...1].
00244 *>                 else D(i, j) = 0.0,            i, j = 1...M
00245 *>
00246 *>
00247 *>             E : if (i <= j) then D(i, j) = [-1...1]
00248 *>                 else E(i, j) = 0.0,            i, j = 1...N
00249 *>
00250 *>                 L, R are chosen from [-10...10],
00251 *>                 which specifies the right hand sides (C, F).
00252 *>
00253 *>  PRTYPE = 4 Full
00254 *>             A(i, j) = [-10...10]
00255 *>             D(i, j) = [-1...1]    i,j = 1...M
00256 *>             B(i, j) = [-10...10]
00257 *>             E(i, j) = [-1...1]    i,j = 1...N
00258 *>             R(i, j) = [-10...10]
00259 *>             L(i, j) = [-1...1]    i = 1..M ,j = 1...N
00260 *>
00261 *>             L, R specifies the right hand sides (C, F).
00262 *>
00263 *>  PRTYPE = 5 special case common and/or close eigs.
00264 *> \endverbatim
00265 *>
00266 *  =====================================================================
00267       SUBROUTINE CLATM5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD,
00268      $                   E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA,
00269      $                   QBLCKB )
00270 *
00271 *  -- LAPACK computational routine (version 3.4.0) --
00272 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00273 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00274 *     November 2011
00275 *
00276 *     .. Scalar Arguments ..
00277       INTEGER            LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N,
00278      $                   PRTYPE, QBLCKA, QBLCKB
00279       REAL               ALPHA
00280 *     ..
00281 *     .. Array Arguments ..
00282       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * ),
00283      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
00284      $                   L( LDL, * ), R( LDR, * )
00285 *     ..
00286 *
00287 *  =====================================================================
00288 *
00289 *     .. Parameters ..
00290       COMPLEX            ONE, TWO, ZERO, HALF, TWENTY
00291       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
00292      $                   TWO = ( 2.0E+0, 0.0E+0 ),
00293      $                   ZERO = ( 0.0E+0, 0.0E+0 ),
00294      $                   HALF = ( 0.5E+0, 0.0E+0 ),
00295      $                   TWENTY = ( 2.0E+1, 0.0E+0 ) )
00296 *     ..
00297 *     .. Local Scalars ..
00298       INTEGER            I, J, K
00299       COMPLEX            IMEPS, REEPS
00300 *     ..
00301 *     .. Intrinsic Functions ..
00302       INTRINSIC          CMPLX, MOD, SIN
00303 *     ..
00304 *     .. External Subroutines ..
00305       EXTERNAL           CGEMM
00306 *     ..
00307 *     .. Executable Statements ..
00308 *
00309       IF( PRTYPE.EQ.1 ) THEN
00310          DO 20 I = 1, M
00311             DO 10 J = 1, M
00312                IF( I.EQ.J ) THEN
00313                   A( I, J ) = ONE
00314                   D( I, J ) = ONE
00315                ELSE IF( I.EQ.J-1 ) THEN
00316                   A( I, J ) = -ONE
00317                   D( I, J ) = ZERO
00318                ELSE
00319                   A( I, J ) = ZERO
00320                   D( I, J ) = ZERO
00321                END IF
00322    10       CONTINUE
00323    20    CONTINUE
00324 *
00325          DO 40 I = 1, N
00326             DO 30 J = 1, N
00327                IF( I.EQ.J ) THEN
00328                   B( I, J ) = ONE - ALPHA
00329                   E( I, J ) = ONE
00330                ELSE IF( I.EQ.J-1 ) THEN
00331                   B( I, J ) = ONE
00332                   E( I, J ) = ZERO
00333                ELSE
00334                   B( I, J ) = ZERO
00335                   E( I, J ) = ZERO
00336                END IF
00337    30       CONTINUE
00338    40    CONTINUE
00339 *
00340          DO 60 I = 1, M
00341             DO 50 J = 1, N
00342                R( I, J ) = ( HALF-SIN( CMPLX( I / J ) ) )*TWENTY
00343                L( I, J ) = R( I, J )
00344    50       CONTINUE
00345    60    CONTINUE
00346 *
00347       ELSE IF( PRTYPE.EQ.2 .OR. PRTYPE.EQ.3 ) THEN
00348          DO 80 I = 1, M
00349             DO 70 J = 1, M
00350                IF( I.LE.J ) THEN
00351                   A( I, J ) = ( HALF-SIN( CMPLX( I ) ) )*TWO
00352                   D( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*TWO
00353                ELSE
00354                   A( I, J ) = ZERO
00355                   D( I, J ) = ZERO
00356                END IF
00357    70       CONTINUE
00358    80    CONTINUE
00359 *
00360          DO 100 I = 1, N
00361             DO 90 J = 1, N
00362                IF( I.LE.J ) THEN
00363                   B( I, J ) = ( HALF-SIN( CMPLX( I+J ) ) )*TWO
00364                   E( I, J ) = ( HALF-SIN( CMPLX( J ) ) )*TWO
00365                ELSE
00366                   B( I, J ) = ZERO
00367                   E( I, J ) = ZERO
00368                END IF
00369    90       CONTINUE
00370   100    CONTINUE
00371 *
00372          DO 120 I = 1, M
00373             DO 110 J = 1, N
00374                R( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*TWENTY
00375                L( I, J ) = ( HALF-SIN( CMPLX( I+J ) ) )*TWENTY
00376   110       CONTINUE
00377   120    CONTINUE
00378 *
00379          IF( PRTYPE.EQ.3 ) THEN
00380             IF( QBLCKA.LE.1 )
00381      $         QBLCKA = 2
00382             DO 130 K = 1, M - 1, QBLCKA
00383                A( K+1, K+1 ) = A( K, K )
00384                A( K+1, K ) = -SIN( A( K, K+1 ) )
00385   130       CONTINUE
00386 *
00387             IF( QBLCKB.LE.1 )
00388      $         QBLCKB = 2
00389             DO 140 K = 1, N - 1, QBLCKB
00390                B( K+1, K+1 ) = B( K, K )
00391                B( K+1, K ) = -SIN( B( K, K+1 ) )
00392   140       CONTINUE
00393          END IF
00394 *
00395       ELSE IF( PRTYPE.EQ.4 ) THEN
00396          DO 160 I = 1, M
00397             DO 150 J = 1, M
00398                A( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*TWENTY
00399                D( I, J ) = ( HALF-SIN( CMPLX( I+J ) ) )*TWO
00400   150       CONTINUE
00401   160    CONTINUE
00402 *
00403          DO 180 I = 1, N
00404             DO 170 J = 1, N
00405                B( I, J ) = ( HALF-SIN( CMPLX( I+J ) ) )*TWENTY
00406                E( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*TWO
00407   170       CONTINUE
00408   180    CONTINUE
00409 *
00410          DO 200 I = 1, M
00411             DO 190 J = 1, N
00412                R( I, J ) = ( HALF-SIN( CMPLX( J / I ) ) )*TWENTY
00413                L( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*TWO
00414   190       CONTINUE
00415   200    CONTINUE
00416 *
00417       ELSE IF( PRTYPE.GE.5 ) THEN
00418          REEPS = HALF*TWO*TWENTY / ALPHA
00419          IMEPS = ( HALF-TWO ) / ALPHA
00420          DO 220 I = 1, M
00421             DO 210 J = 1, N
00422                R( I, J ) = ( HALF-SIN( CMPLX( I*J ) ) )*ALPHA / TWENTY
00423                L( I, J ) = ( HALF-SIN( CMPLX( I+J ) ) )*ALPHA / TWENTY
00424   210       CONTINUE
00425   220    CONTINUE
00426 *
00427          DO 230 I = 1, M
00428             D( I, I ) = ONE
00429   230    CONTINUE
00430 *
00431          DO 240 I = 1, M
00432             IF( I.LE.4 ) THEN
00433                A( I, I ) = ONE
00434                IF( I.GT.2 )
00435      $            A( I, I ) = ONE + REEPS
00436                IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
00437                   A( I, I+1 ) = IMEPS
00438                ELSE IF( I.GT.1 ) THEN
00439                   A( I, I-1 ) = -IMEPS
00440                END IF
00441             ELSE IF( I.LE.8 ) THEN
00442                IF( I.LE.6 ) THEN
00443                   A( I, I ) = REEPS
00444                ELSE
00445                   A( I, I ) = -REEPS
00446                END IF
00447                IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
00448                   A( I, I+1 ) = ONE
00449                ELSE IF( I.GT.1 ) THEN
00450                   A( I, I-1 ) = -ONE
00451                END IF
00452             ELSE
00453                A( I, I ) = ONE
00454                IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
00455                   A( I, I+1 ) = IMEPS*2
00456                ELSE IF( I.GT.1 ) THEN
00457                   A( I, I-1 ) = -IMEPS*2
00458                END IF
00459             END IF
00460   240    CONTINUE
00461 *
00462          DO 250 I = 1, N
00463             E( I, I ) = ONE
00464             IF( I.LE.4 ) THEN
00465                B( I, I ) = -ONE
00466                IF( I.GT.2 )
00467      $            B( I, I ) = ONE - REEPS
00468                IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
00469                   B( I, I+1 ) = IMEPS
00470                ELSE IF( I.GT.1 ) THEN
00471                   B( I, I-1 ) = -IMEPS
00472                END IF
00473             ELSE IF( I.LE.8 ) THEN
00474                IF( I.LE.6 ) THEN
00475                   B( I, I ) = REEPS
00476                ELSE
00477                   B( I, I ) = -REEPS
00478                END IF
00479                IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
00480                   B( I, I+1 ) = ONE + IMEPS
00481                ELSE IF( I.GT.1 ) THEN
00482                   B( I, I-1 ) = -ONE - IMEPS
00483                END IF
00484             ELSE
00485                B( I, I ) = ONE - REEPS
00486                IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
00487                   B( I, I+1 ) = IMEPS*2
00488                ELSE IF( I.GT.1 ) THEN
00489                   B( I, I-1 ) = -IMEPS*2
00490                END IF
00491             END IF
00492   250    CONTINUE
00493       END IF
00494 *
00495 *     Compute rhs (C, F)
00496 *
00497       CALL CGEMM( 'N', 'N', M, N, M, ONE, A, LDA, R, LDR, ZERO, C, LDC )
00498       CALL CGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, B, LDB, ONE, C, LDC )
00499       CALL CGEMM( 'N', 'N', M, N, M, ONE, D, LDD, R, LDR, ZERO, F, LDF )
00500       CALL CGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, E, LDE, ONE, F, LDF )
00501 *
00502 *     End of CLATM5
00503 *
00504       END
 All Files Functions