![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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