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