 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZLATM5( 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 DOUBLE PRECISION ALPHA 00013 * .. 00014 * .. Array Arguments .. 00015 COMPLEX*16 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 * ZLATM5 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*16 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*16 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*16 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*16 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*16 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*16 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*16 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*16 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) DOUBLE PRECISION 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*16 ONE, TWO, ZERO, HALF, TWENTY 00180 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 00181 $ TWO = ( 2.0D+0, 0.0D+0 ), 00182 $ ZERO = ( 0.0D+0, 0.0D+0 ), 00183 $ HALF = ( 0.5D+0, 0.0D+0 ), 00184 $ TWENTY = ( 2.0D+1, 0.0D+0 ) ) 00185 * .. 00186 * .. Local Scalars .. 00187 INTEGER I, J, K 00188 COMPLEX*16 IMEPS, REEPS 00189 * .. 00190 * .. Intrinsic Functions .. 00191 INTRINSIC DCMPLX, MOD, SIN 00192 * .. 00193 * .. External Subroutines .. 00194 EXTERNAL ZGEMM 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( DCMPLX( 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( DCMPLX( I ) ) )*TWO 00241 D( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( I+J ) ) )*TWO 00253 E( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( I*J ) ) )*TWENTY 00264 L( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( I*J ) ) )*TWENTY 00288 D( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( I+J ) ) )*TWENTY 00295 E( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( J / I ) ) )*TWENTY 00302 L( I, J ) = ( HALF-SIN( DCMPLX( 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( DCMPLX( I*J ) ) )*ALPHA / TWENTY 00312 L( I, J ) = ( HALF-SIN( DCMPLX( 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 ZGEMM( 'N', 'N', M, N, M, ONE, A, LDA, R, LDR, ZERO, C, LDC ) 00387 CALL ZGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, B, LDB, ONE, C, LDC ) 00388 CALL ZGEMM( 'N', 'N', M, N, M, ONE, D, LDD, R, LDR, ZERO, F, LDF ) 00389 CALL ZGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, E, LDE, ONE, F, LDF ) 00390 * 00391 * End of ZLATM5 00392 * 00393 END
1.7.3 
