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