![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SDRVEV 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 SDRVEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00012 * NOUNIT, A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, 00013 * VR, LDVR, LRE, LDLRE, RESULT, WORK, NWORK, 00014 * IWORK, INFO ) 00015 * 00016 * .. Scalar Arguments .. 00017 * INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NOUNIT, NSIZES, 00018 * $ NTYPES, NWORK 00019 * REAL THRESH 00020 * .. 00021 * .. Array Arguments .. 00022 * LOGICAL DOTYPE( * ) 00023 * INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 00024 * REAL A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ), 00025 * $ RESULT( 7 ), VL( LDVL, * ), VR( LDVR, * ), 00026 * $ WI( * ), WI1( * ), WORK( * ), WR( * ), WR1( * ) 00027 * .. 00028 * 00029 * 00030 *> \par Purpose: 00031 * ============= 00032 *> 00033 *> \verbatim 00034 *> 00035 *> SDRVEV checks the nonsymmetric eigenvalue problem driver SGEEV. 00036 *> 00037 *> When SDRVEV is called, a number of matrix "sizes" ("n's") and a 00038 *> number of matrix "types" are specified. For each size ("n") 00039 *> and each type of matrix, one matrix will be generated and used 00040 *> to test the nonsymmetric eigenroutines. For each matrix, 7 00041 *> tests will be performed: 00042 *> 00043 *> (1) | A * VR - VR * W | / ( n |A| ulp ) 00044 *> 00045 *> Here VR is the matrix of unit right eigenvectors. 00046 *> W is a block diagonal matrix, with a 1x1 block for each 00047 *> real eigenvalue and a 2x2 block for each complex conjugate 00048 *> pair. If eigenvalues j and j+1 are a complex conjugate pair, 00049 *> so WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 00050 *> 2 x 2 block corresponding to the pair will be: 00051 *> 00052 *> ( wr wi ) 00053 *> ( -wi wr ) 00054 *> 00055 *> Such a block multiplying an n x 2 matrix ( ur ui ) on the 00056 *> right will be the same as multiplying ur + i*ui by wr + i*wi. 00057 *> 00058 *> (2) | A**H * VL - VL * W**H | / ( n |A| ulp ) 00059 *> 00060 *> Here VL is the matrix of unit left eigenvectors, A**H is the 00061 *> conjugate transpose of A, and W is as above. 00062 *> 00063 *> (3) | |VR(i)| - 1 | / ulp and whether largest component real 00064 *> 00065 *> VR(i) denotes the i-th column of VR. 00066 *> 00067 *> (4) | |VL(i)| - 1 | / ulp and whether largest component real 00068 *> 00069 *> VL(i) denotes the i-th column of VL. 00070 *> 00071 *> (5) W(full) = W(partial) 00072 *> 00073 *> W(full) denotes the eigenvalues computed when both VR and VL 00074 *> are also computed, and W(partial) denotes the eigenvalues 00075 *> computed when only W, only W and VR, or only W and VL are 00076 *> computed. 00077 *> 00078 *> (6) VR(full) = VR(partial) 00079 *> 00080 *> VR(full) denotes the right eigenvectors computed when both VR 00081 *> and VL are computed, and VR(partial) denotes the result 00082 *> when only VR is computed. 00083 *> 00084 *> (7) VL(full) = VL(partial) 00085 *> 00086 *> VL(full) denotes the left eigenvectors computed when both VR 00087 *> and VL are also computed, and VL(partial) denotes the result 00088 *> when only VL is computed. 00089 *> 00090 *> The "sizes" are specified by an array NN(1:NSIZES); the value of 00091 *> each element NN(j) specifies one size. 00092 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); 00093 *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 00094 *> Currently, the list of possible types is: 00095 *> 00096 *> (1) The zero matrix. 00097 *> (2) The identity matrix. 00098 *> (3) A (transposed) Jordan block, with 1's on the diagonal. 00099 *> 00100 *> (4) A diagonal matrix with evenly spaced entries 00101 *> 1, ..., ULP and random signs. 00102 *> (ULP = (first number larger than 1) - 1 ) 00103 *> (5) A diagonal matrix with geometrically spaced entries 00104 *> 1, ..., ULP and random signs. 00105 *> (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 00106 *> and random signs. 00107 *> 00108 *> (7) Same as (4), but multiplied by a constant near 00109 *> the overflow threshold 00110 *> (8) Same as (4), but multiplied by a constant near 00111 *> the underflow threshold 00112 *> 00113 *> (9) A matrix of the form U' T U, where U is orthogonal and 00114 *> T has evenly spaced entries 1, ..., ULP with random signs 00115 *> on the diagonal and random O(1) entries in the upper 00116 *> triangle. 00117 *> 00118 *> (10) A matrix of the form U' T U, where U is orthogonal and 00119 *> T has geometrically spaced entries 1, ..., ULP with random 00120 *> signs on the diagonal and random O(1) entries in the upper 00121 *> triangle. 00122 *> 00123 *> (11) A matrix of the form U' T U, where U is orthogonal and 00124 *> T has "clustered" entries 1, ULP,..., ULP with random 00125 *> signs on the diagonal and random O(1) entries in the upper 00126 *> triangle. 00127 *> 00128 *> (12) A matrix of the form U' T U, where U is orthogonal and 00129 *> T has real or complex conjugate paired eigenvalues randomly 00130 *> chosen from ( ULP, 1 ) and random O(1) entries in the upper 00131 *> triangle. 00132 *> 00133 *> (13) A matrix of the form X' T X, where X has condition 00134 *> SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP 00135 *> with random signs on the diagonal and random O(1) entries 00136 *> in the upper triangle. 00137 *> 00138 *> (14) A matrix of the form X' T X, where X has condition 00139 *> SQRT( ULP ) and T has geometrically spaced entries 00140 *> 1, ..., ULP with random signs on the diagonal and random 00141 *> O(1) entries in the upper triangle. 00142 *> 00143 *> (15) A matrix of the form X' T X, where X has condition 00144 *> SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP 00145 *> with random signs on the diagonal and random O(1) entries 00146 *> in the upper triangle. 00147 *> 00148 *> (16) A matrix of the form X' T X, where X has condition 00149 *> SQRT( ULP ) and T has real or complex conjugate paired 00150 *> eigenvalues randomly chosen from ( ULP, 1 ) and random 00151 *> O(1) entries in the upper triangle. 00152 *> 00153 *> (17) Same as (16), but multiplied by a constant 00154 *> near the overflow threshold 00155 *> (18) Same as (16), but multiplied by a constant 00156 *> near the underflow threshold 00157 *> 00158 *> (19) Nonsymmetric matrix with random entries chosen from (-1,1). 00159 *> If N is at least 4, all entries in first two rows and last 00160 *> row, and first column and last two columns are zero. 00161 *> (20) Same as (19), but multiplied by a constant 00162 *> near the overflow threshold 00163 *> (21) Same as (19), but multiplied by a constant 00164 *> near the underflow threshold 00165 *> \endverbatim 00166 * 00167 * Arguments: 00168 * ========== 00169 * 00170 *> \param[in] NSIZES 00171 *> \verbatim 00172 *> NSIZES is INTEGER 00173 *> The number of sizes of matrices to use. If it is zero, 00174 *> SDRVEV does nothing. It must be at least zero. 00175 *> \endverbatim 00176 *> 00177 *> \param[in] NN 00178 *> \verbatim 00179 *> NN is INTEGER array, dimension (NSIZES) 00180 *> An array containing the sizes to be used for the matrices. 00181 *> Zero values will be skipped. The values must be at least 00182 *> zero. 00183 *> \endverbatim 00184 *> 00185 *> \param[in] NTYPES 00186 *> \verbatim 00187 *> NTYPES is INTEGER 00188 *> The number of elements in DOTYPE. If it is zero, SDRVEV 00189 *> does nothing. It must be at least zero. If it is MAXTYP+1 00190 *> and NSIZES is 1, then an additional type, MAXTYP+1 is 00191 *> defined, which is to use whatever matrix is in A. This 00192 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 00193 *> DOTYPE(MAXTYP+1) is .TRUE. . 00194 *> \endverbatim 00195 *> 00196 *> \param[in] DOTYPE 00197 *> \verbatim 00198 *> DOTYPE is LOGICAL array, dimension (NTYPES) 00199 *> If DOTYPE(j) is .TRUE., then for each size in NN a 00200 *> matrix of that size and of type j will be generated. 00201 *> If NTYPES is smaller than the maximum number of types 00202 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through 00203 *> MAXTYP will not be generated. If NTYPES is larger 00204 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 00205 *> will be ignored. 00206 *> \endverbatim 00207 *> 00208 *> \param[in,out] ISEED 00209 *> \verbatim 00210 *> ISEED is INTEGER array, dimension (4) 00211 *> On entry ISEED specifies the seed of the random number 00212 *> generator. The array elements should be between 0 and 4095; 00213 *> if not they will be reduced mod 4096. Also, ISEED(4) must 00214 *> be odd. The random number generator uses a linear 00215 *> congruential sequence limited to small integers, and so 00216 *> should produce machine independent random numbers. The 00217 *> values of ISEED are changed on exit, and can be used in the 00218 *> next call to SDRVEV to continue the same random number 00219 *> sequence. 00220 *> \endverbatim 00221 *> 00222 *> \param[in] THRESH 00223 *> \verbatim 00224 *> THRESH is REAL 00225 *> A test will count as "failed" if the "error", computed as 00226 *> described above, exceeds THRESH. Note that the error 00227 *> is scaled to be O(1), so THRESH should be a reasonably 00228 *> small multiple of 1, e.g., 10 or 100. In particular, 00229 *> it should not depend on the precision (single vs. double) 00230 *> or the size of the matrix. It must be at least zero. 00231 *> \endverbatim 00232 *> 00233 *> \param[in] NOUNIT 00234 *> \verbatim 00235 *> NOUNIT is INTEGER 00236 *> The FORTRAN unit number for printing out error messages 00237 *> (e.g., if a routine returns INFO not equal to 0.) 00238 *> \endverbatim 00239 *> 00240 *> \param[out] A 00241 *> \verbatim 00242 *> A is REAL array, dimension (LDA, max(NN)) 00243 *> Used to hold the matrix whose eigenvalues are to be 00244 *> computed. On exit, A contains the last matrix actually used. 00245 *> \endverbatim 00246 *> 00247 *> \param[in] LDA 00248 *> \verbatim 00249 *> LDA is INTEGER 00250 *> The leading dimension of A, and H. LDA must be at 00251 *> least 1 and at least max(NN). 00252 *> \endverbatim 00253 *> 00254 *> \param[out] H 00255 *> \verbatim 00256 *> H is REAL array, dimension (LDA, max(NN)) 00257 *> Another copy of the test matrix A, modified by SGEEV. 00258 *> \endverbatim 00259 *> 00260 *> \param[out] WR 00261 *> \verbatim 00262 *> WR is REAL array, dimension (max(NN)) 00263 *> \endverbatim 00264 *> 00265 *> \param[out] WI 00266 *> \verbatim 00267 *> WI is REAL array, dimension (max(NN)) 00268 *> 00269 *> The real and imaginary parts of the eigenvalues of A. 00270 *> On exit, WR + WI*i are the eigenvalues of the matrix in A. 00271 *> \endverbatim 00272 *> 00273 *> \param[out] WR1 00274 *> \verbatim 00275 *> WR1 is REAL array, dimension (max(NN)) 00276 *> \endverbatim 00277 *> 00278 *> \param[out] WI1 00279 *> \verbatim 00280 *> WI1 is REAL array, dimension (max(NN)) 00281 *> 00282 *> Like WR, WI, these arrays contain the eigenvalues of A, 00283 *> but those computed when SGEEV only computes a partial 00284 *> eigendecomposition, i.e. not the eigenvalues and left 00285 *> and right eigenvectors. 00286 *> \endverbatim 00287 *> 00288 *> \param[out] VL 00289 *> \verbatim 00290 *> VL is REAL array, dimension (LDVL, max(NN)) 00291 *> VL holds the computed left eigenvectors. 00292 *> \endverbatim 00293 *> 00294 *> \param[in] LDVL 00295 *> \verbatim 00296 *> LDVL is INTEGER 00297 *> Leading dimension of VL. Must be at least max(1,max(NN)). 00298 *> \endverbatim 00299 *> 00300 *> \param[out] VR 00301 *> \verbatim 00302 *> VR is REAL array, dimension (LDVR, max(NN)) 00303 *> VR holds the computed right eigenvectors. 00304 *> \endverbatim 00305 *> 00306 *> \param[in] LDVR 00307 *> \verbatim 00308 *> LDVR is INTEGER 00309 *> Leading dimension of VR. Must be at least max(1,max(NN)). 00310 *> \endverbatim 00311 *> 00312 *> \param[out] LRE 00313 *> \verbatim 00314 *> LRE is REAL array, dimension (LDLRE,max(NN)) 00315 *> LRE holds the computed right or left eigenvectors. 00316 *> \endverbatim 00317 *> 00318 *> \param[in] LDLRE 00319 *> \verbatim 00320 *> LDLRE is INTEGER 00321 *> Leading dimension of LRE. Must be at least max(1,max(NN)). 00322 *> \endverbatim 00323 *> 00324 *> \param[out] RESULT 00325 *> \verbatim 00326 *> RESULT is REAL array, dimension (7) 00327 *> The values computed by the seven tests described above. 00328 *> The values are currently limited to 1/ulp, to avoid overflow. 00329 *> \endverbatim 00330 *> 00331 *> \param[out] WORK 00332 *> \verbatim 00333 *> WORK is REAL array, dimension (NWORK) 00334 *> \endverbatim 00335 *> 00336 *> \param[in] NWORK 00337 *> \verbatim 00338 *> NWORK is INTEGER 00339 *> The number of entries in WORK. This must be at least 00340 *> 5*NN(j)+2*NN(j)**2 for all j. 00341 *> \endverbatim 00342 *> 00343 *> \param[out] IWORK 00344 *> \verbatim 00345 *> IWORK is INTEGER array, dimension (max(NN)) 00346 *> \endverbatim 00347 *> 00348 *> \param[out] INFO 00349 *> \verbatim 00350 *> INFO is INTEGER 00351 *> If 0, then everything ran OK. 00352 *> -1: NSIZES < 0 00353 *> -2: Some NN(j) < 0 00354 *> -3: NTYPES < 0 00355 *> -6: THRESH < 0 00356 *> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). 00357 *> -16: LDVL < 1 or LDVL < NMAX, where NMAX is max( NN(j) ). 00358 *> -18: LDVR < 1 or LDVR < NMAX, where NMAX is max( NN(j) ). 00359 *> -20: LDLRE < 1 or LDLRE < NMAX, where NMAX is max( NN(j) ). 00360 *> -23: NWORK too small. 00361 *> If SLATMR, SLATMS, SLATME or SGEEV returns an error code, 00362 *> the absolute value of it is returned. 00363 *> 00364 *>----------------------------------------------------------------------- 00365 *> 00366 *> Some Local Variables and Parameters: 00367 *> ---- ----- --------- --- ---------- 00368 *> 00369 *> ZERO, ONE Real 0 and 1. 00370 *> MAXTYP The number of types defined. 00371 *> NMAX Largest value in NN. 00372 *> NERRS The number of tests which have exceeded THRESH 00373 *> COND, CONDS, 00374 *> IMODE Values to be passed to the matrix generators. 00375 *> ANORM Norm of A; passed to matrix generators. 00376 *> 00377 *> OVFL, UNFL Overflow and underflow thresholds. 00378 *> ULP, ULPINV Finest relative precision and its inverse. 00379 *> RTULP, RTULPI Square roots of the previous 4 values. 00380 *> 00381 *> The following four arrays decode JTYPE: 00382 *> KTYPE(j) The general type (1-10) for type "j". 00383 *> KMODE(j) The MODE value to be passed to the matrix 00384 *> generator for type "j". 00385 *> KMAGN(j) The order of magnitude ( O(1), 00386 *> O(overflow^(1/2) ), O(underflow^(1/2) ) 00387 *> KCONDS(j) Selectw whether CONDS is to be 1 or 00388 *> 1/sqrt(ulp). (0 means irrelevant.) 00389 *> \endverbatim 00390 * 00391 * Authors: 00392 * ======== 00393 * 00394 *> \author Univ. of Tennessee 00395 *> \author Univ. of California Berkeley 00396 *> \author Univ. of Colorado Denver 00397 *> \author NAG Ltd. 00398 * 00399 *> \date November 2011 00400 * 00401 *> \ingroup single_eig 00402 * 00403 * ===================================================================== 00404 SUBROUTINE SDRVEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00405 $ NOUNIT, A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, 00406 $ VR, LDVR, LRE, LDLRE, RESULT, WORK, NWORK, 00407 $ IWORK, INFO ) 00408 * 00409 * -- LAPACK test routine (version 3.4.0) -- 00410 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00411 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00412 * November 2011 00413 * 00414 * .. Scalar Arguments .. 00415 INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NOUNIT, NSIZES, 00416 $ NTYPES, NWORK 00417 REAL THRESH 00418 * .. 00419 * .. Array Arguments .. 00420 LOGICAL DOTYPE( * ) 00421 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 00422 REAL A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ), 00423 $ RESULT( 7 ), VL( LDVL, * ), VR( LDVR, * ), 00424 $ WI( * ), WI1( * ), WORK( * ), WR( * ), WR1( * ) 00425 * .. 00426 * 00427 * ===================================================================== 00428 * 00429 * .. Parameters .. 00430 REAL ZERO, ONE 00431 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00432 REAL TWO 00433 PARAMETER ( TWO = 2.0E0 ) 00434 INTEGER MAXTYP 00435 PARAMETER ( MAXTYP = 21 ) 00436 * .. 00437 * .. Local Scalars .. 00438 LOGICAL BADNN 00439 CHARACTER*3 PATH 00440 INTEGER IINFO, IMODE, ITYPE, IWK, J, JCOL, JJ, JSIZE, 00441 $ JTYPE, MTYPES, N, NERRS, NFAIL, NMAX, 00442 $ NNWORK, NTEST, NTESTF, NTESTT 00443 REAL ANORM, COND, CONDS, OVFL, RTULP, RTULPI, TNRM, 00444 $ ULP, ULPINV, UNFL, VMX, VRMX, VTST 00445 * .. 00446 * .. Local Arrays .. 00447 CHARACTER ADUMMA( 1 ) 00448 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ), 00449 $ KMAGN( MAXTYP ), KMODE( MAXTYP ), 00450 $ KTYPE( MAXTYP ) 00451 REAL DUM( 1 ), RES( 2 ) 00452 * .. 00453 * .. External Functions .. 00454 REAL SLAMCH, SLAPY2, SNRM2 00455 EXTERNAL SLAMCH, SLAPY2, SNRM2 00456 * .. 00457 * .. External Subroutines .. 00458 EXTERNAL SGEEV, SGET22, SLABAD, SLACPY, SLASUM, SLATME, 00459 $ SLATMR, SLATMS, SLASET, XERBLA 00460 * .. 00461 * .. Intrinsic Functions .. 00462 INTRINSIC ABS, MAX, MIN, SQRT 00463 * .. 00464 * .. Data statements .. 00465 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 / 00466 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2, 00467 $ 3, 1, 2, 3 / 00468 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3, 00469 $ 1, 5, 5, 5, 4, 3, 1 / 00470 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 / 00471 * .. 00472 * .. Executable Statements .. 00473 * 00474 PATH( 1: 1 ) = 'Single precision' 00475 PATH( 2: 3 ) = 'EV' 00476 * 00477 * Check for errors 00478 * 00479 NTESTT = 0 00480 NTESTF = 0 00481 INFO = 0 00482 * 00483 * Important constants 00484 * 00485 BADNN = .FALSE. 00486 NMAX = 0 00487 DO 10 J = 1, NSIZES 00488 NMAX = MAX( NMAX, NN( J ) ) 00489 IF( NN( J ).LT.0 ) 00490 $ BADNN = .TRUE. 00491 10 CONTINUE 00492 * 00493 * Check for errors 00494 * 00495 IF( NSIZES.LT.0 ) THEN 00496 INFO = -1 00497 ELSE IF( BADNN ) THEN 00498 INFO = -2 00499 ELSE IF( NTYPES.LT.0 ) THEN 00500 INFO = -3 00501 ELSE IF( THRESH.LT.ZERO ) THEN 00502 INFO = -6 00503 ELSE IF( NOUNIT.LE.0 ) THEN 00504 INFO = -7 00505 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN 00506 INFO = -9 00507 ELSE IF( LDVL.LT.1 .OR. LDVL.LT.NMAX ) THEN 00508 INFO = -16 00509 ELSE IF( LDVR.LT.1 .OR. LDVR.LT.NMAX ) THEN 00510 INFO = -18 00511 ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.NMAX ) THEN 00512 INFO = -20 00513 ELSE IF( 5*NMAX+2*NMAX**2.GT.NWORK ) THEN 00514 INFO = -23 00515 END IF 00516 * 00517 IF( INFO.NE.0 ) THEN 00518 CALL XERBLA( 'SDRVEV', -INFO ) 00519 RETURN 00520 END IF 00521 * 00522 * Quick return if nothing to do 00523 * 00524 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00525 $ RETURN 00526 * 00527 * More Important constants 00528 * 00529 UNFL = SLAMCH( 'Safe minimum' ) 00530 OVFL = ONE / UNFL 00531 CALL SLABAD( UNFL, OVFL ) 00532 ULP = SLAMCH( 'Precision' ) 00533 ULPINV = ONE / ULP 00534 RTULP = SQRT( ULP ) 00535 RTULPI = ONE / RTULP 00536 * 00537 * Loop over sizes, types 00538 * 00539 NERRS = 0 00540 * 00541 DO 270 JSIZE = 1, NSIZES 00542 N = NN( JSIZE ) 00543 IF( NSIZES.NE.1 ) THEN 00544 MTYPES = MIN( MAXTYP, NTYPES ) 00545 ELSE 00546 MTYPES = MIN( MAXTYP+1, NTYPES ) 00547 END IF 00548 * 00549 DO 260 JTYPE = 1, MTYPES 00550 IF( .NOT.DOTYPE( JTYPE ) ) 00551 $ GO TO 260 00552 * 00553 * Save ISEED in case of an error. 00554 * 00555 DO 20 J = 1, 4 00556 IOLDSD( J ) = ISEED( J ) 00557 20 CONTINUE 00558 * 00559 * Compute "A" 00560 * 00561 * Control parameters: 00562 * 00563 * KMAGN KCONDS KMODE KTYPE 00564 * =1 O(1) 1 clustered 1 zero 00565 * =2 large large clustered 2 identity 00566 * =3 small exponential Jordan 00567 * =4 arithmetic diagonal, (w/ eigenvalues) 00568 * =5 random log symmetric, w/ eigenvalues 00569 * =6 random general, w/ eigenvalues 00570 * =7 random diagonal 00571 * =8 random symmetric 00572 * =9 random general 00573 * =10 random triangular 00574 * 00575 IF( MTYPES.GT.MAXTYP ) 00576 $ GO TO 90 00577 * 00578 ITYPE = KTYPE( JTYPE ) 00579 IMODE = KMODE( JTYPE ) 00580 * 00581 * Compute norm 00582 * 00583 GO TO ( 30, 40, 50 )KMAGN( JTYPE ) 00584 * 00585 30 CONTINUE 00586 ANORM = ONE 00587 GO TO 60 00588 * 00589 40 CONTINUE 00590 ANORM = OVFL*ULP 00591 GO TO 60 00592 * 00593 50 CONTINUE 00594 ANORM = UNFL*ULPINV 00595 GO TO 60 00596 * 00597 60 CONTINUE 00598 * 00599 CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 00600 IINFO = 0 00601 COND = ULPINV 00602 * 00603 * Special Matrices -- Identity & Jordan block 00604 * 00605 * Zero 00606 * 00607 IF( ITYPE.EQ.1 ) THEN 00608 IINFO = 0 00609 * 00610 ELSE IF( ITYPE.EQ.2 ) THEN 00611 * 00612 * Identity 00613 * 00614 DO 70 JCOL = 1, N 00615 A( JCOL, JCOL ) = ANORM 00616 70 CONTINUE 00617 * 00618 ELSE IF( ITYPE.EQ.3 ) THEN 00619 * 00620 * Jordan Block 00621 * 00622 DO 80 JCOL = 1, N 00623 A( JCOL, JCOL ) = ANORM 00624 IF( JCOL.GT.1 ) 00625 $ A( JCOL, JCOL-1 ) = ONE 00626 80 CONTINUE 00627 * 00628 ELSE IF( ITYPE.EQ.4 ) THEN 00629 * 00630 * Diagonal Matrix, [Eigen]values Specified 00631 * 00632 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 00633 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 00634 $ IINFO ) 00635 * 00636 ELSE IF( ITYPE.EQ.5 ) THEN 00637 * 00638 * Symmetric, eigenvalues specified 00639 * 00640 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 00641 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 00642 $ IINFO ) 00643 * 00644 ELSE IF( ITYPE.EQ.6 ) THEN 00645 * 00646 * General, eigenvalues specified 00647 * 00648 IF( KCONDS( JTYPE ).EQ.1 ) THEN 00649 CONDS = ONE 00650 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 00651 CONDS = RTULPI 00652 ELSE 00653 CONDS = ZERO 00654 END IF 00655 * 00656 ADUMMA( 1 ) = ' ' 00657 CALL SLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE, 00658 $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4, 00659 $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ), 00660 $ IINFO ) 00661 * 00662 ELSE IF( ITYPE.EQ.7 ) THEN 00663 * 00664 * Diagonal, random eigenvalues 00665 * 00666 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 00667 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00668 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 00669 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00670 * 00671 ELSE IF( ITYPE.EQ.8 ) THEN 00672 * 00673 * Symmetric, random eigenvalues 00674 * 00675 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 00676 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00677 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00678 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00679 * 00680 ELSE IF( ITYPE.EQ.9 ) THEN 00681 * 00682 * General, random eigenvalues 00683 * 00684 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 00685 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00686 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00687 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00688 IF( N.GE.4 ) THEN 00689 CALL SLASET( 'Full', 2, N, ZERO, ZERO, A, LDA ) 00690 CALL SLASET( 'Full', N-3, 1, ZERO, ZERO, A( 3, 1 ), 00691 $ LDA ) 00692 CALL SLASET( 'Full', N-3, 2, ZERO, ZERO, A( 3, N-1 ), 00693 $ LDA ) 00694 CALL SLASET( 'Full', 1, N, ZERO, ZERO, A( N, 1 ), 00695 $ LDA ) 00696 END IF 00697 * 00698 ELSE IF( ITYPE.EQ.10 ) THEN 00699 * 00700 * Triangular, random eigenvalues 00701 * 00702 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 00703 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00704 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 00705 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00706 * 00707 ELSE 00708 * 00709 IINFO = 1 00710 END IF 00711 * 00712 IF( IINFO.NE.0 ) THEN 00713 WRITE( NOUNIT, FMT = 9993 )'Generator', IINFO, N, JTYPE, 00714 $ IOLDSD 00715 INFO = ABS( IINFO ) 00716 RETURN 00717 END IF 00718 * 00719 90 CONTINUE 00720 * 00721 * Test for minimal and generous workspace 00722 * 00723 DO 250 IWK = 1, 2 00724 IF( IWK.EQ.1 ) THEN 00725 NNWORK = 4*N 00726 ELSE 00727 NNWORK = 5*N + 2*N**2 00728 END IF 00729 NNWORK = MAX( NNWORK, 1 ) 00730 * 00731 * Initialize RESULT 00732 * 00733 DO 100 J = 1, 7 00734 RESULT( J ) = -ONE 00735 100 CONTINUE 00736 * 00737 * Compute eigenvalues and eigenvectors, and test them 00738 * 00739 CALL SLACPY( 'F', N, N, A, LDA, H, LDA ) 00740 CALL SGEEV( 'V', 'V', N, H, LDA, WR, WI, VL, LDVL, VR, 00741 $ LDVR, WORK, NNWORK, IINFO ) 00742 IF( IINFO.NE.0 ) THEN 00743 RESULT( 1 ) = ULPINV 00744 WRITE( NOUNIT, FMT = 9993 )'SGEEV1', IINFO, N, JTYPE, 00745 $ IOLDSD 00746 INFO = ABS( IINFO ) 00747 GO TO 220 00748 END IF 00749 * 00750 * Do Test (1) 00751 * 00752 CALL SGET22( 'N', 'N', 'N', N, A, LDA, VR, LDVR, WR, WI, 00753 $ WORK, RES ) 00754 RESULT( 1 ) = RES( 1 ) 00755 * 00756 * Do Test (2) 00757 * 00758 CALL SGET22( 'T', 'N', 'T', N, A, LDA, VL, LDVL, WR, WI, 00759 $ WORK, RES ) 00760 RESULT( 2 ) = RES( 1 ) 00761 * 00762 * Do Test (3) 00763 * 00764 DO 120 J = 1, N 00765 TNRM = ONE 00766 IF( WI( J ).EQ.ZERO ) THEN 00767 TNRM = SNRM2( N, VR( 1, J ), 1 ) 00768 ELSE IF( WI( J ).GT.ZERO ) THEN 00769 TNRM = SLAPY2( SNRM2( N, VR( 1, J ), 1 ), 00770 $ SNRM2( N, VR( 1, J+1 ), 1 ) ) 00771 END IF 00772 RESULT( 3 ) = MAX( RESULT( 3 ), 00773 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) ) 00774 IF( WI( J ).GT.ZERO ) THEN 00775 VMX = ZERO 00776 VRMX = ZERO 00777 DO 110 JJ = 1, N 00778 VTST = SLAPY2( VR( JJ, J ), VR( JJ, J+1 ) ) 00779 IF( VTST.GT.VMX ) 00780 $ VMX = VTST 00781 IF( VR( JJ, J+1 ).EQ.ZERO .AND. 00782 $ ABS( VR( JJ, J ) ).GT.VRMX ) 00783 $ VRMX = ABS( VR( JJ, J ) ) 00784 110 CONTINUE 00785 IF( VRMX / VMX.LT.ONE-TWO*ULP ) 00786 $ RESULT( 3 ) = ULPINV 00787 END IF 00788 120 CONTINUE 00789 * 00790 * Do Test (4) 00791 * 00792 DO 140 J = 1, N 00793 TNRM = ONE 00794 IF( WI( J ).EQ.ZERO ) THEN 00795 TNRM = SNRM2( N, VL( 1, J ), 1 ) 00796 ELSE IF( WI( J ).GT.ZERO ) THEN 00797 TNRM = SLAPY2( SNRM2( N, VL( 1, J ), 1 ), 00798 $ SNRM2( N, VL( 1, J+1 ), 1 ) ) 00799 END IF 00800 RESULT( 4 ) = MAX( RESULT( 4 ), 00801 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) ) 00802 IF( WI( J ).GT.ZERO ) THEN 00803 VMX = ZERO 00804 VRMX = ZERO 00805 DO 130 JJ = 1, N 00806 VTST = SLAPY2( VL( JJ, J ), VL( JJ, J+1 ) ) 00807 IF( VTST.GT.VMX ) 00808 $ VMX = VTST 00809 IF( VL( JJ, J+1 ).EQ.ZERO .AND. 00810 $ ABS( VL( JJ, J ) ).GT.VRMX ) 00811 $ VRMX = ABS( VL( JJ, J ) ) 00812 130 CONTINUE 00813 IF( VRMX / VMX.LT.ONE-TWO*ULP ) 00814 $ RESULT( 4 ) = ULPINV 00815 END IF 00816 140 CONTINUE 00817 * 00818 * Compute eigenvalues only, and test them 00819 * 00820 CALL SLACPY( 'F', N, N, A, LDA, H, LDA ) 00821 CALL SGEEV( 'N', 'N', N, H, LDA, WR1, WI1, DUM, 1, DUM, 00822 $ 1, WORK, NNWORK, IINFO ) 00823 IF( IINFO.NE.0 ) THEN 00824 RESULT( 1 ) = ULPINV 00825 WRITE( NOUNIT, FMT = 9993 )'SGEEV2', IINFO, N, JTYPE, 00826 $ IOLDSD 00827 INFO = ABS( IINFO ) 00828 GO TO 220 00829 END IF 00830 * 00831 * Do Test (5) 00832 * 00833 DO 150 J = 1, N 00834 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 00835 $ RESULT( 5 ) = ULPINV 00836 150 CONTINUE 00837 * 00838 * Compute eigenvalues and right eigenvectors, and test them 00839 * 00840 CALL SLACPY( 'F', N, N, A, LDA, H, LDA ) 00841 CALL SGEEV( 'N', 'V', N, H, LDA, WR1, WI1, DUM, 1, LRE, 00842 $ LDLRE, WORK, NNWORK, IINFO ) 00843 IF( IINFO.NE.0 ) THEN 00844 RESULT( 1 ) = ULPINV 00845 WRITE( NOUNIT, FMT = 9993 )'SGEEV3', IINFO, N, JTYPE, 00846 $ IOLDSD 00847 INFO = ABS( IINFO ) 00848 GO TO 220 00849 END IF 00850 * 00851 * Do Test (5) again 00852 * 00853 DO 160 J = 1, N 00854 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 00855 $ RESULT( 5 ) = ULPINV 00856 160 CONTINUE 00857 * 00858 * Do Test (6) 00859 * 00860 DO 180 J = 1, N 00861 DO 170 JJ = 1, N 00862 IF( VR( J, JJ ).NE.LRE( J, JJ ) ) 00863 $ RESULT( 6 ) = ULPINV 00864 170 CONTINUE 00865 180 CONTINUE 00866 * 00867 * Compute eigenvalues and left eigenvectors, and test them 00868 * 00869 CALL SLACPY( 'F', N, N, A, LDA, H, LDA ) 00870 CALL SGEEV( 'V', 'N', N, H, LDA, WR1, WI1, LRE, LDLRE, 00871 $ DUM, 1, WORK, NNWORK, IINFO ) 00872 IF( IINFO.NE.0 ) THEN 00873 RESULT( 1 ) = ULPINV 00874 WRITE( NOUNIT, FMT = 9993 )'SGEEV4', IINFO, N, JTYPE, 00875 $ IOLDSD 00876 INFO = ABS( IINFO ) 00877 GO TO 220 00878 END IF 00879 * 00880 * Do Test (5) again 00881 * 00882 DO 190 J = 1, N 00883 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) ) 00884 $ RESULT( 5 ) = ULPINV 00885 190 CONTINUE 00886 * 00887 * Do Test (7) 00888 * 00889 DO 210 J = 1, N 00890 DO 200 JJ = 1, N 00891 IF( VL( J, JJ ).NE.LRE( J, JJ ) ) 00892 $ RESULT( 7 ) = ULPINV 00893 200 CONTINUE 00894 210 CONTINUE 00895 * 00896 * End of Loop -- Check for RESULT(j) > THRESH 00897 * 00898 220 CONTINUE 00899 * 00900 NTEST = 0 00901 NFAIL = 0 00902 DO 230 J = 1, 7 00903 IF( RESULT( J ).GE.ZERO ) 00904 $ NTEST = NTEST + 1 00905 IF( RESULT( J ).GE.THRESH ) 00906 $ NFAIL = NFAIL + 1 00907 230 CONTINUE 00908 * 00909 IF( NFAIL.GT.0 ) 00910 $ NTESTF = NTESTF + 1 00911 IF( NTESTF.EQ.1 ) THEN 00912 WRITE( NOUNIT, FMT = 9999 )PATH 00913 WRITE( NOUNIT, FMT = 9998 ) 00914 WRITE( NOUNIT, FMT = 9997 ) 00915 WRITE( NOUNIT, FMT = 9996 ) 00916 WRITE( NOUNIT, FMT = 9995 )THRESH 00917 NTESTF = 2 00918 END IF 00919 * 00920 DO 240 J = 1, 7 00921 IF( RESULT( J ).GE.THRESH ) THEN 00922 WRITE( NOUNIT, FMT = 9994 )N, IWK, IOLDSD, JTYPE, 00923 $ J, RESULT( J ) 00924 END IF 00925 240 CONTINUE 00926 * 00927 NERRS = NERRS + NFAIL 00928 NTESTT = NTESTT + NTEST 00929 * 00930 250 CONTINUE 00931 260 CONTINUE 00932 270 CONTINUE 00933 * 00934 * Summary 00935 * 00936 CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT ) 00937 * 00938 9999 FORMAT( / 1X, A3, ' -- Real Eigenvalue-Eigenvector Decomposition', 00939 $ ' Driver', / ' Matrix types (see SDRVEV for details): ' ) 00940 * 00941 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ', 00942 $ ' ', ' 5=Diagonal: geometr. spaced entries.', 00943 $ / ' 2=Identity matrix. ', ' 6=Diagona', 00944 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ', 00945 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ', 00946 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s', 00947 $ 'mall, evenly spaced.' ) 00948 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev', 00949 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e', 00950 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ', 00951 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond', 00952 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp', 00953 $ 'lex ', / ' 12=Well-cond., random complex ', 6X, ' ', 00954 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi', 00955 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.', 00956 $ ' complx ' ) 00957 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ', 00958 $ 'with small random entries.', / ' 20=Matrix with large ran', 00959 $ 'dom entries. ', / ) 00960 9995 FORMAT( ' Tests performed with test threshold =', F8.2, 00961 $ / / ' 1 = | A VR - VR W | / ( n |A| ulp ) ', 00962 $ / ' 2 = | transpose(A) VL - VL W | / ( n |A| ulp ) ', 00963 $ / ' 3 = | |VR(i)| - 1 | / ulp ', 00964 $ / ' 4 = | |VL(i)| - 1 | / ulp ', 00965 $ / ' 5 = 0 if W same no matter if VR or VL computed,', 00966 $ ' 1/ulp otherwise', / 00967 $ ' 6 = 0 if VR same no matter if VL computed,', 00968 $ ' 1/ulp otherwise', / 00969 $ ' 7 = 0 if VL same no matter if VR computed,', 00970 $ ' 1/ulp otherwise', / ) 00971 9994 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ), 00972 $ ' type ', I2, ', test(', I2, ')=', G10.3 ) 00973 9993 FORMAT( ' SDRVEV: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00974 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 00975 * 00976 RETURN 00977 * 00978 * End of SDRVEV 00979 * 00980 END