![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SCHKST 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 SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00012 * NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5, 00013 * WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK, 00014 * LWORK, IWORK, LIWORK, RESULT, INFO ) 00015 * 00016 * .. Scalar Arguments .. 00017 * INTEGER INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES, 00018 * $ NTYPES 00019 * REAL THRESH 00020 * .. 00021 * .. Array Arguments .. 00022 * LOGICAL DOTYPE( * ) 00023 * INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 00024 * REAL A( LDA, * ), AP( * ), D1( * ), D2( * ), 00025 * $ D3( * ), D4( * ), D5( * ), RESULT( * ), 00026 * $ SD( * ), SE( * ), TAU( * ), U( LDU, * ), 00027 * $ V( LDU, * ), VP( * ), WA1( * ), WA2( * ), 00028 * $ WA3( * ), WORK( * ), WR( * ), Z( LDU, * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SCHKST checks the symmetric eigenvalue problem routines. 00038 *> 00039 *> SSYTRD factors A as U S U' , where ' means transpose, 00040 *> S is symmetric tridiagonal, and U is orthogonal. 00041 *> SSYTRD can use either just the lower or just the upper triangle 00042 *> of A; SCHKST checks both cases. 00043 *> U is represented as a product of Householder 00044 *> transformations, whose vectors are stored in the first 00045 *> n-1 columns of V, and whose scale factors are in TAU. 00046 *> 00047 *> SSPTRD does the same as SSYTRD, except that A and V are stored 00048 *> in "packed" format. 00049 *> 00050 *> SORGTR constructs the matrix U from the contents of V and TAU. 00051 *> 00052 *> SOPGTR constructs the matrix U from the contents of VP and TAU. 00053 *> 00054 *> SSTEQR factors S as Z D1 Z' , where Z is the orthogonal 00055 *> matrix of eigenvectors and D1 is a diagonal matrix with 00056 *> the eigenvalues on the diagonal. D2 is the matrix of 00057 *> eigenvalues computed when Z is not computed. 00058 *> 00059 *> SSTERF computes D3, the matrix of eigenvalues, by the 00060 *> PWK method, which does not yield eigenvectors. 00061 *> 00062 *> SPTEQR factors S as Z4 D4 Z4' , for a 00063 *> symmetric positive definite tridiagonal matrix. 00064 *> D5 is the matrix of eigenvalues computed when Z is not 00065 *> computed. 00066 *> 00067 *> SSTEBZ computes selected eigenvalues. WA1, WA2, and 00068 *> WA3 will denote eigenvalues computed to high 00069 *> absolute accuracy, with different range options. 00070 *> WR will denote eigenvalues computed to high relative 00071 *> accuracy. 00072 *> 00073 *> SSTEIN computes Y, the eigenvectors of S, given the 00074 *> eigenvalues. 00075 *> 00076 *> SSTEDC factors S as Z D1 Z' , where Z is the orthogonal 00077 *> matrix of eigenvectors and D1 is a diagonal matrix with 00078 *> the eigenvalues on the diagonal ('I' option). It may also 00079 *> update an input orthogonal matrix, usually the output 00080 *> from SSYTRD/SORGTR or SSPTRD/SOPGTR ('V' option). It may 00081 *> also just compute eigenvalues ('N' option). 00082 *> 00083 *> SSTEMR factors S as Z D1 Z' , where Z is the orthogonal 00084 *> matrix of eigenvectors and D1 is a diagonal matrix with 00085 *> the eigenvalues on the diagonal ('I' option). SSTEMR 00086 *> uses the Relatively Robust Representation whenever possible. 00087 *> 00088 *> When SCHKST is called, a number of matrix "sizes" ("n's") and a 00089 *> number of matrix "types" are specified. For each size ("n") 00090 *> and each type of matrix, one matrix will be generated and used 00091 *> to test the symmetric eigenroutines. For each matrix, a number 00092 *> of tests will be performed: 00093 *> 00094 *> (1) | A - V S V' | / ( |A| n ulp ) SSYTRD( UPLO='U', ... ) 00095 *> 00096 *> (2) | I - UV' | / ( n ulp ) SORGTR( UPLO='U', ... ) 00097 *> 00098 *> (3) | A - V S V' | / ( |A| n ulp ) SSYTRD( UPLO='L', ... ) 00099 *> 00100 *> (4) | I - UV' | / ( n ulp ) SORGTR( UPLO='L', ... ) 00101 *> 00102 *> (5-8) Same as 1-4, but for SSPTRD and SOPGTR. 00103 *> 00104 *> (9) | S - Z D Z' | / ( |S| n ulp ) SSTEQR('V',...) 00105 *> 00106 *> (10) | I - ZZ' | / ( n ulp ) SSTEQR('V',...) 00107 *> 00108 *> (11) | D1 - D2 | / ( |D1| ulp ) SSTEQR('N',...) 00109 *> 00110 *> (12) | D1 - D3 | / ( |D1| ulp ) SSTERF 00111 *> 00112 *> (13) 0 if the true eigenvalues (computed by sturm count) 00113 *> of S are within THRESH of 00114 *> those in D1. 2*THRESH if they are not. (Tested using 00115 *> SSTECH) 00116 *> 00117 *> For S positive definite, 00118 *> 00119 *> (14) | S - Z4 D4 Z4' | / ( |S| n ulp ) SPTEQR('V',...) 00120 *> 00121 *> (15) | I - Z4 Z4' | / ( n ulp ) SPTEQR('V',...) 00122 *> 00123 *> (16) | D4 - D5 | / ( 100 |D4| ulp ) SPTEQR('N',...) 00124 *> 00125 *> When S is also diagonally dominant by the factor gamma < 1, 00126 *> 00127 *> (17) max | D4(i) - WR(i) | / ( |D4(i)| omega ) , 00128 *> i 00129 *> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 00130 *> SSTEBZ( 'A', 'E', ...) 00131 *> 00132 *> (18) | WA1 - D3 | / ( |D3| ulp ) SSTEBZ( 'A', 'E', ...) 00133 *> 00134 *> (19) ( max { min | WA2(i)-WA3(j) | } + 00135 *> i j 00136 *> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) 00137 *> i j 00138 *> SSTEBZ( 'I', 'E', ...) 00139 *> 00140 *> (20) | S - Y WA1 Y' | / ( |S| n ulp ) SSTEBZ, SSTEIN 00141 *> 00142 *> (21) | I - Y Y' | / ( n ulp ) SSTEBZ, SSTEIN 00143 *> 00144 *> (22) | S - Z D Z' | / ( |S| n ulp ) SSTEDC('I') 00145 *> 00146 *> (23) | I - ZZ' | / ( n ulp ) SSTEDC('I') 00147 *> 00148 *> (24) | S - Z D Z' | / ( |S| n ulp ) SSTEDC('V') 00149 *> 00150 *> (25) | I - ZZ' | / ( n ulp ) SSTEDC('V') 00151 *> 00152 *> (26) | D1 - D2 | / ( |D1| ulp ) SSTEDC('V') and 00153 *> SSTEDC('N') 00154 *> 00155 *> Test 27 is disabled at the moment because SSTEMR does not 00156 *> guarantee high relatvie accuracy. 00157 *> 00158 *> (27) max | D6(i) - WR(i) | / ( |D6(i)| omega ) , 00159 *> i 00160 *> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 00161 *> SSTEMR('V', 'A') 00162 *> 00163 *> (28) max | D6(i) - WR(i) | / ( |D6(i)| omega ) , 00164 *> i 00165 *> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 00166 *> SSTEMR('V', 'I') 00167 *> 00168 *> Tests 29 through 34 are disable at present because SSTEMR 00169 *> does not handle partial specturm requests. 00170 *> 00171 *> (29) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'I') 00172 *> 00173 *> (30) | I - ZZ' | / ( n ulp ) SSTEMR('V', 'I') 00174 *> 00175 *> (31) ( max { min | WA2(i)-WA3(j) | } + 00176 *> i j 00177 *> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) 00178 *> i j 00179 *> SSTEMR('N', 'I') vs. SSTEMR('V', 'I') 00180 *> 00181 *> (32) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'V') 00182 *> 00183 *> (33) | I - ZZ' | / ( n ulp ) SSTEMR('V', 'V') 00184 *> 00185 *> (34) ( max { min | WA2(i)-WA3(j) | } + 00186 *> i j 00187 *> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) 00188 *> i j 00189 *> SSTEMR('N', 'V') vs. SSTEMR('V', 'V') 00190 *> 00191 *> (35) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'A') 00192 *> 00193 *> (36) | I - ZZ' | / ( n ulp ) SSTEMR('V', 'A') 00194 *> 00195 *> (37) ( max { min | WA2(i)-WA3(j) | } + 00196 *> i j 00197 *> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) 00198 *> i j 00199 *> SSTEMR('N', 'A') vs. SSTEMR('V', 'A') 00200 *> 00201 *> The "sizes" are specified by an array NN(1:NSIZES); the value of 00202 *> each element NN(j) specifies one size. 00203 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); 00204 *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 00205 *> Currently, the list of possible types is: 00206 *> 00207 *> (1) The zero matrix. 00208 *> (2) The identity matrix. 00209 *> 00210 *> (3) A diagonal matrix with evenly spaced entries 00211 *> 1, ..., ULP and random signs. 00212 *> (ULP = (first number larger than 1) - 1 ) 00213 *> (4) A diagonal matrix with geometrically spaced entries 00214 *> 1, ..., ULP and random signs. 00215 *> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 00216 *> and random signs. 00217 *> 00218 *> (6) Same as (4), but multiplied by SQRT( overflow threshold ) 00219 *> (7) Same as (4), but multiplied by SQRT( underflow threshold ) 00220 *> 00221 *> (8) A matrix of the form U' D U, where U is orthogonal and 00222 *> D has evenly spaced entries 1, ..., ULP with random signs 00223 *> on the diagonal. 00224 *> 00225 *> (9) A matrix of the form U' D U, where U is orthogonal and 00226 *> D has geometrically spaced entries 1, ..., ULP with random 00227 *> signs on the diagonal. 00228 *> 00229 *> (10) A matrix of the form U' D U, where U is orthogonal and 00230 *> D has "clustered" entries 1, ULP,..., ULP with random 00231 *> signs on the diagonal. 00232 *> 00233 *> (11) Same as (8), but multiplied by SQRT( overflow threshold ) 00234 *> (12) Same as (8), but multiplied by SQRT( underflow threshold ) 00235 *> 00236 *> (13) Symmetric matrix with random entries chosen from (-1,1). 00237 *> (14) Same as (13), but multiplied by SQRT( overflow threshold ) 00238 *> (15) Same as (13), but multiplied by SQRT( underflow threshold ) 00239 *> (16) Same as (8), but diagonal elements are all positive. 00240 *> (17) Same as (9), but diagonal elements are all positive. 00241 *> (18) Same as (10), but diagonal elements are all positive. 00242 *> (19) Same as (16), but multiplied by SQRT( overflow threshold ) 00243 *> (20) Same as (16), but multiplied by SQRT( underflow threshold ) 00244 *> (21) A diagonally dominant tridiagonal matrix with geometrically 00245 *> spaced diagonal entries 1, ..., ULP. 00246 *> \endverbatim 00247 * 00248 * Arguments: 00249 * ========== 00250 * 00251 *> \param[in] NSIZES 00252 *> \verbatim 00253 *> NSIZES is INTEGER 00254 *> The number of sizes of matrices to use. If it is zero, 00255 *> SCHKST does nothing. It must be at least zero. 00256 *> \endverbatim 00257 *> 00258 *> \param[in] NN 00259 *> \verbatim 00260 *> NN is INTEGER array, dimension (NSIZES) 00261 *> An array containing the sizes to be used for the matrices. 00262 *> Zero values will be skipped. The values must be at least 00263 *> zero. 00264 *> \endverbatim 00265 *> 00266 *> \param[in] NTYPES 00267 *> \verbatim 00268 *> NTYPES is INTEGER 00269 *> The number of elements in DOTYPE. If it is zero, SCHKST 00270 *> does nothing. It must be at least zero. If it is MAXTYP+1 00271 *> and NSIZES is 1, then an additional type, MAXTYP+1 is 00272 *> defined, which is to use whatever matrix is in A. This 00273 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 00274 *> DOTYPE(MAXTYP+1) is .TRUE. . 00275 *> \endverbatim 00276 *> 00277 *> \param[in] DOTYPE 00278 *> \verbatim 00279 *> DOTYPE is LOGICAL array, dimension (NTYPES) 00280 *> If DOTYPE(j) is .TRUE., then for each size in NN a 00281 *> matrix of that size and of type j will be generated. 00282 *> If NTYPES is smaller than the maximum number of types 00283 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through 00284 *> MAXTYP will not be generated. If NTYPES is larger 00285 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 00286 *> will be ignored. 00287 *> \endverbatim 00288 *> 00289 *> \param[in,out] ISEED 00290 *> \verbatim 00291 *> ISEED is INTEGER array, dimension (4) 00292 *> On entry ISEED specifies the seed of the random number 00293 *> generator. The array elements should be between 0 and 4095; 00294 *> if not they will be reduced mod 4096. Also, ISEED(4) must 00295 *> be odd. The random number generator uses a linear 00296 *> congruential sequence limited to small integers, and so 00297 *> should produce machine independent random numbers. The 00298 *> values of ISEED are changed on exit, and can be used in the 00299 *> next call to SCHKST to continue the same random number 00300 *> sequence. 00301 *> \endverbatim 00302 *> 00303 *> \param[in] THRESH 00304 *> \verbatim 00305 *> THRESH is REAL 00306 *> A test will count as "failed" if the "error", computed as 00307 *> described above, exceeds THRESH. Note that the error 00308 *> is scaled to be O(1), so THRESH should be a reasonably 00309 *> small multiple of 1, e.g., 10 or 100. In particular, 00310 *> it should not depend on the precision (single vs. double) 00311 *> or the size of the matrix. It must be at least zero. 00312 *> \endverbatim 00313 *> 00314 *> \param[in] NOUNIT 00315 *> \verbatim 00316 *> NOUNIT is INTEGER 00317 *> The FORTRAN unit number for printing out error messages 00318 *> (e.g., if a routine returns IINFO not equal to 0.) 00319 *> \endverbatim 00320 *> 00321 *> \param[in,out] A 00322 *> \verbatim 00323 *> A is REAL array of 00324 *> dimension ( LDA , max(NN) ) 00325 *> Used to hold the matrix whose eigenvalues are to be 00326 *> computed. On exit, A contains the last matrix actually 00327 *> used. 00328 *> \endverbatim 00329 *> 00330 *> \param[in] LDA 00331 *> \verbatim 00332 *> LDA is INTEGER 00333 *> The leading dimension of A. It must be at 00334 *> least 1 and at least max( NN ). 00335 *> \endverbatim 00336 *> 00337 *> \param[out] AP 00338 *> \verbatim 00339 *> AP is REAL array of 00340 *> dimension( max(NN)*max(NN+1)/2 ) 00341 *> The matrix A stored in packed format. 00342 *> \endverbatim 00343 *> 00344 *> \param[out] SD 00345 *> \verbatim 00346 *> SD is REAL array of 00347 *> dimension( max(NN) ) 00348 *> The diagonal of the tridiagonal matrix computed by SSYTRD. 00349 *> On exit, SD and SE contain the tridiagonal form of the 00350 *> matrix in A. 00351 *> \endverbatim 00352 *> 00353 *> \param[out] SE 00354 *> \verbatim 00355 *> SE is REAL array of 00356 *> dimension( max(NN) ) 00357 *> The off-diagonal of the tridiagonal matrix computed by 00358 *> SSYTRD. On exit, SD and SE contain the tridiagonal form of 00359 *> the matrix in A. 00360 *> \endverbatim 00361 *> 00362 *> \param[out] D1 00363 *> \verbatim 00364 *> D1 is REAL array of 00365 *> dimension( max(NN) ) 00366 *> The eigenvalues of A, as computed by SSTEQR simlutaneously 00367 *> with Z. On exit, the eigenvalues in D1 correspond with the 00368 *> matrix in A. 00369 *> \endverbatim 00370 *> 00371 *> \param[out] D2 00372 *> \verbatim 00373 *> D2 is REAL array of 00374 *> dimension( max(NN) ) 00375 *> The eigenvalues of A, as computed by SSTEQR if Z is not 00376 *> computed. On exit, the eigenvalues in D2 correspond with 00377 *> the matrix in A. 00378 *> \endverbatim 00379 *> 00380 *> \param[out] D3 00381 *> \verbatim 00382 *> D3 is REAL array of 00383 *> dimension( max(NN) ) 00384 *> The eigenvalues of A, as computed by SSTERF. On exit, the 00385 *> eigenvalues in D3 correspond with the matrix in A. 00386 *> \endverbatim 00387 *> 00388 *> \param[out] D4 00389 *> \verbatim 00390 *> D4 is REAL array of 00391 *> dimension( max(NN) ) 00392 *> The eigenvalues of A, as computed by SPTEQR(V). 00393 *> ZPTEQR factors S as Z4 D4 Z4* 00394 *> On exit, the eigenvalues in D4 correspond with the matrix in A. 00395 *> \endverbatim 00396 *> 00397 *> \param[out] D5 00398 *> \verbatim 00399 *> D5 is REAL array of 00400 *> dimension( max(NN) ) 00401 *> The eigenvalues of A, as computed by SPTEQR(N) 00402 *> when Z is not computed. On exit, the 00403 *> eigenvalues in D4 correspond with the matrix in A. 00404 *> \endverbatim 00405 *> 00406 *> \param[out] WA1 00407 *> \verbatim 00408 *> WA1 is REAL array of 00409 *> dimension( max(NN) ) 00410 *> All eigenvalues of A, computed to high 00411 *> absolute accuracy, with different range options. 00412 *> as computed by SSTEBZ. 00413 *> \endverbatim 00414 *> 00415 *> \param[out] WA2 00416 *> \verbatim 00417 *> WA2 is REAL array of 00418 *> dimension( max(NN) ) 00419 *> Selected eigenvalues of A, computed to high 00420 *> absolute accuracy, with different range options. 00421 *> as computed by SSTEBZ. 00422 *> Choose random values for IL and IU, and ask for the 00423 *> IL-th through IU-th eigenvalues. 00424 *> \endverbatim 00425 *> 00426 *> \param[out] WA3 00427 *> \verbatim 00428 *> WA3 is REAL array of 00429 *> dimension( max(NN) ) 00430 *> Selected eigenvalues of A, computed to high 00431 *> absolute accuracy, with different range options. 00432 *> as computed by SSTEBZ. 00433 *> Determine the values VL and VU of the IL-th and IU-th 00434 *> eigenvalues and ask for all eigenvalues in this range. 00435 *> \endverbatim 00436 *> 00437 *> \param[out] WR 00438 *> \verbatim 00439 *> WR is REAL array of 00440 *> dimension( max(NN) ) 00441 *> All eigenvalues of A, computed to high 00442 *> absolute accuracy, with different options. 00443 *> as computed by SSTEBZ. 00444 *> \endverbatim 00445 *> 00446 *> \param[out] U 00447 *> \verbatim 00448 *> U is REAL array of 00449 *> dimension( LDU, max(NN) ). 00450 *> The orthogonal matrix computed by SSYTRD + SORGTR. 00451 *> \endverbatim 00452 *> 00453 *> \param[in] LDU 00454 *> \verbatim 00455 *> LDU is INTEGER 00456 *> The leading dimension of U, Z, and V. It must be at least 1 00457 *> and at least max( NN ). 00458 *> \endverbatim 00459 *> 00460 *> \param[out] V 00461 *> \verbatim 00462 *> V is REAL array of 00463 *> dimension( LDU, max(NN) ). 00464 *> The Housholder vectors computed by SSYTRD in reducing A to 00465 *> tridiagonal form. The vectors computed with UPLO='U' are 00466 *> in the upper triangle, and the vectors computed with UPLO='L' 00467 *> are in the lower triangle. (As described in SSYTRD, the 00468 *> sub- and superdiagonal are not set to 1, although the 00469 *> true Householder vector has a 1 in that position. The 00470 *> routines that use V, such as SORGTR, set those entries to 00471 *> 1 before using them, and then restore them later.) 00472 *> \endverbatim 00473 *> 00474 *> \param[out] VP 00475 *> \verbatim 00476 *> VP is REAL array of 00477 *> dimension( max(NN)*max(NN+1)/2 ) 00478 *> The matrix V stored in packed format. 00479 *> \endverbatim 00480 *> 00481 *> \param[out] TAU 00482 *> \verbatim 00483 *> TAU is REAL array of 00484 *> dimension( max(NN) ) 00485 *> The Householder factors computed by SSYTRD in reducing A 00486 *> to tridiagonal form. 00487 *> \endverbatim 00488 *> 00489 *> \param[out] Z 00490 *> \verbatim 00491 *> Z is REAL array of 00492 *> dimension( LDU, max(NN) ). 00493 *> The orthogonal matrix of eigenvectors computed by SSTEQR, 00494 *> SPTEQR, and SSTEIN. 00495 *> \endverbatim 00496 *> 00497 *> \param[out] WORK 00498 *> \verbatim 00499 *> WORK is REAL array of 00500 *> dimension( LWORK ) 00501 *> \endverbatim 00502 *> 00503 *> \param[in] LWORK 00504 *> \verbatim 00505 *> LWORK is INTEGER 00506 *> The number of entries in WORK. This must be at least 00507 *> 1 + 4 * Nmax + 2 * Nmax * lg Nmax + 3 * Nmax**2 00508 *> where Nmax = max( NN(j), 2 ) and lg = log base 2. 00509 *> \endverbatim 00510 *> 00511 *> \param[out] IWORK 00512 *> \verbatim 00513 *> IWORK is INTEGER array, 00514 *> Workspace. 00515 *> \endverbatim 00516 *> 00517 *> \param[out] LIWORK 00518 *> \verbatim 00519 *> LIWORK is INTEGER 00520 *> The number of entries in IWORK. This must be at least 00521 *> 6 + 6*Nmax + 5 * Nmax * lg Nmax 00522 *> where Nmax = max( NN(j), 2 ) and lg = log base 2. 00523 *> \endverbatim 00524 *> 00525 *> \param[out] RESULT 00526 *> \verbatim 00527 *> RESULT is REAL array, dimension (26) 00528 *> The values computed by the tests described above. 00529 *> The values are currently limited to 1/ulp, to avoid 00530 *> overflow. 00531 *> \endverbatim 00532 *> 00533 *> \param[out] INFO 00534 *> \verbatim 00535 *> INFO is INTEGER 00536 *> If 0, then everything ran OK. 00537 *> -1: NSIZES < 0 00538 *> -2: Some NN(j) < 0 00539 *> -3: NTYPES < 0 00540 *> -5: THRESH < 0 00541 *> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). 00542 *> -23: LDU < 1 or LDU < NMAX. 00543 *> -29: LWORK too small. 00544 *> If SLATMR, SLATMS, SSYTRD, SORGTR, SSTEQR, SSTERF, 00545 *> or SORMC2 returns an error code, the 00546 *> absolute value of it is returned. 00547 *> 00548 *>----------------------------------------------------------------------- 00549 *> 00550 *> Some Local Variables and Parameters: 00551 *> ---- ----- --------- --- ---------- 00552 *> ZERO, ONE Real 0 and 1. 00553 *> MAXTYP The number of types defined. 00554 *> NTEST The number of tests performed, or which can 00555 *> be performed so far, for the current matrix. 00556 *> NTESTT The total number of tests performed so far. 00557 *> NBLOCK Blocksize as returned by ENVIR. 00558 *> NMAX Largest value in NN. 00559 *> NMATS The number of matrices generated so far. 00560 *> NERRS The number of tests which have exceeded THRESH 00561 *> so far. 00562 *> COND, IMODE Values to be passed to the matrix generators. 00563 *> ANORM Norm of A; passed to matrix generators. 00564 *> 00565 *> OVFL, UNFL Overflow and underflow thresholds. 00566 *> ULP, ULPINV Finest relative precision and its inverse. 00567 *> RTOVFL, RTUNFL Square roots of the previous 2 values. 00568 *> The following four arrays decode JTYPE: 00569 *> KTYPE(j) The general type (1-10) for type "j". 00570 *> KMODE(j) The MODE value to be passed to the matrix 00571 *> generator for type "j". 00572 *> KMAGN(j) The order of magnitude ( O(1), 00573 *> O(overflow^(1/2) ), O(underflow^(1/2) ) 00574 *> \endverbatim 00575 * 00576 * Authors: 00577 * ======== 00578 * 00579 *> \author Univ. of Tennessee 00580 *> \author Univ. of California Berkeley 00581 *> \author Univ. of Colorado Denver 00582 *> \author NAG Ltd. 00583 * 00584 *> \date November 2011 00585 * 00586 *> \ingroup single_eig 00587 * 00588 * ===================================================================== 00589 SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 00590 $ NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5, 00591 $ WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK, 00592 $ LWORK, IWORK, LIWORK, RESULT, INFO ) 00593 * 00594 * -- LAPACK test routine (version 3.4.0) -- 00595 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00596 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00597 * November 2011 00598 * 00599 * .. Scalar Arguments .. 00600 INTEGER INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES, 00601 $ NTYPES 00602 REAL THRESH 00603 * .. 00604 * .. Array Arguments .. 00605 LOGICAL DOTYPE( * ) 00606 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 00607 REAL A( LDA, * ), AP( * ), D1( * ), D2( * ), 00608 $ D3( * ), D4( * ), D5( * ), RESULT( * ), 00609 $ SD( * ), SE( * ), TAU( * ), U( LDU, * ), 00610 $ V( LDU, * ), VP( * ), WA1( * ), WA2( * ), 00611 $ WA3( * ), WORK( * ), WR( * ), Z( LDU, * ) 00612 * .. 00613 * 00614 * ===================================================================== 00615 * 00616 * .. Parameters .. 00617 REAL ZERO, ONE, TWO, EIGHT, TEN, HUN 00618 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 00619 $ EIGHT = 8.0E0, TEN = 10.0E0, HUN = 100.0E0 ) 00620 REAL HALF 00621 PARAMETER ( HALF = ONE / TWO ) 00622 INTEGER MAXTYP 00623 PARAMETER ( MAXTYP = 21 ) 00624 LOGICAL SRANGE 00625 PARAMETER ( SRANGE = .FALSE. ) 00626 LOGICAL SREL 00627 PARAMETER ( SREL = .FALSE. ) 00628 * .. 00629 * .. Local Scalars .. 00630 LOGICAL BADNN, TRYRAC 00631 INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC, 00632 $ JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC, 00633 $ M, M2, M3, MTYPES, N, NAP, NBLOCK, NERRS, 00634 $ NMATS, NMAX, NSPLIT, NTEST, NTESTT 00635 REAL ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL, 00636 $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, 00637 $ ULPINV, UNFL, VL, VU 00638 * .. 00639 * .. Local Arrays .. 00640 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), 00641 $ KMAGN( MAXTYP ), KMODE( MAXTYP ), 00642 $ KTYPE( MAXTYP ) 00643 REAL DUMMA( 1 ) 00644 * .. 00645 * .. External Functions .. 00646 INTEGER ILAENV 00647 REAL SLAMCH, SLARND, SSXT1 00648 EXTERNAL ILAENV, SLAMCH, SLARND, SSXT1 00649 * .. 00650 * .. External Subroutines .. 00651 EXTERNAL SCOPY, SLABAD, SLACPY, SLASET, SLASUM, SLATMR, 00652 $ SLATMS, SOPGTR, SORGTR, SPTEQR, SSPT21, SSPTRD, 00653 $ SSTEBZ, SSTECH, SSTEDC, SSTEMR, SSTEIN, SSTEQR, 00654 $ SSTERF, SSTT21, SSTT22, SSYT21, SSYTRD, XERBLA 00655 * .. 00656 * .. Intrinsic Functions .. 00657 INTRINSIC ABS, INT, LOG, MAX, MIN, REAL, SQRT 00658 * .. 00659 * .. Data statements .. 00660 DATA KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8, 00661 $ 8, 8, 9, 9, 9, 9, 9, 10 / 00662 DATA KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1, 00663 $ 2, 3, 1, 1, 1, 2, 3, 1 / 00664 DATA KMODE / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0, 00665 $ 0, 0, 4, 3, 1, 4, 4, 3 / 00666 * .. 00667 * .. Executable Statements .. 00668 * 00669 * Keep ftnchek happy 00670 IDUMMA( 1 ) = 1 00671 * 00672 * Check for errors 00673 * 00674 NTESTT = 0 00675 INFO = 0 00676 * 00677 * Important constants 00678 * 00679 BADNN = .FALSE. 00680 TRYRAC = .TRUE. 00681 NMAX = 1 00682 DO 10 J = 1, NSIZES 00683 NMAX = MAX( NMAX, NN( J ) ) 00684 IF( NN( J ).LT.0 ) 00685 $ BADNN = .TRUE. 00686 10 CONTINUE 00687 * 00688 NBLOCK = ILAENV( 1, 'SSYTRD', 'L', NMAX, -1, -1, -1 ) 00689 NBLOCK = MIN( NMAX, MAX( 1, NBLOCK ) ) 00690 * 00691 * Check for errors 00692 * 00693 IF( NSIZES.LT.0 ) THEN 00694 INFO = -1 00695 ELSE IF( BADNN ) THEN 00696 INFO = -2 00697 ELSE IF( NTYPES.LT.0 ) THEN 00698 INFO = -3 00699 ELSE IF( LDA.LT.NMAX ) THEN 00700 INFO = -9 00701 ELSE IF( LDU.LT.NMAX ) THEN 00702 INFO = -23 00703 ELSE IF( 2*MAX( 2, NMAX )**2.GT.LWORK ) THEN 00704 INFO = -29 00705 END IF 00706 * 00707 IF( INFO.NE.0 ) THEN 00708 CALL XERBLA( 'SCHKST', -INFO ) 00709 RETURN 00710 END IF 00711 * 00712 * Quick return if possible 00713 * 00714 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00715 $ RETURN 00716 * 00717 * More Important constants 00718 * 00719 UNFL = SLAMCH( 'Safe minimum' ) 00720 OVFL = ONE / UNFL 00721 CALL SLABAD( UNFL, OVFL ) 00722 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00723 ULPINV = ONE / ULP 00724 LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) ) 00725 RTUNFL = SQRT( UNFL ) 00726 RTOVFL = SQRT( OVFL ) 00727 * 00728 * Loop over sizes, types 00729 * 00730 DO 20 I = 1, 4 00731 ISEED2( I ) = ISEED( I ) 00732 20 CONTINUE 00733 NERRS = 0 00734 NMATS = 0 00735 * 00736 DO 310 JSIZE = 1, NSIZES 00737 N = NN( JSIZE ) 00738 IF( N.GT.0 ) THEN 00739 LGN = INT( LOG( REAL( N ) ) / LOG( TWO ) ) 00740 IF( 2**LGN.LT.N ) 00741 $ LGN = LGN + 1 00742 IF( 2**LGN.LT.N ) 00743 $ LGN = LGN + 1 00744 LWEDC = 1 + 4*N + 2*N*LGN + 4*N**2 00745 LIWEDC = 6 + 6*N + 5*N*LGN 00746 ELSE 00747 LWEDC = 8 00748 LIWEDC = 12 00749 END IF 00750 NAP = ( N*( N+1 ) ) / 2 00751 ANINV = ONE / REAL( MAX( 1, N ) ) 00752 * 00753 IF( NSIZES.NE.1 ) THEN 00754 MTYPES = MIN( MAXTYP, NTYPES ) 00755 ELSE 00756 MTYPES = MIN( MAXTYP+1, NTYPES ) 00757 END IF 00758 * 00759 DO 300 JTYPE = 1, MTYPES 00760 IF( .NOT.DOTYPE( JTYPE ) ) 00761 $ GO TO 300 00762 NMATS = NMATS + 1 00763 NTEST = 0 00764 * 00765 DO 30 J = 1, 4 00766 IOLDSD( J ) = ISEED( J ) 00767 30 CONTINUE 00768 * 00769 * Compute "A" 00770 * 00771 * Control parameters: 00772 * 00773 * KMAGN KMODE KTYPE 00774 * =1 O(1) clustered 1 zero 00775 * =2 large clustered 2 identity 00776 * =3 small exponential (none) 00777 * =4 arithmetic diagonal, (w/ eigenvalues) 00778 * =5 random log symmetric, w/ eigenvalues 00779 * =6 random (none) 00780 * =7 random diagonal 00781 * =8 random symmetric 00782 * =9 positive definite 00783 * =10 diagonally dominant tridiagonal 00784 * 00785 IF( MTYPES.GT.MAXTYP ) 00786 $ GO TO 100 00787 * 00788 ITYPE = KTYPE( JTYPE ) 00789 IMODE = KMODE( JTYPE ) 00790 * 00791 * Compute norm 00792 * 00793 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 00794 * 00795 40 CONTINUE 00796 ANORM = ONE 00797 GO TO 70 00798 * 00799 50 CONTINUE 00800 ANORM = ( RTOVFL*ULP )*ANINV 00801 GO TO 70 00802 * 00803 60 CONTINUE 00804 ANORM = RTUNFL*N*ULPINV 00805 GO TO 70 00806 * 00807 70 CONTINUE 00808 * 00809 CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 00810 IINFO = 0 00811 IF( JTYPE.LE.15 ) THEN 00812 COND = ULPINV 00813 ELSE 00814 COND = ULPINV*ANINV / TEN 00815 END IF 00816 * 00817 * Special Matrices -- Identity & Jordan block 00818 * 00819 * Zero 00820 * 00821 IF( ITYPE.EQ.1 ) THEN 00822 IINFO = 0 00823 * 00824 ELSE IF( ITYPE.EQ.2 ) THEN 00825 * 00826 * Identity 00827 * 00828 DO 80 JC = 1, N 00829 A( JC, JC ) = ANORM 00830 80 CONTINUE 00831 * 00832 ELSE IF( ITYPE.EQ.4 ) THEN 00833 * 00834 * Diagonal Matrix, [Eigen]values Specified 00835 * 00836 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 00837 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ), 00838 $ IINFO ) 00839 * 00840 * 00841 ELSE IF( ITYPE.EQ.5 ) THEN 00842 * 00843 * Symmetric, eigenvalues specified 00844 * 00845 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 00846 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 00847 $ IINFO ) 00848 * 00849 ELSE IF( ITYPE.EQ.7 ) THEN 00850 * 00851 * Diagonal, random eigenvalues 00852 * 00853 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 00854 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00855 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 00856 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00857 * 00858 ELSE IF( ITYPE.EQ.8 ) THEN 00859 * 00860 * Symmetric, random eigenvalues 00861 * 00862 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 00863 $ 'T', 'N', WORK( N+1 ), 1, ONE, 00864 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 00865 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00866 * 00867 ELSE IF( ITYPE.EQ.9 ) THEN 00868 * 00869 * Positive definite, eigenvalues specified. 00870 * 00871 CALL SLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND, 00872 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ), 00873 $ IINFO ) 00874 * 00875 ELSE IF( ITYPE.EQ.10 ) THEN 00876 * 00877 * Positive definite tridiagonal, eigenvalues specified. 00878 * 00879 CALL SLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND, 00880 $ ANORM, 1, 1, 'N', A, LDA, WORK( N+1 ), 00881 $ IINFO ) 00882 DO 90 I = 2, N 00883 TEMP1 = ABS( A( I-1, I ) ) / 00884 $ SQRT( ABS( A( I-1, I-1 )*A( I, I ) ) ) 00885 IF( TEMP1.GT.HALF ) THEN 00886 A( I-1, I ) = HALF*SQRT( ABS( A( I-1, I-1 )*A( I, 00887 $ I ) ) ) 00888 A( I, I-1 ) = A( I-1, I ) 00889 END IF 00890 90 CONTINUE 00891 * 00892 ELSE 00893 * 00894 IINFO = 1 00895 END IF 00896 * 00897 IF( IINFO.NE.0 ) THEN 00898 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 00899 $ IOLDSD 00900 INFO = ABS( IINFO ) 00901 RETURN 00902 END IF 00903 * 00904 100 CONTINUE 00905 * 00906 * Call SSYTRD and SORGTR to compute S and U from 00907 * upper triangle. 00908 * 00909 CALL SLACPY( 'U', N, N, A, LDA, V, LDU ) 00910 * 00911 NTEST = 1 00912 CALL SSYTRD( 'U', N, V, LDU, SD, SE, TAU, WORK, LWORK, 00913 $ IINFO ) 00914 * 00915 IF( IINFO.NE.0 ) THEN 00916 WRITE( NOUNIT, FMT = 9999 )'SSYTRD(U)', IINFO, N, JTYPE, 00917 $ IOLDSD 00918 INFO = ABS( IINFO ) 00919 IF( IINFO.LT.0 ) THEN 00920 RETURN 00921 ELSE 00922 RESULT( 1 ) = ULPINV 00923 GO TO 280 00924 END IF 00925 END IF 00926 * 00927 CALL SLACPY( 'U', N, N, V, LDU, U, LDU ) 00928 * 00929 NTEST = 2 00930 CALL SORGTR( 'U', N, U, LDU, TAU, WORK, LWORK, IINFO ) 00931 IF( IINFO.NE.0 ) THEN 00932 WRITE( NOUNIT, FMT = 9999 )'SORGTR(U)', IINFO, N, JTYPE, 00933 $ IOLDSD 00934 INFO = ABS( IINFO ) 00935 IF( IINFO.LT.0 ) THEN 00936 RETURN 00937 ELSE 00938 RESULT( 2 ) = ULPINV 00939 GO TO 280 00940 END IF 00941 END IF 00942 * 00943 * Do tests 1 and 2 00944 * 00945 CALL SSYT21( 2, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V, 00946 $ LDU, TAU, WORK, RESULT( 1 ) ) 00947 CALL SSYT21( 3, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V, 00948 $ LDU, TAU, WORK, RESULT( 2 ) ) 00949 * 00950 * Call SSYTRD and SORGTR to compute S and U from 00951 * lower triangle, do tests. 00952 * 00953 CALL SLACPY( 'L', N, N, A, LDA, V, LDU ) 00954 * 00955 NTEST = 3 00956 CALL SSYTRD( 'L', N, V, LDU, SD, SE, TAU, WORK, LWORK, 00957 $ IINFO ) 00958 * 00959 IF( IINFO.NE.0 ) THEN 00960 WRITE( NOUNIT, FMT = 9999 )'SSYTRD(L)', IINFO, N, JTYPE, 00961 $ IOLDSD 00962 INFO = ABS( IINFO ) 00963 IF( IINFO.LT.0 ) THEN 00964 RETURN 00965 ELSE 00966 RESULT( 3 ) = ULPINV 00967 GO TO 280 00968 END IF 00969 END IF 00970 * 00971 CALL SLACPY( 'L', N, N, V, LDU, U, LDU ) 00972 * 00973 NTEST = 4 00974 CALL SORGTR( 'L', N, U, LDU, TAU, WORK, LWORK, IINFO ) 00975 IF( IINFO.NE.0 ) THEN 00976 WRITE( NOUNIT, FMT = 9999 )'SORGTR(L)', IINFO, N, JTYPE, 00977 $ IOLDSD 00978 INFO = ABS( IINFO ) 00979 IF( IINFO.LT.0 ) THEN 00980 RETURN 00981 ELSE 00982 RESULT( 4 ) = ULPINV 00983 GO TO 280 00984 END IF 00985 END IF 00986 * 00987 CALL SSYT21( 2, 'Lower', N, 1, A, LDA, SD, SE, U, LDU, V, 00988 $ LDU, TAU, WORK, RESULT( 3 ) ) 00989 CALL SSYT21( 3, 'Lower', N, 1, A, LDA, SD, SE, U, LDU, V, 00990 $ LDU, TAU, WORK, RESULT( 4 ) ) 00991 * 00992 * Store the upper triangle of A in AP 00993 * 00994 I = 0 00995 DO 120 JC = 1, N 00996 DO 110 JR = 1, JC 00997 I = I + 1 00998 AP( I ) = A( JR, JC ) 00999 110 CONTINUE 01000 120 CONTINUE 01001 * 01002 * Call SSPTRD and SOPGTR to compute S and U from AP 01003 * 01004 CALL SCOPY( NAP, AP, 1, VP, 1 ) 01005 * 01006 NTEST = 5 01007 CALL SSPTRD( 'U', N, VP, SD, SE, TAU, IINFO ) 01008 * 01009 IF( IINFO.NE.0 ) THEN 01010 WRITE( NOUNIT, FMT = 9999 )'SSPTRD(U)', IINFO, N, JTYPE, 01011 $ IOLDSD 01012 INFO = ABS( IINFO ) 01013 IF( IINFO.LT.0 ) THEN 01014 RETURN 01015 ELSE 01016 RESULT( 5 ) = ULPINV 01017 GO TO 280 01018 END IF 01019 END IF 01020 * 01021 NTEST = 6 01022 CALL SOPGTR( 'U', N, VP, TAU, U, LDU, WORK, IINFO ) 01023 IF( IINFO.NE.0 ) THEN 01024 WRITE( NOUNIT, FMT = 9999 )'SOPGTR(U)', IINFO, N, JTYPE, 01025 $ IOLDSD 01026 INFO = ABS( IINFO ) 01027 IF( IINFO.LT.0 ) THEN 01028 RETURN 01029 ELSE 01030 RESULT( 6 ) = ULPINV 01031 GO TO 280 01032 END IF 01033 END IF 01034 * 01035 * Do tests 5 and 6 01036 * 01037 CALL SSPT21( 2, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU, 01038 $ WORK, RESULT( 5 ) ) 01039 CALL SSPT21( 3, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU, 01040 $ WORK, RESULT( 6 ) ) 01041 * 01042 * Store the lower triangle of A in AP 01043 * 01044 I = 0 01045 DO 140 JC = 1, N 01046 DO 130 JR = JC, N 01047 I = I + 1 01048 AP( I ) = A( JR, JC ) 01049 130 CONTINUE 01050 140 CONTINUE 01051 * 01052 * Call SSPTRD and SOPGTR to compute S and U from AP 01053 * 01054 CALL SCOPY( NAP, AP, 1, VP, 1 ) 01055 * 01056 NTEST = 7 01057 CALL SSPTRD( 'L', N, VP, SD, SE, TAU, IINFO ) 01058 * 01059 IF( IINFO.NE.0 ) THEN 01060 WRITE( NOUNIT, FMT = 9999 )'SSPTRD(L)', IINFO, N, JTYPE, 01061 $ IOLDSD 01062 INFO = ABS( IINFO ) 01063 IF( IINFO.LT.0 ) THEN 01064 RETURN 01065 ELSE 01066 RESULT( 7 ) = ULPINV 01067 GO TO 280 01068 END IF 01069 END IF 01070 * 01071 NTEST = 8 01072 CALL SOPGTR( 'L', N, VP, TAU, U, LDU, WORK, IINFO ) 01073 IF( IINFO.NE.0 ) THEN 01074 WRITE( NOUNIT, FMT = 9999 )'SOPGTR(L)', IINFO, N, JTYPE, 01075 $ IOLDSD 01076 INFO = ABS( IINFO ) 01077 IF( IINFO.LT.0 ) THEN 01078 RETURN 01079 ELSE 01080 RESULT( 8 ) = ULPINV 01081 GO TO 280 01082 END IF 01083 END IF 01084 * 01085 CALL SSPT21( 2, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU, 01086 $ WORK, RESULT( 7 ) ) 01087 CALL SSPT21( 3, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU, 01088 $ WORK, RESULT( 8 ) ) 01089 * 01090 * Call SSTEQR to compute D1, D2, and Z, do tests. 01091 * 01092 * Compute D1 and Z 01093 * 01094 CALL SCOPY( N, SD, 1, D1, 1 ) 01095 IF( N.GT.0 ) 01096 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01097 CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) 01098 * 01099 NTEST = 9 01100 CALL SSTEQR( 'V', N, D1, WORK, Z, LDU, WORK( N+1 ), IINFO ) 01101 IF( IINFO.NE.0 ) THEN 01102 WRITE( NOUNIT, FMT = 9999 )'SSTEQR(V)', IINFO, N, JTYPE, 01103 $ IOLDSD 01104 INFO = ABS( IINFO ) 01105 IF( IINFO.LT.0 ) THEN 01106 RETURN 01107 ELSE 01108 RESULT( 9 ) = ULPINV 01109 GO TO 280 01110 END IF 01111 END IF 01112 * 01113 * Compute D2 01114 * 01115 CALL SCOPY( N, SD, 1, D2, 1 ) 01116 IF( N.GT.0 ) 01117 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01118 * 01119 NTEST = 11 01120 CALL SSTEQR( 'N', N, D2, WORK, WORK( N+1 ), LDU, 01121 $ WORK( N+1 ), IINFO ) 01122 IF( IINFO.NE.0 ) THEN 01123 WRITE( NOUNIT, FMT = 9999 )'SSTEQR(N)', IINFO, N, JTYPE, 01124 $ IOLDSD 01125 INFO = ABS( IINFO ) 01126 IF( IINFO.LT.0 ) THEN 01127 RETURN 01128 ELSE 01129 RESULT( 11 ) = ULPINV 01130 GO TO 280 01131 END IF 01132 END IF 01133 * 01134 * Compute D3 (using PWK method) 01135 * 01136 CALL SCOPY( N, SD, 1, D3, 1 ) 01137 IF( N.GT.0 ) 01138 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01139 * 01140 NTEST = 12 01141 CALL SSTERF( N, D3, WORK, IINFO ) 01142 IF( IINFO.NE.0 ) THEN 01143 WRITE( NOUNIT, FMT = 9999 )'SSTERF', IINFO, N, JTYPE, 01144 $ IOLDSD 01145 INFO = ABS( IINFO ) 01146 IF( IINFO.LT.0 ) THEN 01147 RETURN 01148 ELSE 01149 RESULT( 12 ) = ULPINV 01150 GO TO 280 01151 END IF 01152 END IF 01153 * 01154 * Do Tests 9 and 10 01155 * 01156 CALL SSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, 01157 $ RESULT( 9 ) ) 01158 * 01159 * Do Tests 11 and 12 01160 * 01161 TEMP1 = ZERO 01162 TEMP2 = ZERO 01163 TEMP3 = ZERO 01164 TEMP4 = ZERO 01165 * 01166 DO 150 J = 1, N 01167 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) ) 01168 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 01169 TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) ) 01170 TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) ) 01171 150 CONTINUE 01172 * 01173 RESULT( 11 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 01174 RESULT( 12 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) ) 01175 * 01176 * Do Test 13 -- Sturm Sequence Test of Eigenvalues 01177 * Go up by factors of two until it succeeds 01178 * 01179 NTEST = 13 01180 TEMP1 = THRESH*( HALF-ULP ) 01181 * 01182 DO 160 J = 0, LOG2UI 01183 CALL SSTECH( N, SD, SE, D1, TEMP1, WORK, IINFO ) 01184 IF( IINFO.EQ.0 ) 01185 $ GO TO 170 01186 TEMP1 = TEMP1*TWO 01187 160 CONTINUE 01188 * 01189 170 CONTINUE 01190 RESULT( 13 ) = TEMP1 01191 * 01192 * For positive definite matrices ( JTYPE.GT.15 ) call SPTEQR 01193 * and do tests 14, 15, and 16 . 01194 * 01195 IF( JTYPE.GT.15 ) THEN 01196 * 01197 * Compute D4 and Z4 01198 * 01199 CALL SCOPY( N, SD, 1, D4, 1 ) 01200 IF( N.GT.0 ) 01201 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01202 CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) 01203 * 01204 NTEST = 14 01205 CALL SPTEQR( 'V', N, D4, WORK, Z, LDU, WORK( N+1 ), 01206 $ IINFO ) 01207 IF( IINFO.NE.0 ) THEN 01208 WRITE( NOUNIT, FMT = 9999 )'SPTEQR(V)', IINFO, N, 01209 $ JTYPE, IOLDSD 01210 INFO = ABS( IINFO ) 01211 IF( IINFO.LT.0 ) THEN 01212 RETURN 01213 ELSE 01214 RESULT( 14 ) = ULPINV 01215 GO TO 280 01216 END IF 01217 END IF 01218 * 01219 * Do Tests 14 and 15 01220 * 01221 CALL SSTT21( N, 0, SD, SE, D4, DUMMA, Z, LDU, WORK, 01222 $ RESULT( 14 ) ) 01223 * 01224 * Compute D5 01225 * 01226 CALL SCOPY( N, SD, 1, D5, 1 ) 01227 IF( N.GT.0 ) 01228 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01229 * 01230 NTEST = 16 01231 CALL SPTEQR( 'N', N, D5, WORK, Z, LDU, WORK( N+1 ), 01232 $ IINFO ) 01233 IF( IINFO.NE.0 ) THEN 01234 WRITE( NOUNIT, FMT = 9999 )'SPTEQR(N)', IINFO, N, 01235 $ JTYPE, IOLDSD 01236 INFO = ABS( IINFO ) 01237 IF( IINFO.LT.0 ) THEN 01238 RETURN 01239 ELSE 01240 RESULT( 16 ) = ULPINV 01241 GO TO 280 01242 END IF 01243 END IF 01244 * 01245 * Do Test 16 01246 * 01247 TEMP1 = ZERO 01248 TEMP2 = ZERO 01249 DO 180 J = 1, N 01250 TEMP1 = MAX( TEMP1, ABS( D4( J ) ), ABS( D5( J ) ) ) 01251 TEMP2 = MAX( TEMP2, ABS( D4( J )-D5( J ) ) ) 01252 180 CONTINUE 01253 * 01254 RESULT( 16 ) = TEMP2 / MAX( UNFL, 01255 $ HUN*ULP*MAX( TEMP1, TEMP2 ) ) 01256 ELSE 01257 RESULT( 14 ) = ZERO 01258 RESULT( 15 ) = ZERO 01259 RESULT( 16 ) = ZERO 01260 END IF 01261 * 01262 * Call SSTEBZ with different options and do tests 17-18. 01263 * 01264 * If S is positive definite and diagonally dominant, 01265 * ask for all eigenvalues with high relative accuracy. 01266 * 01267 VL = ZERO 01268 VU = ZERO 01269 IL = 0 01270 IU = 0 01271 IF( JTYPE.EQ.21 ) THEN 01272 NTEST = 17 01273 ABSTOL = UNFL + UNFL 01274 CALL SSTEBZ( 'A', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, 01275 $ M, NSPLIT, WR, IWORK( 1 ), IWORK( N+1 ), 01276 $ WORK, IWORK( 2*N+1 ), IINFO ) 01277 IF( IINFO.NE.0 ) THEN 01278 WRITE( NOUNIT, FMT = 9999 )'SSTEBZ(A,rel)', IINFO, N, 01279 $ JTYPE, IOLDSD 01280 INFO = ABS( IINFO ) 01281 IF( IINFO.LT.0 ) THEN 01282 RETURN 01283 ELSE 01284 RESULT( 17 ) = ULPINV 01285 GO TO 280 01286 END IF 01287 END IF 01288 * 01289 * Do test 17 01290 * 01291 TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) / 01292 $ ( ONE-HALF )**4 01293 * 01294 TEMP1 = ZERO 01295 DO 190 J = 1, N 01296 TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) / 01297 $ ( ABSTOL+ABS( D4( J ) ) ) ) 01298 190 CONTINUE 01299 * 01300 RESULT( 17 ) = TEMP1 / TEMP2 01301 ELSE 01302 RESULT( 17 ) = ZERO 01303 END IF 01304 * 01305 * Now ask for all eigenvalues with high absolute accuracy. 01306 * 01307 NTEST = 18 01308 ABSTOL = UNFL + UNFL 01309 CALL SSTEBZ( 'A', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, M, 01310 $ NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), WORK, 01311 $ IWORK( 2*N+1 ), IINFO ) 01312 IF( IINFO.NE.0 ) THEN 01313 WRITE( NOUNIT, FMT = 9999 )'SSTEBZ(A)', IINFO, N, JTYPE, 01314 $ IOLDSD 01315 INFO = ABS( IINFO ) 01316 IF( IINFO.LT.0 ) THEN 01317 RETURN 01318 ELSE 01319 RESULT( 18 ) = ULPINV 01320 GO TO 280 01321 END IF 01322 END IF 01323 * 01324 * Do test 18 01325 * 01326 TEMP1 = ZERO 01327 TEMP2 = ZERO 01328 DO 200 J = 1, N 01329 TEMP1 = MAX( TEMP1, ABS( D3( J ) ), ABS( WA1( J ) ) ) 01330 TEMP2 = MAX( TEMP2, ABS( D3( J )-WA1( J ) ) ) 01331 200 CONTINUE 01332 * 01333 RESULT( 18 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 01334 * 01335 * Choose random values for IL and IU, and ask for the 01336 * IL-th through IU-th eigenvalues. 01337 * 01338 NTEST = 19 01339 IF( N.LE.1 ) THEN 01340 IL = 1 01341 IU = N 01342 ELSE 01343 IL = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) 01344 IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) 01345 IF( IU.LT.IL ) THEN 01346 ITEMP = IU 01347 IU = IL 01348 IL = ITEMP 01349 END IF 01350 END IF 01351 * 01352 CALL SSTEBZ( 'I', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, 01353 $ M2, NSPLIT, WA2, IWORK( 1 ), IWORK( N+1 ), 01354 $ WORK, IWORK( 2*N+1 ), IINFO ) 01355 IF( IINFO.NE.0 ) THEN 01356 WRITE( NOUNIT, FMT = 9999 )'SSTEBZ(I)', IINFO, N, JTYPE, 01357 $ IOLDSD 01358 INFO = ABS( IINFO ) 01359 IF( IINFO.LT.0 ) THEN 01360 RETURN 01361 ELSE 01362 RESULT( 19 ) = ULPINV 01363 GO TO 280 01364 END IF 01365 END IF 01366 * 01367 * Determine the values VL and VU of the IL-th and IU-th 01368 * eigenvalues and ask for all eigenvalues in this range. 01369 * 01370 IF( N.GT.0 ) THEN 01371 IF( IL.NE.1 ) THEN 01372 VL = WA1( IL ) - MAX( HALF*( WA1( IL )-WA1( IL-1 ) ), 01373 $ ULP*ANORM, TWO*RTUNFL ) 01374 ELSE 01375 VL = WA1( 1 ) - MAX( HALF*( WA1( N )-WA1( 1 ) ), 01376 $ ULP*ANORM, TWO*RTUNFL ) 01377 END IF 01378 IF( IU.NE.N ) THEN 01379 VU = WA1( IU ) + MAX( HALF*( WA1( IU+1 )-WA1( IU ) ), 01380 $ ULP*ANORM, TWO*RTUNFL ) 01381 ELSE 01382 VU = WA1( N ) + MAX( HALF*( WA1( N )-WA1( 1 ) ), 01383 $ ULP*ANORM, TWO*RTUNFL ) 01384 END IF 01385 ELSE 01386 VL = ZERO 01387 VU = ONE 01388 END IF 01389 * 01390 CALL SSTEBZ( 'V', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, 01391 $ M3, NSPLIT, WA3, IWORK( 1 ), IWORK( N+1 ), 01392 $ WORK, IWORK( 2*N+1 ), IINFO ) 01393 IF( IINFO.NE.0 ) THEN 01394 WRITE( NOUNIT, FMT = 9999 )'SSTEBZ(V)', IINFO, N, JTYPE, 01395 $ IOLDSD 01396 INFO = ABS( IINFO ) 01397 IF( IINFO.LT.0 ) THEN 01398 RETURN 01399 ELSE 01400 RESULT( 19 ) = ULPINV 01401 GO TO 280 01402 END IF 01403 END IF 01404 * 01405 IF( M3.EQ.0 .AND. N.NE.0 ) THEN 01406 RESULT( 19 ) = ULPINV 01407 GO TO 280 01408 END IF 01409 * 01410 * Do test 19 01411 * 01412 TEMP1 = SSXT1( 1, WA2, M2, WA3, M3, ABSTOL, ULP, UNFL ) 01413 TEMP2 = SSXT1( 1, WA3, M3, WA2, M2, ABSTOL, ULP, UNFL ) 01414 IF( N.GT.0 ) THEN 01415 TEMP3 = MAX( ABS( WA1( N ) ), ABS( WA1( 1 ) ) ) 01416 ELSE 01417 TEMP3 = ZERO 01418 END IF 01419 * 01420 RESULT( 19 ) = ( TEMP1+TEMP2 ) / MAX( UNFL, TEMP3*ULP ) 01421 * 01422 * Call SSTEIN to compute eigenvectors corresponding to 01423 * eigenvalues in WA1. (First call SSTEBZ again, to make sure 01424 * it returns these eigenvalues in the correct order.) 01425 * 01426 NTEST = 21 01427 CALL SSTEBZ( 'A', 'B', N, VL, VU, IL, IU, ABSTOL, SD, SE, M, 01428 $ NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), WORK, 01429 $ IWORK( 2*N+1 ), IINFO ) 01430 IF( IINFO.NE.0 ) THEN 01431 WRITE( NOUNIT, FMT = 9999 )'SSTEBZ(A,B)', IINFO, N, 01432 $ JTYPE, IOLDSD 01433 INFO = ABS( IINFO ) 01434 IF( IINFO.LT.0 ) THEN 01435 RETURN 01436 ELSE 01437 RESULT( 20 ) = ULPINV 01438 RESULT( 21 ) = ULPINV 01439 GO TO 280 01440 END IF 01441 END IF 01442 * 01443 CALL SSTEIN( N, SD, SE, M, WA1, IWORK( 1 ), IWORK( N+1 ), Z, 01444 $ LDU, WORK, IWORK( 2*N+1 ), IWORK( 3*N+1 ), 01445 $ IINFO ) 01446 IF( IINFO.NE.0 ) THEN 01447 WRITE( NOUNIT, FMT = 9999 )'SSTEIN', IINFO, N, JTYPE, 01448 $ IOLDSD 01449 INFO = ABS( IINFO ) 01450 IF( IINFO.LT.0 ) THEN 01451 RETURN 01452 ELSE 01453 RESULT( 20 ) = ULPINV 01454 RESULT( 21 ) = ULPINV 01455 GO TO 280 01456 END IF 01457 END IF 01458 * 01459 * Do tests 20 and 21 01460 * 01461 CALL SSTT21( N, 0, SD, SE, WA1, DUMMA, Z, LDU, WORK, 01462 $ RESULT( 20 ) ) 01463 * 01464 * Call SSTEDC(I) to compute D1 and Z, do tests. 01465 * 01466 * Compute D1 and Z 01467 * 01468 CALL SCOPY( N, SD, 1, D1, 1 ) 01469 IF( N.GT.0 ) 01470 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01471 CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) 01472 * 01473 NTEST = 22 01474 CALL SSTEDC( 'I', N, D1, WORK, Z, LDU, WORK( N+1 ), LWEDC-N, 01475 $ IWORK, LIWEDC, IINFO ) 01476 IF( IINFO.NE.0 ) THEN 01477 WRITE( NOUNIT, FMT = 9999 )'SSTEDC(I)', IINFO, N, JTYPE, 01478 $ IOLDSD 01479 INFO = ABS( IINFO ) 01480 IF( IINFO.LT.0 ) THEN 01481 RETURN 01482 ELSE 01483 RESULT( 22 ) = ULPINV 01484 GO TO 280 01485 END IF 01486 END IF 01487 * 01488 * Do Tests 22 and 23 01489 * 01490 CALL SSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, 01491 $ RESULT( 22 ) ) 01492 * 01493 * Call SSTEDC(V) to compute D1 and Z, do tests. 01494 * 01495 * Compute D1 and Z 01496 * 01497 CALL SCOPY( N, SD, 1, D1, 1 ) 01498 IF( N.GT.0 ) 01499 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01500 CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) 01501 * 01502 NTEST = 24 01503 CALL SSTEDC( 'V', N, D1, WORK, Z, LDU, WORK( N+1 ), LWEDC-N, 01504 $ IWORK, LIWEDC, IINFO ) 01505 IF( IINFO.NE.0 ) THEN 01506 WRITE( NOUNIT, FMT = 9999 )'SSTEDC(V)', IINFO, N, JTYPE, 01507 $ IOLDSD 01508 INFO = ABS( IINFO ) 01509 IF( IINFO.LT.0 ) THEN 01510 RETURN 01511 ELSE 01512 RESULT( 24 ) = ULPINV 01513 GO TO 280 01514 END IF 01515 END IF 01516 * 01517 * Do Tests 24 and 25 01518 * 01519 CALL SSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, 01520 $ RESULT( 24 ) ) 01521 * 01522 * Call SSTEDC(N) to compute D2, do tests. 01523 * 01524 * Compute D2 01525 * 01526 CALL SCOPY( N, SD, 1, D2, 1 ) 01527 IF( N.GT.0 ) 01528 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01529 CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) 01530 * 01531 NTEST = 26 01532 CALL SSTEDC( 'N', N, D2, WORK, Z, LDU, WORK( N+1 ), LWEDC-N, 01533 $ IWORK, LIWEDC, IINFO ) 01534 IF( IINFO.NE.0 ) THEN 01535 WRITE( NOUNIT, FMT = 9999 )'SSTEDC(N)', IINFO, N, JTYPE, 01536 $ IOLDSD 01537 INFO = ABS( IINFO ) 01538 IF( IINFO.LT.0 ) THEN 01539 RETURN 01540 ELSE 01541 RESULT( 26 ) = ULPINV 01542 GO TO 280 01543 END IF 01544 END IF 01545 * 01546 * Do Test 26 01547 * 01548 TEMP1 = ZERO 01549 TEMP2 = ZERO 01550 * 01551 DO 210 J = 1, N 01552 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) ) 01553 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 01554 210 CONTINUE 01555 * 01556 RESULT( 26 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 01557 * 01558 * Only test SSTEMR if IEEE compliant 01559 * 01560 IF( ILAENV( 10, 'SSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND. 01561 $ ILAENV( 11, 'SSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN 01562 * 01563 * Call SSTEMR, do test 27 (relative eigenvalue accuracy) 01564 * 01565 * If S is positive definite and diagonally dominant, 01566 * ask for all eigenvalues with high relative accuracy. 01567 * 01568 VL = ZERO 01569 VU = ZERO 01570 IL = 0 01571 IU = 0 01572 IF( JTYPE.EQ.21 .AND. SREL ) THEN 01573 NTEST = 27 01574 ABSTOL = UNFL + UNFL 01575 CALL SSTEMR( 'V', 'A', N, SD, SE, VL, VU, IL, IU, 01576 $ M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC, 01577 $ WORK, LWORK, IWORK( 2*N+1 ), LWORK-2*N, 01578 $ IINFO ) 01579 IF( IINFO.NE.0 ) THEN 01580 WRITE( NOUNIT, FMT = 9999 )'SSTEMR(V,A,rel)', 01581 $ IINFO, N, JTYPE, IOLDSD 01582 INFO = ABS( IINFO ) 01583 IF( IINFO.LT.0 ) THEN 01584 RETURN 01585 ELSE 01586 RESULT( 27 ) = ULPINV 01587 GO TO 270 01588 END IF 01589 END IF 01590 * 01591 * Do test 27 01592 * 01593 TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) / 01594 $ ( ONE-HALF )**4 01595 * 01596 TEMP1 = ZERO 01597 DO 220 J = 1, N 01598 TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) / 01599 $ ( ABSTOL+ABS( D4( J ) ) ) ) 01600 220 CONTINUE 01601 * 01602 RESULT( 27 ) = TEMP1 / TEMP2 01603 * 01604 IL = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) 01605 IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) 01606 IF( IU.LT.IL ) THEN 01607 ITEMP = IU 01608 IU = IL 01609 IL = ITEMP 01610 END IF 01611 * 01612 IF( SRANGE ) THEN 01613 NTEST = 28 01614 ABSTOL = UNFL + UNFL 01615 CALL SSTEMR( 'V', 'I', N, SD, SE, VL, VU, IL, IU, 01616 $ M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC, 01617 $ WORK, LWORK, IWORK( 2*N+1 ), 01618 $ LWORK-2*N, IINFO ) 01619 * 01620 IF( IINFO.NE.0 ) THEN 01621 WRITE( NOUNIT, FMT = 9999 )'SSTEMR(V,I,rel)', 01622 $ IINFO, N, JTYPE, IOLDSD 01623 INFO = ABS( IINFO ) 01624 IF( IINFO.LT.0 ) THEN 01625 RETURN 01626 ELSE 01627 RESULT( 28 ) = ULPINV 01628 GO TO 270 01629 END IF 01630 END IF 01631 * 01632 * 01633 * Do test 28 01634 * 01635 TEMP2 = TWO*( TWO*N-ONE )*ULP* 01636 $ ( ONE+EIGHT*HALF**2 ) / ( ONE-HALF )**4 01637 * 01638 TEMP1 = ZERO 01639 DO 230 J = IL, IU 01640 TEMP1 = MAX( TEMP1, ABS( WR( J-IL+1 )-D4( N-J+ 01641 $ 1 ) ) / ( ABSTOL+ABS( WR( J-IL+1 ) ) ) ) 01642 230 CONTINUE 01643 * 01644 RESULT( 28 ) = TEMP1 / TEMP2 01645 ELSE 01646 RESULT( 28 ) = ZERO 01647 END IF 01648 ELSE 01649 RESULT( 27 ) = ZERO 01650 RESULT( 28 ) = ZERO 01651 END IF 01652 * 01653 * Call SSTEMR(V,I) to compute D1 and Z, do tests. 01654 * 01655 * Compute D1 and Z 01656 * 01657 CALL SCOPY( N, SD, 1, D5, 1 ) 01658 IF( N.GT.0 ) 01659 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01660 CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) 01661 * 01662 IF( SRANGE ) THEN 01663 NTEST = 29 01664 IL = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) 01665 IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) 01666 IF( IU.LT.IL ) THEN 01667 ITEMP = IU 01668 IU = IL 01669 IL = ITEMP 01670 END IF 01671 CALL SSTEMR( 'V', 'I', N, D5, WORK, VL, VU, IL, IU, 01672 $ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC, 01673 $ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ), 01674 $ LIWORK-2*N, IINFO ) 01675 IF( IINFO.NE.0 ) THEN 01676 WRITE( NOUNIT, FMT = 9999 )'SSTEMR(V,I)', IINFO, 01677 $ N, JTYPE, IOLDSD 01678 INFO = ABS( IINFO ) 01679 IF( IINFO.LT.0 ) THEN 01680 RETURN 01681 ELSE 01682 RESULT( 29 ) = ULPINV 01683 GO TO 280 01684 END IF 01685 END IF 01686 * 01687 * Do Tests 29 and 30 01688 * 01689 CALL SSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, 01690 $ M, RESULT( 29 ) ) 01691 * 01692 * Call SSTEMR to compute D2, do tests. 01693 * 01694 * Compute D2 01695 * 01696 CALL SCOPY( N, SD, 1, D5, 1 ) 01697 IF( N.GT.0 ) 01698 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01699 * 01700 NTEST = 31 01701 CALL SSTEMR( 'N', 'I', N, D5, WORK, VL, VU, IL, IU, 01702 $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC, 01703 $ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ), 01704 $ LIWORK-2*N, IINFO ) 01705 IF( IINFO.NE.0 ) THEN 01706 WRITE( NOUNIT, FMT = 9999 )'SSTEMR(N,I)', IINFO, 01707 $ N, JTYPE, IOLDSD 01708 INFO = ABS( IINFO ) 01709 IF( IINFO.LT.0 ) THEN 01710 RETURN 01711 ELSE 01712 RESULT( 31 ) = ULPINV 01713 GO TO 280 01714 END IF 01715 END IF 01716 * 01717 * Do Test 31 01718 * 01719 TEMP1 = ZERO 01720 TEMP2 = ZERO 01721 * 01722 DO 240 J = 1, IU - IL + 1 01723 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), 01724 $ ABS( D2( J ) ) ) 01725 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 01726 240 CONTINUE 01727 * 01728 RESULT( 31 ) = TEMP2 / MAX( UNFL, 01729 $ ULP*MAX( TEMP1, TEMP2 ) ) 01730 * 01731 * 01732 * Call SSTEMR(V,V) to compute D1 and Z, do tests. 01733 * 01734 * Compute D1 and Z 01735 * 01736 CALL SCOPY( N, SD, 1, D5, 1 ) 01737 IF( N.GT.0 ) 01738 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01739 CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) 01740 * 01741 NTEST = 32 01742 * 01743 IF( N.GT.0 ) THEN 01744 IF( IL.NE.1 ) THEN 01745 VL = D2( IL ) - MAX( HALF* 01746 $ ( D2( IL )-D2( IL-1 ) ), ULP*ANORM, 01747 $ TWO*RTUNFL ) 01748 ELSE 01749 VL = D2( 1 ) - MAX( HALF*( D2( N )-D2( 1 ) ), 01750 $ ULP*ANORM, TWO*RTUNFL ) 01751 END IF 01752 IF( IU.NE.N ) THEN 01753 VU = D2( IU ) + MAX( HALF* 01754 $ ( D2( IU+1 )-D2( IU ) ), ULP*ANORM, 01755 $ TWO*RTUNFL ) 01756 ELSE 01757 VU = D2( N ) + MAX( HALF*( D2( N )-D2( 1 ) ), 01758 $ ULP*ANORM, TWO*RTUNFL ) 01759 END IF 01760 ELSE 01761 VL = ZERO 01762 VU = ONE 01763 END IF 01764 * 01765 CALL SSTEMR( 'V', 'V', N, D5, WORK, VL, VU, IL, IU, 01766 $ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC, 01767 $ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ), 01768 $ LIWORK-2*N, IINFO ) 01769 IF( IINFO.NE.0 ) THEN 01770 WRITE( NOUNIT, FMT = 9999 )'SSTEMR(V,V)', IINFO, 01771 $ N, JTYPE, IOLDSD 01772 INFO = ABS( IINFO ) 01773 IF( IINFO.LT.0 ) THEN 01774 RETURN 01775 ELSE 01776 RESULT( 32 ) = ULPINV 01777 GO TO 280 01778 END IF 01779 END IF 01780 * 01781 * Do Tests 32 and 33 01782 * 01783 CALL SSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, 01784 $ M, RESULT( 32 ) ) 01785 * 01786 * Call SSTEMR to compute D2, do tests. 01787 * 01788 * Compute D2 01789 * 01790 CALL SCOPY( N, SD, 1, D5, 1 ) 01791 IF( N.GT.0 ) 01792 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01793 * 01794 NTEST = 34 01795 CALL SSTEMR( 'N', 'V', N, D5, WORK, VL, VU, IL, IU, 01796 $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC, 01797 $ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ), 01798 $ LIWORK-2*N, IINFO ) 01799 IF( IINFO.NE.0 ) THEN 01800 WRITE( NOUNIT, FMT = 9999 )'SSTEMR(N,V)', IINFO, 01801 $ N, JTYPE, IOLDSD 01802 INFO = ABS( IINFO ) 01803 IF( IINFO.LT.0 ) THEN 01804 RETURN 01805 ELSE 01806 RESULT( 34 ) = ULPINV 01807 GO TO 280 01808 END IF 01809 END IF 01810 * 01811 * Do Test 34 01812 * 01813 TEMP1 = ZERO 01814 TEMP2 = ZERO 01815 * 01816 DO 250 J = 1, IU - IL + 1 01817 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), 01818 $ ABS( D2( J ) ) ) 01819 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 01820 250 CONTINUE 01821 * 01822 RESULT( 34 ) = TEMP2 / MAX( UNFL, 01823 $ ULP*MAX( TEMP1, TEMP2 ) ) 01824 ELSE 01825 RESULT( 29 ) = ZERO 01826 RESULT( 30 ) = ZERO 01827 RESULT( 31 ) = ZERO 01828 RESULT( 32 ) = ZERO 01829 RESULT( 33 ) = ZERO 01830 RESULT( 34 ) = ZERO 01831 END IF 01832 * 01833 * 01834 * Call SSTEMR(V,A) to compute D1 and Z, do tests. 01835 * 01836 * Compute D1 and Z 01837 * 01838 CALL SCOPY( N, SD, 1, D5, 1 ) 01839 IF( N.GT.0 ) 01840 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01841 * 01842 NTEST = 35 01843 * 01844 CALL SSTEMR( 'V', 'A', N, D5, WORK, VL, VU, IL, IU, 01845 $ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC, 01846 $ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ), 01847 $ LIWORK-2*N, IINFO ) 01848 IF( IINFO.NE.0 ) THEN 01849 WRITE( NOUNIT, FMT = 9999 )'SSTEMR(V,A)', IINFO, N, 01850 $ JTYPE, IOLDSD 01851 INFO = ABS( IINFO ) 01852 IF( IINFO.LT.0 ) THEN 01853 RETURN 01854 ELSE 01855 RESULT( 35 ) = ULPINV 01856 GO TO 280 01857 END IF 01858 END IF 01859 * 01860 * Do Tests 35 and 36 01861 * 01862 CALL SSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, M, 01863 $ RESULT( 35 ) ) 01864 * 01865 * Call SSTEMR to compute D2, do tests. 01866 * 01867 * Compute D2 01868 * 01869 CALL SCOPY( N, SD, 1, D5, 1 ) 01870 IF( N.GT.0 ) 01871 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 01872 * 01873 NTEST = 37 01874 CALL SSTEMR( 'N', 'A', N, D5, WORK, VL, VU, IL, IU, 01875 $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC, 01876 $ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ), 01877 $ LIWORK-2*N, IINFO ) 01878 IF( IINFO.NE.0 ) THEN 01879 WRITE( NOUNIT, FMT = 9999 )'SSTEMR(N,A)', IINFO, N, 01880 $ JTYPE, IOLDSD 01881 INFO = ABS( IINFO ) 01882 IF( IINFO.LT.0 ) THEN 01883 RETURN 01884 ELSE 01885 RESULT( 37 ) = ULPINV 01886 GO TO 280 01887 END IF 01888 END IF 01889 * 01890 * Do Test 34 01891 * 01892 TEMP1 = ZERO 01893 TEMP2 = ZERO 01894 * 01895 DO 260 J = 1, N 01896 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) ) 01897 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 01898 260 CONTINUE 01899 * 01900 RESULT( 37 ) = TEMP2 / MAX( UNFL, 01901 $ ULP*MAX( TEMP1, TEMP2 ) ) 01902 END IF 01903 270 CONTINUE 01904 280 CONTINUE 01905 NTESTT = NTESTT + NTEST 01906 * 01907 * End of Loop -- Check for RESULT(j) > THRESH 01908 * 01909 * 01910 * Print out tests which fail. 01911 * 01912 DO 290 JR = 1, NTEST 01913 IF( RESULT( JR ).GE.THRESH ) THEN 01914 * 01915 * If this is the first test to fail, 01916 * print a header to the data file. 01917 * 01918 IF( NERRS.EQ.0 ) THEN 01919 WRITE( NOUNIT, FMT = 9998 )'SST' 01920 WRITE( NOUNIT, FMT = 9997 ) 01921 WRITE( NOUNIT, FMT = 9996 ) 01922 WRITE( NOUNIT, FMT = 9995 )'Symmetric' 01923 WRITE( NOUNIT, FMT = 9994 ) 01924 * 01925 * Tests performed 01926 * 01927 WRITE( NOUNIT, FMT = 9988 ) 01928 END IF 01929 NERRS = NERRS + 1 01930 WRITE( NOUNIT, FMT = 9990 )N, IOLDSD, JTYPE, JR, 01931 $ RESULT( JR ) 01932 END IF 01933 290 CONTINUE 01934 300 CONTINUE 01935 310 CONTINUE 01936 * 01937 * Summary 01938 * 01939 CALL SLASUM( 'SST', NOUNIT, NERRS, NTESTT ) 01940 RETURN 01941 * 01942 9999 FORMAT( ' SCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 01943 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 01944 * 01945 9998 FORMAT( / 1X, A3, ' -- Real Symmetric eigenvalue problem' ) 01946 9997 FORMAT( ' Matrix types (see SCHKST for details): ' ) 01947 * 01948 9996 FORMAT( / ' Special Matrices:', 01949 $ / ' 1=Zero matrix. ', 01950 $ ' 5=Diagonal: clustered entries.', 01951 $ / ' 2=Identity matrix. ', 01952 $ ' 6=Diagonal: large, evenly spaced.', 01953 $ / ' 3=Diagonal: evenly spaced entries. ', 01954 $ ' 7=Diagonal: small, evenly spaced.', 01955 $ / ' 4=Diagonal: geometr. spaced entries.' ) 01956 9995 FORMAT( ' Dense ', A, ' Matrices:', 01957 $ / ' 8=Evenly spaced eigenvals. ', 01958 $ ' 12=Small, evenly spaced eigenvals.', 01959 $ / ' 9=Geometrically spaced eigenvals. ', 01960 $ ' 13=Matrix with random O(1) entries.', 01961 $ / ' 10=Clustered eigenvalues. ', 01962 $ ' 14=Matrix with large random entries.', 01963 $ / ' 11=Large, evenly spaced eigenvals. ', 01964 $ ' 15=Matrix with small random entries.' ) 01965 9994 FORMAT( ' 16=Positive definite, evenly spaced eigenvalues', 01966 $ / ' 17=Positive definite, geometrically spaced eigenvlaues', 01967 $ / ' 18=Positive definite, clustered eigenvalues', 01968 $ / ' 19=Positive definite, small evenly spaced eigenvalues', 01969 $ / ' 20=Positive definite, large evenly spaced eigenvalues', 01970 $ / ' 21=Diagonally dominant tridiagonal, geometrically', 01971 $ ' spaced eigenvalues' ) 01972 * 01973 9990 FORMAT( ' N=', I5, ', seed=', 4( I4, ',' ), ' type ', I2, 01974 $ ', test(', I2, ')=', G10.3 ) 01975 * 01976 9988 FORMAT( / 'Test performed: see SCHKST for details.', / ) 01977 * End of SCHKST 01978 * 01979 END