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