![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SDRVSX 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 SDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00012 * NIUNIT, NOUNIT, A, LDA, H, HT, WR, WI, WRT, 00013 * WIT, WRTMP, WITMP, VS, LDVS, VS1, RESULT, WORK, 00014 * LWORK, IWORK, BWORK, INFO ) 00015 * 00016 * .. Scalar Arguments .. 00017 * INTEGER INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES, 00018 * $ NTYPES 00019 * REAL THRESH 00020 * .. 00021 * .. Array Arguments .. 00022 * LOGICAL BWORK( * ), DOTYPE( * ) 00023 * INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 00024 * REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ), 00025 * $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ), 00026 * $ WI( * ), WIT( * ), WITMP( * ), WORK( * ), 00027 * $ WR( * ), WRT( * ), WRTMP( * ) 00028 * .. 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> SDRVSX checks the nonsymmetric eigenvalue (Schur form) problem 00037 *> expert driver SGEESX. 00038 *> 00039 *> SDRVSX 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 SDRVSX 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 WR+sqrt(-1)*WI 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 WR+sqrt(-1)*WI are eigenvalues of T 00083 *> 1/ulp otherwise 00084 *> If workspace sufficient, also compare WR, WI 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 signs. 00121 *> (ULP = (first number larger than 1) - 1 ) 00122 *> (5) A diagonal matrix with geometrically spaced entries 00123 *> 1, ..., ULP and random signs. 00124 *> (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 00125 *> and random signs. 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 orthogonal and 00133 *> T has evenly spaced entries 1, ..., ULP with random signs 00134 *> on the diagonal and random O(1) entries in the upper 00135 *> triangle. 00136 *> 00137 *> (10) A matrix of the form U' T U, where U is orthogonal and 00138 *> T has geometrically spaced entries 1, ..., ULP with random 00139 *> signs on the diagonal and random O(1) entries in the upper 00140 *> 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 *> signs on the diagonal and random O(1) entries in the upper 00145 *> triangle. 00146 *> 00147 *> (12) A matrix of the form U' T U, where U is orthogonal and 00148 *> T has real or complex conjugate paired eigenvalues randomly 00149 *> chosen from ( ULP, 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 signs on the diagonal and random O(1) entries 00155 *> 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 signs on the diagonal and random 00160 *> 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 signs on the diagonal and random O(1) entries 00165 *> 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 real or complex conjugate paired 00169 *> eigenvalues randomly chosen from ( ULP, 1 ) and random 00170 *> O(1) entries in the upper 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 SGEESX 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 SGEESX 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 SDRVSX to continue the same random number 00266 *> sequence. 00267 *> \endverbatim 00268 *> 00269 *> \param[in] THRESH 00270 *> \verbatim 00271 *> THRESH is REAL 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 REAL 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 REAL array, dimension (LDA, max(NN)) 00311 *> Another copy of the test matrix A, modified by SGEESX. 00312 *> \endverbatim 00313 *> 00314 *> \param[out] HT 00315 *> \verbatim 00316 *> HT is REAL array, dimension (LDA, max(NN)) 00317 *> Yet another copy of the test matrix A, modified by SGEESX. 00318 *> \endverbatim 00319 *> 00320 *> \param[out] WR 00321 *> \verbatim 00322 *> WR is REAL array, dimension (max(NN)) 00323 *> \endverbatim 00324 *> 00325 *> \param[out] WI 00326 *> \verbatim 00327 *> WI is REAL array, dimension (max(NN)) 00328 *> 00329 *> The real and imaginary parts of the eigenvalues of A. 00330 *> On exit, WR + WI*i are the eigenvalues of the matrix in A. 00331 *> \endverbatim 00332 *> 00333 *> \param[out] WRT 00334 *> \verbatim 00335 *> WRT is REAL array, dimension (max(NN)) 00336 *> \endverbatim 00337 *> 00338 *> \param[out] WIT 00339 *> \verbatim 00340 *> WIT is REAL array, dimension (max(NN)) 00341 *> 00342 *> Like WR, WI, these arrays contain the eigenvalues of A, 00343 *> but those computed when SGEESX only computes a partial 00344 *> eigendecomposition, i.e. not Schur vectors 00345 *> \endverbatim 00346 *> 00347 *> \param[out] WRTMP 00348 *> \verbatim 00349 *> WRTMP is REAL array, dimension (max(NN)) 00350 *> \endverbatim 00351 *> 00352 *> \param[out] WITMP 00353 *> \verbatim 00354 *> WITMP is REAL array, dimension (max(NN)) 00355 *> 00356 *> More temporary storage for eigenvalues. 00357 *> \endverbatim 00358 *> 00359 *> \param[out] VS 00360 *> \verbatim 00361 *> VS is REAL array, dimension (LDVS, max(NN)) 00362 *> VS holds the computed Schur vectors. 00363 *> \endverbatim 00364 *> 00365 *> \param[in] LDVS 00366 *> \verbatim 00367 *> LDVS is INTEGER 00368 *> Leading dimension of VS. Must be at least max(1,max(NN)). 00369 *> \endverbatim 00370 *> 00371 *> \param[out] VS1 00372 *> \verbatim 00373 *> VS1 is REAL array, dimension (LDVS, max(NN)) 00374 *> VS1 holds another copy of the computed Schur vectors. 00375 *> \endverbatim 00376 *> 00377 *> \param[out] RESULT 00378 *> \verbatim 00379 *> RESULT is REAL array, dimension (17) 00380 *> The values computed by the 17 tests described above. 00381 *> The values are currently limited to 1/ulp, to avoid overflow. 00382 *> \endverbatim 00383 *> 00384 *> \param[out] WORK 00385 *> \verbatim 00386 *> WORK is REAL array, dimension (LWORK) 00387 *> \endverbatim 00388 *> 00389 *> \param[in] LWORK 00390 *> \verbatim 00391 *> LWORK is INTEGER 00392 *> The number of entries in WORK. This must be at least 00393 *> max(3*NN(j),2*NN(j)**2) for all j. 00394 *> \endverbatim 00395 *> 00396 *> \param[out] IWORK 00397 *> \verbatim 00398 *> IWORK is INTEGER array, dimension (max(NN)*max(NN)) 00399 *> \endverbatim 00400 *> 00401 *> \param[out] BWORK 00402 *> \verbatim 00403 *> BWORK is LOGICAL array, dimension (max(NN)) 00404 *> \endverbatim 00405 *> 00406 *> \param[out] INFO 00407 *> \verbatim 00408 *> INFO is INTEGER 00409 *> If 0, successful exit. 00410 *> <0, input parameter -INFO is incorrect 00411 *> >0, SLATMR, SLATMS, SLATME or SGET24 returned an error 00412 *> code and INFO is its absolute value 00413 *> 00414 *>----------------------------------------------------------------------- 00415 *> 00416 *> Some Local Variables and Parameters: 00417 *> ---- ----- --------- --- ---------- 00418 *> ZERO, ONE Real 0 and 1. 00419 *> MAXTYP The number of types defined. 00420 *> NMAX Largest value in NN. 00421 *> NERRS The number of tests which have exceeded THRESH 00422 *> COND, CONDS, 00423 *> IMODE Values to be passed to the matrix generators. 00424 *> ANORM Norm of A; passed to matrix generators. 00425 *> 00426 *> OVFL, UNFL Overflow and underflow thresholds. 00427 *> ULP, ULPINV Finest relative precision and its inverse. 00428 *> RTULP, RTULPI Square roots of the previous 4 values. 00429 *> The following four arrays decode JTYPE: 00430 *> KTYPE(j) The general type (1-10) for type "j". 00431 *> KMODE(j) The MODE value to be passed to the matrix 00432 *> generator for type "j". 00433 *> KMAGN(j) The order of magnitude ( O(1), 00434 *> O(overflow^(1/2) ), O(underflow^(1/2) ) 00435 *> KCONDS(j) Selectw whether CONDS is to be 1 or 00436 *> 1/sqrt(ulp). (0 means irrelevant.) 00437 *> \endverbatim 00438 * 00439 * Authors: 00440 * ======== 00441 * 00442 *> \author Univ. of Tennessee 00443 *> \author Univ. of California Berkeley 00444 *> \author Univ. of Colorado Denver 00445 *> \author NAG Ltd. 00446 * 00447 *> \date November 2011 00448 * 00449 *> \ingroup single_eig 00450 * 00451 * ===================================================================== 00452 SUBROUTINE SDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00453 $ NIUNIT, NOUNIT, A, LDA, H, HT, WR, WI, WRT, 00454 $ WIT, WRTMP, WITMP, VS, LDVS, VS1, RESULT, WORK, 00455 $ LWORK, IWORK, BWORK, INFO ) 00456 * 00457 * -- LAPACK test routine (version 3.4.0) -- 00458 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00459 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00460 * November 2011 00461 * 00462 * .. Scalar Arguments .. 00463 INTEGER INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES, 00464 $ NTYPES 00465 REAL THRESH 00466 * .. 00467 * .. Array Arguments .. 00468 LOGICAL BWORK( * ), DOTYPE( * ) 00469 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 00470 REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ), 00471 $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ), 00472 $ WI( * ), WIT( * ), WITMP( * ), WORK( * ), 00473 $ WR( * ), WRT( * ), WRTMP( * ) 00474 * .. 00475 * 00476 * ===================================================================== 00477 * 00478 * .. Parameters .. 00479 REAL ZERO, ONE 00480 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00481 INTEGER MAXTYP 00482 PARAMETER ( MAXTYP = 21 ) 00483 * .. 00484 * .. Local Scalars .. 00485 LOGICAL BADNN 00486 CHARACTER*3 PATH 00487 INTEGER I, IINFO, IMODE, ITYPE, IWK, J, JCOL, JSIZE, 00488 $ JTYPE, MTYPES, N, NERRS, NFAIL, NMAX, 00489 $ NNWORK, NSLCT, NTEST, NTESTF, NTESTT 00490 REAL ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN, 00491 $ RTULP, RTULPI, ULP, ULPINV, UNFL 00492 * .. 00493 * .. Local Arrays .. 00494 CHARACTER ADUMMA( 1 ) 00495 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ), 00496 $ KCONDS( MAXTYP ), KMAGN( MAXTYP ), 00497 $ KMODE( MAXTYP ), KTYPE( MAXTYP ) 00498 * .. 00499 * .. Arrays in Common .. 00500 LOGICAL SELVAL( 20 ) 00501 REAL SELWI( 20 ), SELWR( 20 ) 00502 * .. 00503 * .. Scalars in Common .. 00504 INTEGER SELDIM, SELOPT 00505 * .. 00506 * .. Common blocks .. 00507 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI 00508 * .. 00509 * .. External Functions .. 00510 REAL SLAMCH 00511 EXTERNAL SLAMCH 00512 * .. 00513 * .. External Subroutines .. 00514 EXTERNAL SGET24, SLABAD, SLASUM, SLATME, SLATMR, SLATMS, 00515 $ SLASET, XERBLA 00516 * .. 00517 * .. Intrinsic Functions .. 00518 INTRINSIC ABS, MAX, MIN, SQRT 00519 * .. 00520 * .. Data statements .. 00521 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 / 00522 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2, 00523 $ 3, 1, 2, 3 / 00524 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3, 00525 $ 1, 5, 5, 5, 4, 3, 1 / 00526 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 / 00527 * .. 00528 * .. Executable Statements .. 00529 * 00530 PATH( 1: 1 ) = 'Single precision' 00531 PATH( 2: 3 ) = 'SX' 00532 * 00533 * Check for errors 00534 * 00535 NTESTT = 0 00536 NTESTF = 0 00537 INFO = 0 00538 * 00539 * Important constants 00540 * 00541 BADNN = .FALSE. 00542 * 00543 * 12 is the largest dimension in the input file of precomputed 00544 * problems 00545 * 00546 NMAX = 12 00547 DO 10 J = 1, NSIZES 00548 NMAX = MAX( NMAX, NN( J ) ) 00549 IF( NN( J ).LT.0 ) 00550 $ BADNN = .TRUE. 00551 10 CONTINUE 00552 * 00553 * Check for errors 00554 * 00555 IF( NSIZES.LT.0 ) THEN 00556 INFO = -1 00557 ELSE IF( BADNN ) THEN 00558 INFO = -2 00559 ELSE IF( NTYPES.LT.0 ) THEN 00560 INFO = -3 00561 ELSE IF( THRESH.LT.ZERO ) THEN 00562 INFO = -6 00563 ELSE IF( NIUNIT.LE.0 ) THEN 00564 INFO = -7 00565 ELSE IF( NOUNIT.LE.0 ) THEN 00566 INFO = -8 00567 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN 00568 INFO = -10 00569 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.NMAX ) THEN 00570 INFO = -20 00571 ELSE IF( MAX( 3*NMAX, 2*NMAX**2 ).GT.LWORK ) THEN 00572 INFO = -24 00573 END IF 00574 * 00575 IF( INFO.NE.0 ) THEN 00576 CALL XERBLA( 'SDRVSX', -INFO ) 00577 RETURN 00578 END IF 00579 * 00580 * If nothing to do check on NIUNIT 00581 * 00582 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00583 $ GO TO 150 00584 * 00585 * More Important constants 00586 * 00587 UNFL = SLAMCH( 'Safe minimum' ) 00588 OVFL = ONE / UNFL 00589 CALL SLABAD( UNFL, OVFL ) 00590 ULP = SLAMCH( 'Precision' ) 00591 ULPINV = ONE / ULP 00592 RTULP = SQRT( ULP ) 00593 RTULPI = ONE / RTULP 00594 * 00595 * Loop over sizes, types 00596 * 00597 NERRS = 0 00598 * 00599 DO 140 JSIZE = 1, NSIZES 00600 N = NN( JSIZE ) 00601 IF( NSIZES.NE.1 ) THEN 00602 MTYPES = MIN( MAXTYP, NTYPES ) 00603 ELSE 00604 MTYPES = MIN( MAXTYP+1, NTYPES ) 00605 END IF 00606 * 00607 DO 130 JTYPE = 1, MTYPES 00608 IF( .NOT.DOTYPE( JTYPE ) ) 00609 $ GO TO 130 00610 * 00611 * Save ISEED in case of an error. 00612 * 00613 DO 20 J = 1, 4 00614 IOLDSD( J ) = ISEED( J ) 00615 20 CONTINUE 00616 * 00617 * Compute "A" 00618 * 00619 * Control parameters: 00620 * 00621 * KMAGN KCONDS KMODE KTYPE 00622 * =1 O(1) 1 clustered 1 zero 00623 * =2 large large clustered 2 identity 00624 * =3 small exponential Jordan 00625 * =4 arithmetic diagonal, (w/ eigenvalues) 00626 * =5 random log symmetric, w/ eigenvalues 00627 * =6 random general, w/ eigenvalues 00628 * =7 random diagonal 00629 * =8 random symmetric 00630 * =9 random general 00631 * =10 random triangular 00632 * 00633 IF( MTYPES.GT.MAXTYP ) 00634 $ GO TO 90 00635 * 00636 ITYPE = KTYPE( JTYPE ) 00637 IMODE = KMODE( JTYPE ) 00638 * 00639 * Compute norm 00640 * 00641 GO TO ( 30, 40, 50 )KMAGN( JTYPE ) 00642 * 00643 30 CONTINUE 00644 ANORM = ONE 00645 GO TO 60 00646 * 00647 40 CONTINUE 00648 ANORM = OVFL*ULP 00649 GO TO 60 00650 * 00651 50 CONTINUE 00652 ANORM = UNFL*ULPINV 00653 GO TO 60 00654 * 00655 60 CONTINUE 00656 * 00657 CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 00658 IINFO = 0 00659 COND = ULPINV 00660 * 00661 * Special Matrices -- Identity & Jordan block 00662 * 00663 * Zero 00664 * 00665 IF( ITYPE.EQ.1 ) THEN 00666 IINFO = 0 00667 * 00668 ELSE IF( ITYPE.EQ.2 ) THEN 00669 * 00670 * Identity 00671 * 00672 DO 70 JCOL = 1, N 00673 A( JCOL, JCOL ) = ANORM 00674 70 CONTINUE 00675 * 00676 ELSE IF( ITYPE.EQ.3 ) THEN 00677 * 00678 * Jordan Block 00679 * 00680 DO 80 JCOL = 1, N 00681 A( JCOL, JCOL ) = ANORM 00682 IF( JCOL.GT.1 ) 00683 $ A( JCOL, JCOL-1 ) = ONE 00684 80 CONTINUE 00685 * 00686 ELSE IF( ITYPE.EQ.4 ) THEN 00687 * 00688 * Diagonal Matrix, [Eigen]values Specified 00689 * 00690 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 00691 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 00692 $ IINFO ) 00693 * 00694 ELSE IF( ITYPE.EQ.5 ) THEN 00695 * 00696 * Symmetric, eigenvalues specified 00697 * 00698 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 00699 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 00700 $ IINFO ) 00701 * 00702 ELSE IF( ITYPE.EQ.6 ) THEN 00703 * 00704 * General, eigenvalues specified 00705 * 00706 IF( KCONDS( JTYPE ).EQ.1 ) THEN 00707 CONDS = ONE 00708 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN 00709 CONDS = RTULPI 00710 ELSE 00711 CONDS = ZERO 00712 END IF 00713 * 00714 ADUMMA( 1 ) = ' ' 00715 CALL SLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE, 00716 $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4, 00717 $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ), 00718 $ IINFO ) 00719 * 00720 ELSE IF( ITYPE.EQ.7 ) THEN 00721 * 00722 * Diagonal, random eigenvalues 00723 * 00724 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 00725 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00726 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 00727 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00728 * 00729 ELSE IF( ITYPE.EQ.8 ) THEN 00730 * 00731 * Symmetric, random eigenvalues 00732 * 00733 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 00734 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00735 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00736 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00737 * 00738 ELSE IF( ITYPE.EQ.9 ) THEN 00739 * 00740 * General, random eigenvalues 00741 * 00742 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 00743 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00744 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00745 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00746 IF( N.GE.4 ) THEN 00747 CALL SLASET( 'Full', 2, N, ZERO, ZERO, A, LDA ) 00748 CALL SLASET( 'Full', N-3, 1, ZERO, ZERO, A( 3, 1 ), 00749 $ LDA ) 00750 CALL SLASET( 'Full', N-3, 2, ZERO, ZERO, A( 3, N-1 ), 00751 $ LDA ) 00752 CALL SLASET( 'Full', 1, N, ZERO, ZERO, A( N, 1 ), 00753 $ LDA ) 00754 END IF 00755 * 00756 ELSE IF( ITYPE.EQ.10 ) THEN 00757 * 00758 * Triangular, random eigenvalues 00759 * 00760 CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 00761 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00762 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0, 00763 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00764 * 00765 ELSE 00766 * 00767 IINFO = 1 00768 END IF 00769 * 00770 IF( IINFO.NE.0 ) THEN 00771 WRITE( NOUNIT, FMT = 9991 )'Generator', IINFO, N, JTYPE, 00772 $ IOLDSD 00773 INFO = ABS( IINFO ) 00774 RETURN 00775 END IF 00776 * 00777 90 CONTINUE 00778 * 00779 * Test for minimal and generous workspace 00780 * 00781 DO 120 IWK = 1, 2 00782 IF( IWK.EQ.1 ) THEN 00783 NNWORK = 3*N 00784 ELSE 00785 NNWORK = MAX( 3*N, 2*N*N ) 00786 END IF 00787 NNWORK = MAX( NNWORK, 1 ) 00788 * 00789 CALL SGET24( .FALSE., JTYPE, THRESH, IOLDSD, NOUNIT, N, 00790 $ A, LDA, H, HT, WR, WI, WRT, WIT, WRTMP, 00791 $ WITMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT, 00792 $ ISLCT, RESULT, WORK, NNWORK, IWORK, BWORK, 00793 $ INFO ) 00794 * 00795 * Check for RESULT(j) > THRESH 00796 * 00797 NTEST = 0 00798 NFAIL = 0 00799 DO 100 J = 1, 15 00800 IF( RESULT( J ).GE.ZERO ) 00801 $ NTEST = NTEST + 1 00802 IF( RESULT( J ).GE.THRESH ) 00803 $ NFAIL = NFAIL + 1 00804 100 CONTINUE 00805 * 00806 IF( NFAIL.GT.0 ) 00807 $ NTESTF = NTESTF + 1 00808 IF( NTESTF.EQ.1 ) THEN 00809 WRITE( NOUNIT, FMT = 9999 )PATH 00810 WRITE( NOUNIT, FMT = 9998 ) 00811 WRITE( NOUNIT, FMT = 9997 ) 00812 WRITE( NOUNIT, FMT = 9996 ) 00813 WRITE( NOUNIT, FMT = 9995 )THRESH 00814 WRITE( NOUNIT, FMT = 9994 ) 00815 NTESTF = 2 00816 END IF 00817 * 00818 DO 110 J = 1, 15 00819 IF( RESULT( J ).GE.THRESH ) THEN 00820 WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE, 00821 $ J, RESULT( J ) 00822 END IF 00823 110 CONTINUE 00824 * 00825 NERRS = NERRS + NFAIL 00826 NTESTT = NTESTT + NTEST 00827 * 00828 120 CONTINUE 00829 130 CONTINUE 00830 140 CONTINUE 00831 * 00832 150 CONTINUE 00833 * 00834 * Read in data from file to check accuracy of condition estimation 00835 * Read input data until N=0 00836 * 00837 JTYPE = 0 00838 160 CONTINUE 00839 READ( NIUNIT, FMT = *, END = 200 )N, NSLCT 00840 IF( N.EQ.0 ) 00841 $ GO TO 200 00842 JTYPE = JTYPE + 1 00843 ISEED( 1 ) = JTYPE 00844 IF( NSLCT.GT.0 ) 00845 $ READ( NIUNIT, FMT = * )( ISLCT( I ), I = 1, NSLCT ) 00846 DO 170 I = 1, N 00847 READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N ) 00848 170 CONTINUE 00849 READ( NIUNIT, FMT = * )RCDEIN, RCDVIN 00850 * 00851 CALL SGET24( .TRUE., 22, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT, 00852 $ WR, WI, WRT, WIT, WRTMP, WITMP, VS, LDVS, VS1, 00853 $ RCDEIN, RCDVIN, NSLCT, ISLCT, RESULT, WORK, LWORK, 00854 $ IWORK, BWORK, INFO ) 00855 * 00856 * Check for RESULT(j) > THRESH 00857 * 00858 NTEST = 0 00859 NFAIL = 0 00860 DO 180 J = 1, 17 00861 IF( RESULT( J ).GE.ZERO ) 00862 $ NTEST = NTEST + 1 00863 IF( RESULT( J ).GE.THRESH ) 00864 $ NFAIL = NFAIL + 1 00865 180 CONTINUE 00866 * 00867 IF( NFAIL.GT.0 ) 00868 $ NTESTF = NTESTF + 1 00869 IF( NTESTF.EQ.1 ) THEN 00870 WRITE( NOUNIT, FMT = 9999 )PATH 00871 WRITE( NOUNIT, FMT = 9998 ) 00872 WRITE( NOUNIT, FMT = 9997 ) 00873 WRITE( NOUNIT, FMT = 9996 ) 00874 WRITE( NOUNIT, FMT = 9995 )THRESH 00875 WRITE( NOUNIT, FMT = 9994 ) 00876 NTESTF = 2 00877 END IF 00878 DO 190 J = 1, 17 00879 IF( RESULT( J ).GE.THRESH ) THEN 00880 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, J, RESULT( J ) 00881 END IF 00882 190 CONTINUE 00883 * 00884 NERRS = NERRS + NFAIL 00885 NTESTT = NTESTT + NTEST 00886 GO TO 160 00887 200 CONTINUE 00888 * 00889 * Summary 00890 * 00891 CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT ) 00892 * 00893 9999 FORMAT( / 1X, A3, ' -- Real Schur Form Decomposition Expert ', 00894 $ 'Driver', / ' Matrix types (see SDRVSX for details):' ) 00895 * 00896 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ', 00897 $ ' ', ' 5=Diagonal: geometr. spaced entries.', 00898 $ / ' 2=Identity matrix. ', ' 6=Diagona', 00899 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ', 00900 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ', 00901 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s', 00902 $ 'mall, evenly spaced.' ) 00903 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev', 00904 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e', 00905 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ', 00906 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond', 00907 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp', 00908 $ 'lex ', / ' 12=Well-cond., random complex ', ' ', 00909 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi', 00910 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.', 00911 $ ' complx ' ) 00912 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ', 00913 $ 'with small random entries.', / ' 20=Matrix with large ran', 00914 $ 'dom entries. ', / ) 00915 9995 FORMAT( ' Tests performed with test threshold =', F8.2, 00916 $ / ' ( A denotes A on input and T denotes A on output)', 00917 $ / / ' 1 = 0 if T in Schur form (no sort), ', 00918 $ ' 1/ulp otherwise', / 00919 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)', 00920 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', / 00921 $ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),', 00922 $ ' 1/ulp otherwise', / 00923 $ ' 5 = 0 if T same no matter if VS computed (no sort),', 00924 $ ' 1/ulp otherwise', / 00925 $ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)', 00926 $ ', 1/ulp otherwise' ) 00927 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise', 00928 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)', 00929 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ', 00930 $ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),', 00931 $ ' 1/ulp otherwise', / 00932 $ ' 11 = 0 if T same no matter what else computed (sort),', 00933 $ ' 1/ulp otherwise', / 00934 $ ' 12 = 0 if WR, WI same no matter what else computed ', 00935 $ '(sort), 1/ulp otherwise', / 00936 $ ' 13 = 0 if sorting succesful, 1/ulp otherwise', 00937 $ / ' 14 = 0 if RCONDE same no matter what else computed,', 00938 $ ' 1/ulp otherwise', / 00939 $ ' 15 = 0 if RCONDv same no matter what else computed,', 00940 $ ' 1/ulp otherwise', / 00941 $ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),', 00942 $ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' ) 00943 9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ), 00944 $ ' type ', I2, ', test(', I2, ')=', G10.3 ) 00945 9992 FORMAT( ' N=', I5, ', input example =', I3, ', test(', I2, ')=', 00946 $ G10.3 ) 00947 9991 FORMAT( ' SDRVSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 00948 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 00949 * 00950 RETURN 00951 * 00952 * End of SDRVSX 00953 * 00954 END