![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLATM5 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 ZLATM5( 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 * DOUBLE PRECISION ALPHA 00019 * .. 00020 * .. Array Arguments .. 00021 * COMPLEX*16 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 *> ZLATM5 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*16 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*16 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*16 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*16 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*16 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*16 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*16 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*16 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 DOUBLE PRECISION 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 complex16_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 ZLATM5( 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 DOUBLE PRECISION ALPHA 00280 * .. 00281 * .. Array Arguments .. 00282 COMPLEX*16 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*16 ONE, TWO, ZERO, HALF, TWENTY 00291 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 00292 $ TWO = ( 2.0D+0, 0.0D+0 ), 00293 $ ZERO = ( 0.0D+0, 0.0D+0 ), 00294 $ HALF = ( 0.5D+0, 0.0D+0 ), 00295 $ TWENTY = ( 2.0D+1, 0.0D+0 ) ) 00296 * .. 00297 * .. Local Scalars .. 00298 INTEGER I, J, K 00299 COMPLEX*16 IMEPS, REEPS 00300 * .. 00301 * .. Intrinsic Functions .. 00302 INTRINSIC DCMPLX, MOD, SIN 00303 * .. 00304 * .. External Subroutines .. 00305 EXTERNAL ZGEMM 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( DCMPLX( 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( DCMPLX( I ) ) )*TWO 00352 D( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( I+J ) ) )*TWO 00364 E( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( I*J ) ) )*TWENTY 00375 L( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( I*J ) ) )*TWENTY 00399 D( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( I+J ) ) )*TWENTY 00406 E( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( J / I ) ) )*TWENTY 00413 L( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( I*J ) ) )*ALPHA / TWENTY 00423 L( I, J ) = ( HALF-SIN( DCMPLX( 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 ZGEMM( 'N', 'N', M, N, M, ONE, A, LDA, R, LDR, ZERO, C, LDC ) 00498 CALL ZGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, B, LDB, ONE, C, LDC ) 00499 CALL ZGEMM( 'N', 'N', M, N, M, ONE, D, LDD, R, LDR, ZERO, F, LDF ) 00500 CALL ZGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, E, LDE, ONE, F, LDF ) 00501 * 00502 * End of ZLATM5 00503 * 00504 END