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