 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CLATME( N, DIST, ISEED, D, MODE, COND, DMAX, EI, 00002 $ RSIGN, 00003 $ UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, 00004 $ A, 00005 $ LDA, WORK, INFO ) 00006 * 00007 * -- LAPACK test routine (version 3.1) -- 00008 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00009 * June 2010 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER DIST, EI, RSIGN, SIM, UPPER 00013 INTEGER INFO, KL, KU, LDA, MODE, MODES, N 00014 REAL ANORM, COND, CONDS 00015 COMPLEX DMAX 00016 * .. 00017 * .. Array Arguments .. 00018 INTEGER ISEED( 4 ) 00019 REAL DS( * ) 00020 COMPLEX A( LDA, * ), D( * ), WORK( * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * CLATME generates random non-symmetric square matrices with 00027 * specified eigenvalues for testing LAPACK programs. 00028 * 00029 * CLATME operates by applying the following sequence of 00030 * operations: 00031 * 00032 * 1. Set the diagonal to D, where D may be input or 00033 * computed according to MODE, COND, DMAX, and RSIGN 00034 * as described below. 00035 * 00036 * 2. If UPPER='T', the upper triangle of A is set to random values 00037 * out of distribution DIST. 00038 * 00039 * 3. If SIM='T', A is multiplied on the left by a random matrix 00040 * X, whose singular values are specified by DS, MODES, and 00041 * CONDS, and on the right by X inverse. 00042 * 00043 * 4. If KL < N-1, the lower bandwidth is reduced to KL using 00044 * Householder transformations. If KU < N-1, the upper 00045 * bandwidth is reduced to KU. 00046 * 00047 * 5. If ANORM is not negative, the matrix is scaled to have 00048 * maximum-element-norm ANORM. 00049 * 00050 * (Note: since the matrix cannot be reduced beyond Hessenberg form, 00051 * no packing options are available.) 00052 * 00053 * Arguments 00054 * ========= 00055 * 00056 * N (input) INTEGER 00057 * The number of columns (or rows) of A. Not modified. 00058 * 00059 * DIST (input) CHARACTER*1 00060 * On entry, DIST specifies the type of distribution to be used 00061 * to generate the random eigen-/singular values, and on the 00062 * upper triangle (see UPPER). 00063 * 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform ) 00064 * 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric ) 00065 * 'N' => NORMAL( 0, 1 ) ( 'N' for normal ) 00066 * 'D' => uniform on the complex disc |z| < 1. 00067 * Not modified. 00068 * 00069 * ISEED (input/output) INTEGER array, dimension ( 4 ) 00070 * On entry ISEED specifies the seed of the random number 00071 * generator. They should lie between 0 and 4095 inclusive, 00072 * and ISEED(4) should be odd. The random number generator 00073 * uses a linear congruential sequence limited to small 00074 * integers, and so should produce machine independent 00075 * random numbers. The values of ISEED are changed on 00076 * exit, and can be used in the next call to CLATME 00077 * to continue the same random number sequence. 00078 * Changed on exit. 00079 * 00080 * D (input/output) COMPLEX array, dimension ( N ) 00081 * This array is used to specify the eigenvalues of A. If 00082 * MODE=0, then D is assumed to contain the eigenvalues 00083 * otherwise they will be computed according to MODE, COND, 00084 * DMAX, and RSIGN and placed in D. 00085 * Modified if MODE is nonzero. 00086 * 00087 * MODE (input) INTEGER 00088 * On entry this describes how the eigenvalues are to 00089 * be specified: 00090 * MODE = 0 means use D as input 00091 * MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND 00092 * MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND 00093 * MODE = 3 sets D(I)=COND**(-(I-1)/(N-1)) 00094 * MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND) 00095 * MODE = 5 sets D to random numbers in the range 00096 * ( 1/COND , 1 ) such that their logarithms 00097 * are uniformly distributed. 00098 * MODE = 6 set D to random numbers from same distribution 00099 * as the rest of the matrix. 00100 * MODE < 0 has the same meaning as ABS(MODE), except that 00101 * the order of the elements of D is reversed. 00102 * Thus if MODE is between 1 and 4, D has entries ranging 00103 * from 1 to 1/COND, if between -1 and -4, D has entries 00104 * ranging from 1/COND to 1, 00105 * Not modified. 00106 * 00107 * COND (input) REAL 00108 * On entry, this is used as described under MODE above. 00109 * If used, it must be >= 1. Not modified. 00110 * 00111 * DMAX (input) COMPLEX 00112 * If MODE is neither -6, 0 nor 6, the contents of D, as 00113 * computed according to MODE and COND, will be scaled by 00114 * DMAX / max(abs(D(i))). Note that DMAX need not be 00115 * positive or real: if DMAX is negative or complex (or zero), 00116 * D will be scaled by a negative or complex number (or zero). 00117 * If RSIGN='F' then the largest (absolute) eigenvalue will be 00118 * equal to DMAX. 00119 * Not modified. 00120 * 00121 * EI (input) CHARACTER*1 array, dimension ( N ) 00122 * (ignored) 00123 * Not modified. 00124 * 00125 * RSIGN (input) CHARACTER*1 00126 * If MODE is not 0, 6, or -6, and RSIGN='T', then the 00127 * elements of D, as computed according to MODE and COND, will 00128 * be multiplied by a random complex number from the unit 00129 * circle |z| = 1. If RSIGN='F', they will not be. RSIGN may 00130 * only have the values 'T' or 'F'. 00131 * Not modified. 00132 * 00133 * UPPER (input) CHARACTER*1 00134 * If UPPER='T', then the elements of A above the diagonal 00135 * will be set to random numbers out of DIST. If UPPER='F', 00136 * they will not. UPPER may only have the values 'T' or 'F'. 00137 * Not modified. 00138 * 00139 * SIM (input) CHARACTER*1 00140 * If SIM='T', then A will be operated on by a "similarity 00141 * transform", i.e., multiplied on the left by a matrix X and 00142 * on the right by X inverse. X = U S V, where U and V are 00143 * random unitary matrices and S is a (diagonal) matrix of 00144 * singular values specified by DS, MODES, and CONDS. If 00145 * SIM='F', then A will not be transformed. 00146 * Not modified. 00147 * 00148 * DS (input/output) REAL array, dimension ( N ) 00149 * This array is used to specify the singular values of X, 00150 * in the same way that D specifies the eigenvalues of A. 00151 * If MODE=0, the DS contains the singular values, which 00152 * may not be zero. 00153 * Modified if MODE is nonzero. 00154 * 00155 * MODES (input) INTEGER 00156 * 00157 * CONDS (input) REAL 00158 * Similar to MODE and COND, but for specifying the diagonal 00159 * of S. MODES=-6 and +6 are not allowed (since they would 00160 * result in randomly ill-conditioned eigenvalues.) 00161 * 00162 * KL (input) INTEGER 00163 * This specifies the lower bandwidth of the matrix. KL=1 00164 * specifies upper Hessenberg form. If KL is at least N-1, 00165 * then A will have full lower bandwidth. 00166 * Not modified. 00167 * 00168 * KU (input) INTEGER 00169 * This specifies the upper bandwidth of the matrix. KU=1 00170 * specifies lower Hessenberg form. If KU is at least N-1, 00171 * then A will have full upper bandwidth; if KU and KL 00172 * are both at least N-1, then A will be dense. Only one of 00173 * KU and KL may be less than N-1. 00174 * Not modified. 00175 * 00176 * ANORM (input) REAL 00177 * If ANORM is not negative, then A will be scaled by a non- 00178 * negative real number to make the maximum-element-norm of A 00179 * to be ANORM. 00180 * Not modified. 00181 * 00182 * A (output) COMPLEX array, dimension ( LDA, N ) 00183 * On exit A is the desired test matrix. 00184 * Modified. 00185 * 00186 * LDA (input) INTEGER 00187 * LDA specifies the first dimension of A as declared in the 00188 * calling program. LDA must be at least M. 00189 * Not modified. 00190 * 00191 * WORK (workspace) COMPLEX array, dimension ( 3*N ) 00192 * Workspace. 00193 * Modified. 00194 * 00195 * INFO (output) INTEGER 00196 * Error code. On exit, INFO will be set to one of the 00197 * following values: 00198 * 0 => normal return 00199 * -1 => N negative 00200 * -2 => DIST illegal string 00201 * -5 => MODE not in range -6 to 6 00202 * -6 => COND less than 1.0, and MODE neither -6, 0 nor 6 00203 * -9 => RSIGN is not 'T' or 'F' 00204 * -10 => UPPER is not 'T' or 'F' 00205 * -11 => SIM is not 'T' or 'F' 00206 * -12 => MODES=0 and DS has a zero singular value. 00207 * -13 => MODES is not in the range -5 to 5. 00208 * -14 => MODES is nonzero and CONDS is less than 1. 00209 * -15 => KL is less than 1. 00210 * -16 => KU is less than 1, or KL and KU are both less than 00211 * N-1. 00212 * -19 => LDA is less than M. 00213 * 1 => Error return from CLATM1 (computing D) 00214 * 2 => Cannot scale to DMAX (max. eigenvalue is 0) 00215 * 3 => Error return from SLATM1 (computing DS) 00216 * 4 => Error return from CLARGE 00217 * 5 => Zero singular value from SLATM1. 00218 * 00219 * ===================================================================== 00220 * 00221 * .. Parameters .. 00222 REAL ZERO 00223 PARAMETER ( ZERO = 0.0E+0 ) 00224 REAL ONE 00225 PARAMETER ( ONE = 1.0E+0 ) 00226 COMPLEX CZERO 00227 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) 00228 COMPLEX CONE 00229 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00230 * .. 00231 * .. Local Scalars .. 00232 LOGICAL BADS 00233 INTEGER I, IC, ICOLS, IDIST, IINFO, IR, IROWS, IRSIGN, 00234 $ ISIM, IUPPER, J, JC, JCR 00235 REAL RALPHA, TEMP 00236 COMPLEX ALPHA, TAU, XNORMS 00237 * .. 00238 * .. Local Arrays .. 00239 REAL TEMPA( 1 ) 00240 * .. 00241 * .. External Functions .. 00242 LOGICAL LSAME 00243 REAL CLANGE 00244 COMPLEX CLARND 00245 EXTERNAL LSAME, CLANGE, CLARND 00246 * .. 00247 * .. External Subroutines .. 00248 EXTERNAL CCOPY, CGEMV, CGERC, CLACGV, CLARFG, CLARGE, 00249 $ CLARNV, CLATM1, CLASET, CSCAL, CSSCAL, SLATM1, 00250 $ XERBLA 00251 * .. 00252 * .. Intrinsic Functions .. 00253 INTRINSIC ABS, CONJG, MAX, MOD 00254 * .. 00255 * .. Executable Statements .. 00256 * 00257 * 1) Decode and Test the input parameters. 00258 * Initialize flags & seed. 00259 * 00260 INFO = 0 00261 * 00262 * Quick return if possible 00263 * 00264 IF( N.EQ.0 ) 00265 $ RETURN 00266 * 00267 * Decode DIST 00268 * 00269 IF( LSAME( DIST, 'U' ) ) THEN 00270 IDIST = 1 00271 ELSE IF( LSAME( DIST, 'S' ) ) THEN 00272 IDIST = 2 00273 ELSE IF( LSAME( DIST, 'N' ) ) THEN 00274 IDIST = 3 00275 ELSE IF( LSAME( DIST, 'D' ) ) THEN 00276 IDIST = 4 00277 ELSE 00278 IDIST = -1 00279 END IF 00280 * 00281 * Decode RSIGN 00282 * 00283 IF( LSAME( RSIGN, 'T' ) ) THEN 00284 IRSIGN = 1 00285 ELSE IF( LSAME( RSIGN, 'F' ) ) THEN 00286 IRSIGN = 0 00287 ELSE 00288 IRSIGN = -1 00289 END IF 00290 * 00291 * Decode UPPER 00292 * 00293 IF( LSAME( UPPER, 'T' ) ) THEN 00294 IUPPER = 1 00295 ELSE IF( LSAME( UPPER, 'F' ) ) THEN 00296 IUPPER = 0 00297 ELSE 00298 IUPPER = -1 00299 END IF 00300 * 00301 * Decode SIM 00302 * 00303 IF( LSAME( SIM, 'T' ) ) THEN 00304 ISIM = 1 00305 ELSE IF( LSAME( SIM, 'F' ) ) THEN 00306 ISIM = 0 00307 ELSE 00308 ISIM = -1 00309 END IF 00310 * 00311 * Check DS, if MODES=0 and ISIM=1 00312 * 00313 BADS = .FALSE. 00314 IF( MODES.EQ.0 .AND. ISIM.EQ.1 ) THEN 00315 DO 10 J = 1, N 00316 IF( DS( J ).EQ.ZERO ) 00317 $ BADS = .TRUE. 00318 10 CONTINUE 00319 END IF 00320 * 00321 * Set INFO if an error 00322 * 00323 IF( N.LT.0 ) THEN 00324 INFO = -1 00325 ELSE IF( IDIST.EQ.-1 ) THEN 00326 INFO = -2 00327 ELSE IF( ABS( MODE ).GT.6 ) THEN 00328 INFO = -5 00329 ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE ) 00330 $ THEN 00331 INFO = -6 00332 ELSE IF( IRSIGN.EQ.-1 ) THEN 00333 INFO = -9 00334 ELSE IF( IUPPER.EQ.-1 ) THEN 00335 INFO = -10 00336 ELSE IF( ISIM.EQ.-1 ) THEN 00337 INFO = -11 00338 ELSE IF( BADS ) THEN 00339 INFO = -12 00340 ELSE IF( ISIM.EQ.1 .AND. ABS( MODES ).GT.5 ) THEN 00341 INFO = -13 00342 ELSE IF( ISIM.EQ.1 .AND. MODES.NE.0 .AND. CONDS.LT.ONE ) THEN 00343 INFO = -14 00344 ELSE IF( KL.LT.1 ) THEN 00345 INFO = -15 00346 ELSE IF( KU.LT.1 .OR. ( KU.LT.N-1 .AND. KL.LT.N-1 ) ) THEN 00347 INFO = -16 00348 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00349 INFO = -19 00350 END IF 00351 * 00352 IF( INFO.NE.0 ) THEN 00353 CALL XERBLA( 'CLATME', -INFO ) 00354 RETURN 00355 END IF 00356 * 00357 * Initialize random number generator 00358 * 00359 DO 20 I = 1, 4 00360 ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 ) 00361 20 CONTINUE 00362 * 00363 IF( MOD( ISEED( 4 ), 2 ).NE.1 ) 00364 $ ISEED( 4 ) = ISEED( 4 ) + 1 00365 * 00366 * 2) Set up diagonal of A 00367 * 00368 * Compute D according to COND and MODE 00369 * 00370 CALL CLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, N, IINFO ) 00371 IF( IINFO.NE.0 ) THEN 00372 INFO = 1 00373 RETURN 00374 END IF 00375 IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN 00376 * 00377 * Scale by DMAX 00378 * 00379 TEMP = ABS( D( 1 ) ) 00380 DO 30 I = 2, N 00381 TEMP = MAX( TEMP, ABS( D( I ) ) ) 00382 30 CONTINUE 00383 * 00384 IF( TEMP.GT.ZERO ) THEN 00385 ALPHA = DMAX / TEMP 00386 ELSE 00387 INFO = 2 00388 RETURN 00389 END IF 00390 * 00391 CALL CSCAL( N, ALPHA, D, 1 ) 00392 * 00393 END IF 00394 * 00395 CALL CLASET( 'Full', N, N, CZERO, CZERO, A, LDA ) 00396 CALL CCOPY( N, D, 1, A, LDA+1 ) 00397 * 00398 * 3) If UPPER='T', set upper triangle of A to random numbers. 00399 * 00400 IF( IUPPER.NE.0 ) THEN 00401 DO 40 JC = 2, N 00402 CALL CLARNV( IDIST, ISEED, JC-1, A( 1, JC ) ) 00403 40 CONTINUE 00404 END IF 00405 * 00406 * 4) If SIM='T', apply similarity transformation. 00407 * 00408 * -1 00409 * Transform is X A X , where X = U S V, thus 00410 * 00411 * it is U S V A V' (1/S) U' 00412 * 00413 IF( ISIM.NE.0 ) THEN 00414 * 00415 * Compute S (singular values of the eigenvector matrix) 00416 * according to CONDS and MODES 00417 * 00418 CALL SLATM1( MODES, CONDS, 0, 0, ISEED, DS, N, IINFO ) 00419 IF( IINFO.NE.0 ) THEN 00420 INFO = 3 00421 RETURN 00422 END IF 00423 * 00424 * Multiply by V and V' 00425 * 00426 CALL CLARGE( N, A, LDA, ISEED, WORK, IINFO ) 00427 IF( IINFO.NE.0 ) THEN 00428 INFO = 4 00429 RETURN 00430 END IF 00431 * 00432 * Multiply by S and (1/S) 00433 * 00434 DO 50 J = 1, N 00435 CALL CSSCAL( N, DS( J ), A( J, 1 ), LDA ) 00436 IF( DS( J ).NE.ZERO ) THEN 00437 CALL CSSCAL( N, ONE / DS( J ), A( 1, J ), 1 ) 00438 ELSE 00439 INFO = 5 00440 RETURN 00441 END IF 00442 50 CONTINUE 00443 * 00444 * Multiply by U and U' 00445 * 00446 CALL CLARGE( N, A, LDA, ISEED, WORK, IINFO ) 00447 IF( IINFO.NE.0 ) THEN 00448 INFO = 4 00449 RETURN 00450 END IF 00451 END IF 00452 * 00453 * 5) Reduce the bandwidth. 00454 * 00455 IF( KL.LT.N-1 ) THEN 00456 * 00457 * Reduce bandwidth -- kill column 00458 * 00459 DO 60 JCR = KL + 1, N - 1 00460 IC = JCR - KL 00461 IROWS = N + 1 - JCR 00462 ICOLS = N + KL - JCR 00463 * 00464 CALL CCOPY( IROWS, A( JCR, IC ), 1, WORK, 1 ) 00465 XNORMS = WORK( 1 ) 00466 CALL CLARFG( IROWS, XNORMS, WORK( 2 ), 1, TAU ) 00467 TAU = CONJG( TAU ) 00468 WORK( 1 ) = CONE 00469 ALPHA = CLARND( 5, ISEED ) 00470 * 00471 CALL CGEMV( 'C', IROWS, ICOLS, CONE, A( JCR, IC+1 ), LDA, 00472 $ WORK, 1, CZERO, WORK( IROWS+1 ), 1 ) 00473 CALL CGERC( IROWS, ICOLS, -TAU, WORK, 1, WORK( IROWS+1 ), 1, 00474 $ A( JCR, IC+1 ), LDA ) 00475 * 00476 CALL CGEMV( 'N', N, IROWS, CONE, A( 1, JCR ), LDA, WORK, 1, 00477 $ CZERO, WORK( IROWS+1 ), 1 ) 00478 CALL CGERC( N, IROWS, -CONJG( TAU ), WORK( IROWS+1 ), 1, 00479 $ WORK, 1, A( 1, JCR ), LDA ) 00480 * 00481 A( JCR, IC ) = XNORMS 00482 CALL CLASET( 'Full', IROWS-1, 1, CZERO, CZERO, 00483 $ A( JCR+1, IC ), LDA ) 00484 * 00485 CALL CSCAL( ICOLS+1, ALPHA, A( JCR, IC ), LDA ) 00486 CALL CSCAL( N, CONJG( ALPHA ), A( 1, JCR ), 1 ) 00487 60 CONTINUE 00488 ELSE IF( KU.LT.N-1 ) THEN 00489 * 00490 * Reduce upper bandwidth -- kill a row at a time. 00491 * 00492 DO 70 JCR = KU + 1, N - 1 00493 IR = JCR - KU 00494 IROWS = N + KU - JCR 00495 ICOLS = N + 1 - JCR 00496 * 00497 CALL CCOPY( ICOLS, A( IR, JCR ), LDA, WORK, 1 ) 00498 XNORMS = WORK( 1 ) 00499 CALL CLARFG( ICOLS, XNORMS, WORK( 2 ), 1, TAU ) 00500 TAU = CONJG( TAU ) 00501 WORK( 1 ) = CONE 00502 CALL CLACGV( ICOLS-1, WORK( 2 ), 1 ) 00503 ALPHA = CLARND( 5, ISEED ) 00504 * 00505 CALL CGEMV( 'N', IROWS, ICOLS, CONE, A( IR+1, JCR ), LDA, 00506 $ WORK, 1, CZERO, WORK( ICOLS+1 ), 1 ) 00507 CALL CGERC( IROWS, ICOLS, -TAU, WORK( ICOLS+1 ), 1, WORK, 1, 00508 $ A( IR+1, JCR ), LDA ) 00509 * 00510 CALL CGEMV( 'C', ICOLS, N, CONE, A( JCR, 1 ), LDA, WORK, 1, 00511 $ CZERO, WORK( ICOLS+1 ), 1 ) 00512 CALL CGERC( ICOLS, N, -CONJG( TAU ), WORK, 1, 00513 $ WORK( ICOLS+1 ), 1, A( JCR, 1 ), LDA ) 00514 * 00515 A( IR, JCR ) = XNORMS 00516 CALL CLASET( 'Full', 1, ICOLS-1, CZERO, CZERO, 00517 $ A( IR, JCR+1 ), LDA ) 00518 * 00519 CALL CSCAL( IROWS+1, ALPHA, A( IR, JCR ), 1 ) 00520 CALL CSCAL( N, CONJG( ALPHA ), A( JCR, 1 ), LDA ) 00521 70 CONTINUE 00522 END IF 00523 * 00524 * Scale the matrix to have norm ANORM 00525 * 00526 IF( ANORM.GE.ZERO ) THEN 00527 TEMP = CLANGE( 'M', N, N, A, LDA, TEMPA ) 00528 IF( TEMP.GT.ZERO ) THEN 00529 RALPHA = ANORM / TEMP 00530 DO 80 J = 1, N 00531 CALL CSSCAL( N, RALPHA, A( 1, J ), 1 ) 00532 80 CONTINUE 00533 END IF 00534 END IF 00535 * 00536 RETURN 00537 * 00538 * End of CLATME 00539 * 00540 END
1.7.3 
