LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slatm5.f
Go to the documentation of this file.
00001 *> \brief \b SLATM5
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 SLATM5( 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 *       REAL               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 *> SLATM5 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 real_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 SLATM5( 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       REAL               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       REAL               ONE, ZERO, TWENTY, HALF, TWO
00291       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0, TWENTY = 2.0E+1,
00292      $                   HALF = 0.5E+0, TWO = 2.0E+0 )
00293 *     ..
00294 *     .. Local Scalars ..
00295       INTEGER            I, J, K
00296       REAL               IMEPS, REEPS
00297 *     ..
00298 *     .. Intrinsic Functions ..
00299       INTRINSIC          MOD, REAL, SIN
00300 *     ..
00301 *     .. External Subroutines ..
00302       EXTERNAL           SGEMM
00303 *     ..
00304 *     .. Executable Statements ..
00305 *
00306       IF( PRTYPE.EQ.1 ) THEN
00307          DO 20 I = 1, M
00308             DO 10 J = 1, M
00309                IF( I.EQ.J ) THEN
00310                   A( I, J ) = ONE
00311                   D( I, J ) = ONE
00312                ELSE IF( I.EQ.J-1 ) THEN
00313                   A( I, J ) = -ONE
00314                   D( I, J ) = ZERO
00315                ELSE
00316                   A( I, J ) = ZERO
00317                   D( I, J ) = ZERO
00318                END IF
00319    10       CONTINUE
00320    20    CONTINUE
00321 *
00322          DO 40 I = 1, N
00323             DO 30 J = 1, N
00324                IF( I.EQ.J ) THEN
00325                   B( I, J ) = ONE - ALPHA
00326                   E( I, J ) = ONE
00327                ELSE IF( I.EQ.J-1 ) THEN
00328                   B( I, J ) = ONE
00329                   E( I, J ) = ZERO
00330                ELSE
00331                   B( I, J ) = ZERO
00332                   E( I, J ) = ZERO
00333                END IF
00334    30       CONTINUE
00335    40    CONTINUE
00336 *
00337          DO 60 I = 1, M
00338             DO 50 J = 1, N
00339                R( I, J ) = ( HALF-SIN( REAL( I / J ) ) )*TWENTY
00340                L( I, J ) = R( I, J )
00341    50       CONTINUE
00342    60    CONTINUE
00343 *
00344       ELSE IF( PRTYPE.EQ.2 .OR. PRTYPE.EQ.3 ) THEN
00345          DO 80 I = 1, M
00346             DO 70 J = 1, M
00347                IF( I.LE.J ) THEN
00348                   A( I, J ) = ( HALF-SIN( REAL( I ) ) )*TWO
00349                   D( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*TWO
00350                ELSE
00351                   A( I, J ) = ZERO
00352                   D( I, J ) = ZERO
00353                END IF
00354    70       CONTINUE
00355    80    CONTINUE
00356 *
00357          DO 100 I = 1, N
00358             DO 90 J = 1, N
00359                IF( I.LE.J ) THEN
00360                   B( I, J ) = ( HALF-SIN( REAL( I+J ) ) )*TWO
00361                   E( I, J ) = ( HALF-SIN( REAL( J ) ) )*TWO
00362                ELSE
00363                   B( I, J ) = ZERO
00364                   E( I, J ) = ZERO
00365                END IF
00366    90       CONTINUE
00367   100    CONTINUE
00368 *
00369          DO 120 I = 1, M
00370             DO 110 J = 1, N
00371                R( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*TWENTY
00372                L( I, J ) = ( HALF-SIN( REAL( I+J ) ) )*TWENTY
00373   110       CONTINUE
00374   120    CONTINUE
00375 *
00376          IF( PRTYPE.EQ.3 ) THEN
00377             IF( QBLCKA.LE.1 )
00378      $         QBLCKA = 2
00379             DO 130 K = 1, M - 1, QBLCKA
00380                A( K+1, K+1 ) = A( K, K )
00381                A( K+1, K ) = -SIN( A( K, K+1 ) )
00382   130       CONTINUE
00383 *
00384             IF( QBLCKB.LE.1 )
00385      $         QBLCKB = 2
00386             DO 140 K = 1, N - 1, QBLCKB
00387                B( K+1, K+1 ) = B( K, K )
00388                B( K+1, K ) = -SIN( B( K, K+1 ) )
00389   140       CONTINUE
00390          END IF
00391 *
00392       ELSE IF( PRTYPE.EQ.4 ) THEN
00393          DO 160 I = 1, M
00394             DO 150 J = 1, M
00395                A( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*TWENTY
00396                D( I, J ) = ( HALF-SIN( REAL( I+J ) ) )*TWO
00397   150       CONTINUE
00398   160    CONTINUE
00399 *
00400          DO 180 I = 1, N
00401             DO 170 J = 1, N
00402                B( I, J ) = ( HALF-SIN( REAL( I+J ) ) )*TWENTY
00403                E( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*TWO
00404   170       CONTINUE
00405   180    CONTINUE
00406 *
00407          DO 200 I = 1, M
00408             DO 190 J = 1, N
00409                R( I, J ) = ( HALF-SIN( REAL( J / I ) ) )*TWENTY
00410                L( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*TWO
00411   190       CONTINUE
00412   200    CONTINUE
00413 *
00414       ELSE IF( PRTYPE.GE.5 ) THEN
00415          REEPS = HALF*TWO*TWENTY / ALPHA
00416          IMEPS = ( HALF-TWO ) / ALPHA
00417          DO 220 I = 1, M
00418             DO 210 J = 1, N
00419                R( I, J ) = ( HALF-SIN( REAL( I*J ) ) )*ALPHA / TWENTY
00420                L( I, J ) = ( HALF-SIN( REAL( I+J ) ) )*ALPHA / TWENTY
00421   210       CONTINUE
00422   220    CONTINUE
00423 *
00424          DO 230 I = 1, M
00425             D( I, I ) = ONE
00426   230    CONTINUE
00427 *
00428          DO 240 I = 1, M
00429             IF( I.LE.4 ) THEN
00430                A( I, I ) = ONE
00431                IF( I.GT.2 )
00432      $            A( I, I ) = ONE + REEPS
00433                IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
00434                   A( I, I+1 ) = IMEPS
00435                ELSE IF( I.GT.1 ) THEN
00436                   A( I, I-1 ) = -IMEPS
00437                END IF
00438             ELSE IF( I.LE.8 ) THEN
00439                IF( I.LE.6 ) THEN
00440                   A( I, I ) = REEPS
00441                ELSE
00442                   A( I, I ) = -REEPS
00443                END IF
00444                IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
00445                   A( I, I+1 ) = ONE
00446                ELSE IF( I.GT.1 ) THEN
00447                   A( I, I-1 ) = -ONE
00448                END IF
00449             ELSE
00450                A( I, I ) = ONE
00451                IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN
00452                   A( I, I+1 ) = IMEPS*2
00453                ELSE IF( I.GT.1 ) THEN
00454                   A( I, I-1 ) = -IMEPS*2
00455                END IF
00456             END IF
00457   240    CONTINUE
00458 *
00459          DO 250 I = 1, N
00460             E( I, I ) = ONE
00461             IF( I.LE.4 ) THEN
00462                B( I, I ) = -ONE
00463                IF( I.GT.2 )
00464      $            B( I, I ) = ONE - REEPS
00465                IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
00466                   B( I, I+1 ) = IMEPS
00467                ELSE IF( I.GT.1 ) THEN
00468                   B( I, I-1 ) = -IMEPS
00469                END IF
00470             ELSE IF( I.LE.8 ) THEN
00471                IF( I.LE.6 ) THEN
00472                   B( I, I ) = REEPS
00473                ELSE
00474                   B( I, I ) = -REEPS
00475                END IF
00476                IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
00477                   B( I, I+1 ) = ONE + IMEPS
00478                ELSE IF( I.GT.1 ) THEN
00479                   B( I, I-1 ) = -ONE - IMEPS
00480                END IF
00481             ELSE
00482                B( I, I ) = ONE - REEPS
00483                IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN
00484                   B( I, I+1 ) = IMEPS*2
00485                ELSE IF( I.GT.1 ) THEN
00486                   B( I, I-1 ) = -IMEPS*2
00487                END IF
00488             END IF
00489   250    CONTINUE
00490       END IF
00491 *
00492 *     Compute rhs (C, F)
00493 *
00494       CALL SGEMM( 'N', 'N', M, N, M, ONE, A, LDA, R, LDR, ZERO, C, LDC )
00495       CALL SGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, B, LDB, ONE, C, LDC )
00496       CALL SGEMM( 'N', 'N', M, N, M, ONE, D, LDD, R, LDR, ZERO, F, LDF )
00497       CALL SGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, E, LDE, ONE, F, LDF )
00498 *
00499 *     End of SLATM5
00500 *
00501       END
 All Files Functions