 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLATM5( 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 DOUBLE PRECISION 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 * DLATM5 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DOUBLE PRECISION ONE, ZERO, TWENTY, HALF, TWO 00180 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0, TWENTY = 2.0D+1, 00181 $ HALF = 0.5D+0, TWO = 2.0D+0 ) 00182 * .. 00183 * .. Local Scalars .. 00184 INTEGER I, J, K 00185 DOUBLE PRECISION IMEPS, REEPS 00186 * .. 00187 * .. Intrinsic Functions .. 00188 INTRINSIC DBLE, MOD, SIN 00189 * .. 00190 * .. External Subroutines .. 00191 EXTERNAL DGEMM 00192 * .. 00193 * .. Executable Statements .. 00194 * 00195 IF( PRTYPE.EQ.1 ) THEN 00196 DO 20 I = 1, M 00197 DO 10 J = 1, M 00198 IF( I.EQ.J ) THEN 00199 A( I, J ) = ONE 00200 D( I, J ) = ONE 00201 ELSE IF( I.EQ.J-1 ) THEN 00202 A( I, J ) = -ONE 00203 D( I, J ) = ZERO 00204 ELSE 00205 A( I, J ) = ZERO 00206 D( I, J ) = ZERO 00207 END IF 00208 10 CONTINUE 00209 20 CONTINUE 00210 * 00211 DO 40 I = 1, N 00212 DO 30 J = 1, N 00213 IF( I.EQ.J ) THEN 00214 B( I, J ) = ONE - ALPHA 00215 E( I, J ) = ONE 00216 ELSE IF( I.EQ.J-1 ) THEN 00217 B( I, J ) = ONE 00218 E( I, J ) = ZERO 00219 ELSE 00220 B( I, J ) = ZERO 00221 E( I, J ) = ZERO 00222 END IF 00223 30 CONTINUE 00224 40 CONTINUE 00225 * 00226 DO 60 I = 1, M 00227 DO 50 J = 1, N 00228 R( I, J ) = ( HALF-SIN( DBLE( I / J ) ) )*TWENTY 00229 L( I, J ) = R( I, J ) 00230 50 CONTINUE 00231 60 CONTINUE 00232 * 00233 ELSE IF( PRTYPE.EQ.2 .OR. PRTYPE.EQ.3 ) THEN 00234 DO 80 I = 1, M 00235 DO 70 J = 1, M 00236 IF( I.LE.J ) THEN 00237 A( I, J ) = ( HALF-SIN( DBLE( I ) ) )*TWO 00238 D( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*TWO 00239 ELSE 00240 A( I, J ) = ZERO 00241 D( I, J ) = ZERO 00242 END IF 00243 70 CONTINUE 00244 80 CONTINUE 00245 * 00246 DO 100 I = 1, N 00247 DO 90 J = 1, N 00248 IF( I.LE.J ) THEN 00249 B( I, J ) = ( HALF-SIN( DBLE( I+J ) ) )*TWO 00250 E( I, J ) = ( HALF-SIN( DBLE( J ) ) )*TWO 00251 ELSE 00252 B( I, J ) = ZERO 00253 E( I, J ) = ZERO 00254 END IF 00255 90 CONTINUE 00256 100 CONTINUE 00257 * 00258 DO 120 I = 1, M 00259 DO 110 J = 1, N 00260 R( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*TWENTY 00261 L( I, J ) = ( HALF-SIN( DBLE( I+J ) ) )*TWENTY 00262 110 CONTINUE 00263 120 CONTINUE 00264 * 00265 IF( PRTYPE.EQ.3 ) THEN 00266 IF( QBLCKA.LE.1 ) 00267 $ QBLCKA = 2 00268 DO 130 K = 1, M - 1, QBLCKA 00269 A( K+1, K+1 ) = A( K, K ) 00270 A( K+1, K ) = -SIN( A( K, K+1 ) ) 00271 130 CONTINUE 00272 * 00273 IF( QBLCKB.LE.1 ) 00274 $ QBLCKB = 2 00275 DO 140 K = 1, N - 1, QBLCKB 00276 B( K+1, K+1 ) = B( K, K ) 00277 B( K+1, K ) = -SIN( B( K, K+1 ) ) 00278 140 CONTINUE 00279 END IF 00280 * 00281 ELSE IF( PRTYPE.EQ.4 ) THEN 00282 DO 160 I = 1, M 00283 DO 150 J = 1, M 00284 A( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*TWENTY 00285 D( I, J ) = ( HALF-SIN( DBLE( I+J ) ) )*TWO 00286 150 CONTINUE 00287 160 CONTINUE 00288 * 00289 DO 180 I = 1, N 00290 DO 170 J = 1, N 00291 B( I, J ) = ( HALF-SIN( DBLE( I+J ) ) )*TWENTY 00292 E( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*TWO 00293 170 CONTINUE 00294 180 CONTINUE 00295 * 00296 DO 200 I = 1, M 00297 DO 190 J = 1, N 00298 R( I, J ) = ( HALF-SIN( DBLE( J / I ) ) )*TWENTY 00299 L( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*TWO 00300 190 CONTINUE 00301 200 CONTINUE 00302 * 00303 ELSE IF( PRTYPE.GE.5 ) THEN 00304 REEPS = HALF*TWO*TWENTY / ALPHA 00305 IMEPS = ( HALF-TWO ) / ALPHA 00306 DO 220 I = 1, M 00307 DO 210 J = 1, N 00308 R( I, J ) = ( HALF-SIN( DBLE( I*J ) ) )*ALPHA / TWENTY 00309 L( I, J ) = ( HALF-SIN( DBLE( I+J ) ) )*ALPHA / TWENTY 00310 210 CONTINUE 00311 220 CONTINUE 00312 * 00313 DO 230 I = 1, M 00314 D( I, I ) = ONE 00315 230 CONTINUE 00316 * 00317 DO 240 I = 1, M 00318 IF( I.LE.4 ) THEN 00319 A( I, I ) = ONE 00320 IF( I.GT.2 ) 00321 $ A( I, I ) = ONE + REEPS 00322 IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN 00323 A( I, I+1 ) = IMEPS 00324 ELSE IF( I.GT.1 ) THEN 00325 A( I, I-1 ) = -IMEPS 00326 END IF 00327 ELSE IF( I.LE.8 ) THEN 00328 IF( I.LE.6 ) THEN 00329 A( I, I ) = REEPS 00330 ELSE 00331 A( I, I ) = -REEPS 00332 END IF 00333 IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN 00334 A( I, I+1 ) = ONE 00335 ELSE IF( I.GT.1 ) THEN 00336 A( I, I-1 ) = -ONE 00337 END IF 00338 ELSE 00339 A( I, I ) = ONE 00340 IF( MOD( I, 2 ).NE.0 .AND. I.LT.M ) THEN 00341 A( I, I+1 ) = IMEPS*2 00342 ELSE IF( I.GT.1 ) THEN 00343 A( I, I-1 ) = -IMEPS*2 00344 END IF 00345 END IF 00346 240 CONTINUE 00347 * 00348 DO 250 I = 1, N 00349 E( I, I ) = ONE 00350 IF( I.LE.4 ) THEN 00351 B( I, I ) = -ONE 00352 IF( I.GT.2 ) 00353 $ B( I, I ) = ONE - REEPS 00354 IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN 00355 B( I, I+1 ) = IMEPS 00356 ELSE IF( I.GT.1 ) THEN 00357 B( I, I-1 ) = -IMEPS 00358 END IF 00359 ELSE IF( I.LE.8 ) THEN 00360 IF( I.LE.6 ) THEN 00361 B( I, I ) = REEPS 00362 ELSE 00363 B( I, I ) = -REEPS 00364 END IF 00365 IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN 00366 B( I, I+1 ) = ONE + IMEPS 00367 ELSE IF( I.GT.1 ) THEN 00368 B( I, I-1 ) = -ONE - IMEPS 00369 END IF 00370 ELSE 00371 B( I, I ) = ONE - REEPS 00372 IF( MOD( I, 2 ).NE.0 .AND. I.LT.N ) THEN 00373 B( I, I+1 ) = IMEPS*2 00374 ELSE IF( I.GT.1 ) THEN 00375 B( I, I-1 ) = -IMEPS*2 00376 END IF 00377 END IF 00378 250 CONTINUE 00379 END IF 00380 * 00381 * Compute rhs (C, F) 00382 * 00383 CALL DGEMM( 'N', 'N', M, N, M, ONE, A, LDA, R, LDR, ZERO, C, LDC ) 00384 CALL DGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, B, LDB, ONE, C, LDC ) 00385 CALL DGEMM( 'N', 'N', M, N, M, ONE, D, LDD, R, LDR, ZERO, F, LDF ) 00386 CALL DGEMM( 'N', 'N', M, N, N, -ONE, L, LDL, E, LDE, ONE, F, LDF ) 00387 * 00388 * End of DLATM5 00389 * 00390 END
1.7.3 
