1*> \brief \b ZCHKST2STG 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE ZCHKST2STG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 12* NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5, 13* WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK, 14* LWORK, RWORK, LRWORK, IWORK, LIWORK, RESULT, 15* INFO ) 16* 17* .. Scalar Arguments .. 18* INTEGER INFO, LDA, LDU, LIWORK, LRWORK, LWORK, NOUNIT, 19* $ NSIZES, NTYPES 20* DOUBLE PRECISION THRESH 21* .. 22* .. Array Arguments .. 23* LOGICAL DOTYPE( * ) 24* INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 25* DOUBLE PRECISION D1( * ), D2( * ), D3( * ), D4( * ), D5( * ), 26* $ RESULT( * ), RWORK( * ), SD( * ), SE( * ), 27* $ WA1( * ), WA2( * ), WA3( * ), WR( * ) 28* COMPLEX*16 A( LDA, * ), AP( * ), TAU( * ), U( LDU, * ), 29* $ V( LDU, * ), VP( * ), WORK( * ), Z( LDU, * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> ZCHKST2STG checks the Hermitian eigenvalue problem routines 39*> using the 2-stage reduction techniques. Since the generation 40*> of Q or the vectors is not available in this release, we only 41*> compare the eigenvalue resulting when using the 2-stage to the 42*> one considered as reference using the standard 1-stage reduction 43*> ZHETRD. For that, we call the standard ZHETRD and compute D1 using 44*> DSTEQR, then we call the 2-stage ZHETRD_2STAGE with Upper and Lower 45*> and we compute D2 and D3 using DSTEQR and then we replaced tests 46*> 3 and 4 by tests 11 and 12. test 1 and 2 remain to verify that 47*> the 1-stage results are OK and can be trusted. 48*> This testing routine will converge to the ZCHKST in the next 49*> release when vectors and generation of Q will be implemented. 50*> 51*> ZHETRD factors A as U S U* , where * means conjugate transpose, 52*> S is real symmetric tridiagonal, and U is unitary. 53*> ZHETRD can use either just the lower or just the upper triangle 54*> of A; ZCHKST2STG checks both cases. 55*> U is represented as a product of Householder 56*> transformations, whose vectors are stored in the first 57*> n-1 columns of V, and whose scale factors are in TAU. 58*> 59*> ZHPTRD does the same as ZHETRD, except that A and V are stored 60*> in "packed" format. 61*> 62*> ZUNGTR constructs the matrix U from the contents of V and TAU. 63*> 64*> ZUPGTR constructs the matrix U from the contents of VP and TAU. 65*> 66*> ZSTEQR factors S as Z D1 Z* , where Z is the unitary 67*> matrix of eigenvectors and D1 is a diagonal matrix with 68*> the eigenvalues on the diagonal. D2 is the matrix of 69*> eigenvalues computed when Z is not computed. 70*> 71*> DSTERF computes D3, the matrix of eigenvalues, by the 72*> PWK method, which does not yield eigenvectors. 73*> 74*> ZPTEQR factors S as Z4 D4 Z4* , for a 75*> Hermitian positive definite tridiagonal matrix. 76*> D5 is the matrix of eigenvalues computed when Z is not 77*> computed. 78*> 79*> DSTEBZ computes selected eigenvalues. WA1, WA2, and 80*> WA3 will denote eigenvalues computed to high 81*> absolute accuracy, with different range options. 82*> WR will denote eigenvalues computed to high relative 83*> accuracy. 84*> 85*> ZSTEIN computes Y, the eigenvectors of S, given the 86*> eigenvalues. 87*> 88*> ZSTEDC factors S as Z D1 Z* , where Z is the unitary 89*> matrix of eigenvectors and D1 is a diagonal matrix with 90*> the eigenvalues on the diagonal ('I' option). It may also 91*> update an input unitary matrix, usually the output 92*> from ZHETRD/ZUNGTR or ZHPTRD/ZUPGTR ('V' option). It may 93*> also just compute eigenvalues ('N' option). 94*> 95*> ZSTEMR factors S as Z D1 Z* , where Z is the unitary 96*> matrix of eigenvectors and D1 is a diagonal matrix with 97*> the eigenvalues on the diagonal ('I' option). ZSTEMR 98*> uses the Relatively Robust Representation whenever possible. 99*> 100*> When ZCHKST2STG is called, a number of matrix "sizes" ("n's") and a 101*> number of matrix "types" are specified. For each size ("n") 102*> and each type of matrix, one matrix will be generated and used 103*> to test the Hermitian eigenroutines. For each matrix, a number 104*> of tests will be performed: 105*> 106*> (1) | A - V S V* | / ( |A| n ulp ) ZHETRD( UPLO='U', ... ) 107*> 108*> (2) | I - UV* | / ( n ulp ) ZUNGTR( UPLO='U', ... ) 109*> 110*> (3) | A - V S V* | / ( |A| n ulp ) ZHETRD( UPLO='L', ... ) 111*> replaced by | D1 - D2 | / ( |D1| ulp ) where D1 is the 112*> eigenvalue matrix computed using S and D2 is the 113*> eigenvalue matrix computed using S_2stage the output of 114*> ZHETRD_2STAGE("N", "U",....). D1 and D2 are computed 115*> via DSTEQR('N',...) 116*> 117*> (4) | I - UV* | / ( n ulp ) ZUNGTR( UPLO='L', ... ) 118*> replaced by | D1 - D3 | / ( |D1| ulp ) where D1 is the 119*> eigenvalue matrix computed using S and D3 is the 120*> eigenvalue matrix computed using S_2stage the output of 121*> ZHETRD_2STAGE("N", "L",....). D1 and D3 are computed 122*> via DSTEQR('N',...) 123*> 124*> (5-8) Same as 1-4, but for ZHPTRD and ZUPGTR. 125*> 126*> (9) | S - Z D Z* | / ( |S| n ulp ) ZSTEQR('V',...) 127*> 128*> (10) | I - ZZ* | / ( n ulp ) ZSTEQR('V',...) 129*> 130*> (11) | D1 - D2 | / ( |D1| ulp ) ZSTEQR('N',...) 131*> 132*> (12) | D1 - D3 | / ( |D1| ulp ) DSTERF 133*> 134*> (13) 0 if the true eigenvalues (computed by sturm count) 135*> of S are within THRESH of 136*> those in D1. 2*THRESH if they are not. (Tested using 137*> DSTECH) 138*> 139*> For S positive definite, 140*> 141*> (14) | S - Z4 D4 Z4* | / ( |S| n ulp ) ZPTEQR('V',...) 142*> 143*> (15) | I - Z4 Z4* | / ( n ulp ) ZPTEQR('V',...) 144*> 145*> (16) | D4 - D5 | / ( 100 |D4| ulp ) ZPTEQR('N',...) 146*> 147*> When S is also diagonally dominant by the factor gamma < 1, 148*> 149*> (17) max | D4(i) - WR(i) | / ( |D4(i)| omega ) , 150*> i 151*> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 152*> DSTEBZ( 'A', 'E', ...) 153*> 154*> (18) | WA1 - D3 | / ( |D3| ulp ) DSTEBZ( 'A', 'E', ...) 155*> 156*> (19) ( max { min | WA2(i)-WA3(j) | } + 157*> i j 158*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) 159*> i j 160*> DSTEBZ( 'I', 'E', ...) 161*> 162*> (20) | S - Y WA1 Y* | / ( |S| n ulp ) DSTEBZ, ZSTEIN 163*> 164*> (21) | I - Y Y* | / ( n ulp ) DSTEBZ, ZSTEIN 165*> 166*> (22) | S - Z D Z* | / ( |S| n ulp ) ZSTEDC('I') 167*> 168*> (23) | I - ZZ* | / ( n ulp ) ZSTEDC('I') 169*> 170*> (24) | S - Z D Z* | / ( |S| n ulp ) ZSTEDC('V') 171*> 172*> (25) | I - ZZ* | / ( n ulp ) ZSTEDC('V') 173*> 174*> (26) | D1 - D2 | / ( |D1| ulp ) ZSTEDC('V') and 175*> ZSTEDC('N') 176*> 177*> Test 27 is disabled at the moment because ZSTEMR does not 178*> guarantee high relatvie accuracy. 179*> 180*> (27) max | D6(i) - WR(i) | / ( |D6(i)| omega ) , 181*> i 182*> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 183*> ZSTEMR('V', 'A') 184*> 185*> (28) max | D6(i) - WR(i) | / ( |D6(i)| omega ) , 186*> i 187*> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 188*> ZSTEMR('V', 'I') 189*> 190*> Tests 29 through 34 are disable at present because ZSTEMR 191*> does not handle partial spectrum requests. 192*> 193*> (29) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'I') 194*> 195*> (30) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'I') 196*> 197*> (31) ( max { min | WA2(i)-WA3(j) | } + 198*> i j 199*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) 200*> i j 201*> ZSTEMR('N', 'I') vs. CSTEMR('V', 'I') 202*> 203*> (32) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'V') 204*> 205*> (33) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'V') 206*> 207*> (34) ( max { min | WA2(i)-WA3(j) | } + 208*> i j 209*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) 210*> i j 211*> ZSTEMR('N', 'V') vs. CSTEMR('V', 'V') 212*> 213*> (35) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'A') 214*> 215*> (36) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'A') 216*> 217*> (37) ( max { min | WA2(i)-WA3(j) | } + 218*> i j 219*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) 220*> i j 221*> ZSTEMR('N', 'A') vs. CSTEMR('V', 'A') 222*> 223*> The "sizes" are specified by an array NN(1:NSIZES); the value of 224*> each element NN(j) specifies one size. 225*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); 226*> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 227*> Currently, the list of possible types is: 228*> 229*> (1) The zero matrix. 230*> (2) The identity matrix. 231*> 232*> (3) A diagonal matrix with evenly spaced entries 233*> 1, ..., ULP and random signs. 234*> (ULP = (first number larger than 1) - 1 ) 235*> (4) A diagonal matrix with geometrically spaced entries 236*> 1, ..., ULP and random signs. 237*> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 238*> and random signs. 239*> 240*> (6) Same as (4), but multiplied by SQRT( overflow threshold ) 241*> (7) Same as (4), but multiplied by SQRT( underflow threshold ) 242*> 243*> (8) A matrix of the form U* D U, where U is unitary and 244*> D has evenly spaced entries 1, ..., ULP with random signs 245*> on the diagonal. 246*> 247*> (9) A matrix of the form U* D U, where U is unitary and 248*> D has geometrically spaced entries 1, ..., ULP with random 249*> signs on the diagonal. 250*> 251*> (10) A matrix of the form U* D U, where U is unitary and 252*> D has "clustered" entries 1, ULP,..., ULP with random 253*> signs on the diagonal. 254*> 255*> (11) Same as (8), but multiplied by SQRT( overflow threshold ) 256*> (12) Same as (8), but multiplied by SQRT( underflow threshold ) 257*> 258*> (13) Hermitian matrix with random entries chosen from (-1,1). 259*> (14) Same as (13), but multiplied by SQRT( overflow threshold ) 260*> (15) Same as (13), but multiplied by SQRT( underflow threshold ) 261*> (16) Same as (8), but diagonal elements are all positive. 262*> (17) Same as (9), but diagonal elements are all positive. 263*> (18) Same as (10), but diagonal elements are all positive. 264*> (19) Same as (16), but multiplied by SQRT( overflow threshold ) 265*> (20) Same as (16), but multiplied by SQRT( underflow threshold ) 266*> (21) A diagonally dominant tridiagonal matrix with geometrically 267*> spaced diagonal entries 1, ..., ULP. 268*> \endverbatim 269* 270* Arguments: 271* ========== 272* 273*> \param[in] NSIZES 274*> \verbatim 275*> NSIZES is INTEGER 276*> The number of sizes of matrices to use. If it is zero, 277*> ZCHKST2STG does nothing. It must be at least zero. 278*> \endverbatim 279*> 280*> \param[in] NN 281*> \verbatim 282*> NN is INTEGER array, dimension (NSIZES) 283*> An array containing the sizes to be used for the matrices. 284*> Zero values will be skipped. The values must be at least 285*> zero. 286*> \endverbatim 287*> 288*> \param[in] NTYPES 289*> \verbatim 290*> NTYPES is INTEGER 291*> The number of elements in DOTYPE. If it is zero, ZCHKST2STG 292*> does nothing. It must be at least zero. If it is MAXTYP+1 293*> and NSIZES is 1, then an additional type, MAXTYP+1 is 294*> defined, which is to use whatever matrix is in A. This 295*> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 296*> DOTYPE(MAXTYP+1) is .TRUE. . 297*> \endverbatim 298*> 299*> \param[in] DOTYPE 300*> \verbatim 301*> DOTYPE is LOGICAL array, dimension (NTYPES) 302*> If DOTYPE(j) is .TRUE., then for each size in NN a 303*> matrix of that size and of type j will be generated. 304*> If NTYPES is smaller than the maximum number of types 305*> defined (PARAMETER MAXTYP), then types NTYPES+1 through 306*> MAXTYP will not be generated. If NTYPES is larger 307*> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 308*> will be ignored. 309*> \endverbatim 310*> 311*> \param[in,out] ISEED 312*> \verbatim 313*> ISEED is INTEGER array, dimension (4) 314*> On entry ISEED specifies the seed of the random number 315*> generator. The array elements should be between 0 and 4095; 316*> if not they will be reduced mod 4096. Also, ISEED(4) must 317*> be odd. The random number generator uses a linear 318*> congruential sequence limited to small integers, and so 319*> should produce machine independent random numbers. The 320*> values of ISEED are changed on exit, and can be used in the 321*> next call to ZCHKST2STG to continue the same random number 322*> sequence. 323*> \endverbatim 324*> 325*> \param[in] THRESH 326*> \verbatim 327*> THRESH is DOUBLE PRECISION 328*> A test will count as "failed" if the "error", computed as 329*> described above, exceeds THRESH. Note that the error 330*> is scaled to be O(1), so THRESH should be a reasonably 331*> small multiple of 1, e.g., 10 or 100. In particular, 332*> it should not depend on the precision (single vs. double) 333*> or the size of the matrix. It must be at least zero. 334*> \endverbatim 335*> 336*> \param[in] NOUNIT 337*> \verbatim 338*> NOUNIT is INTEGER 339*> The FORTRAN unit number for printing out error messages 340*> (e.g., if a routine returns IINFO not equal to 0.) 341*> \endverbatim 342*> 343*> \param[in,out] A 344*> \verbatim 345*> A is COMPLEX*16 array of 346*> dimension ( LDA , max(NN) ) 347*> Used to hold the matrix whose eigenvalues are to be 348*> computed. On exit, A contains the last matrix actually 349*> used. 350*> \endverbatim 351*> 352*> \param[in] LDA 353*> \verbatim 354*> LDA is INTEGER 355*> The leading dimension of A. It must be at 356*> least 1 and at least max( NN ). 357*> \endverbatim 358*> 359*> \param[out] AP 360*> \verbatim 361*> AP is COMPLEX*16 array of 362*> dimension( max(NN)*max(NN+1)/2 ) 363*> The matrix A stored in packed format. 364*> \endverbatim 365*> 366*> \param[out] SD 367*> \verbatim 368*> SD is DOUBLE PRECISION array of 369*> dimension( max(NN) ) 370*> The diagonal of the tridiagonal matrix computed by ZHETRD. 371*> On exit, SD and SE contain the tridiagonal form of the 372*> matrix in A. 373*> \endverbatim 374*> 375*> \param[out] SE 376*> \verbatim 377*> SE is DOUBLE PRECISION array of 378*> dimension( max(NN) ) 379*> The off-diagonal of the tridiagonal matrix computed by 380*> ZHETRD. On exit, SD and SE contain the tridiagonal form of 381*> the matrix in A. 382*> \endverbatim 383*> 384*> \param[out] D1 385*> \verbatim 386*> D1 is DOUBLE PRECISION array of 387*> dimension( max(NN) ) 388*> The eigenvalues of A, as computed by ZSTEQR simlutaneously 389*> with Z. On exit, the eigenvalues in D1 correspond with the 390*> matrix in A. 391*> \endverbatim 392*> 393*> \param[out] D2 394*> \verbatim 395*> D2 is DOUBLE PRECISION array of 396*> dimension( max(NN) ) 397*> The eigenvalues of A, as computed by ZSTEQR if Z is not 398*> computed. On exit, the eigenvalues in D2 correspond with 399*> the matrix in A. 400*> \endverbatim 401*> 402*> \param[out] D3 403*> \verbatim 404*> D3 is DOUBLE PRECISION array of 405*> dimension( max(NN) ) 406*> The eigenvalues of A, as computed by DSTERF. On exit, the 407*> eigenvalues in D3 correspond with the matrix in A. 408*> \endverbatim 409*> 410*> \param[out] D4 411*> \verbatim 412*> D4 is DOUBLE PRECISION array of 413*> dimension( max(NN) ) 414*> The eigenvalues of A, as computed by ZPTEQR(V). 415*> ZPTEQR factors S as Z4 D4 Z4* 416*> On exit, the eigenvalues in D4 correspond with the matrix in A. 417*> \endverbatim 418*> 419*> \param[out] D5 420*> \verbatim 421*> D5 is DOUBLE PRECISION array of 422*> dimension( max(NN) ) 423*> The eigenvalues of A, as computed by ZPTEQR(N) 424*> when Z is not computed. On exit, the 425*> eigenvalues in D4 correspond with the matrix in A. 426*> \endverbatim 427*> 428*> \param[out] WA1 429*> \verbatim 430*> WA1 is DOUBLE PRECISION array of 431*> dimension( max(NN) ) 432*> All eigenvalues of A, computed to high 433*> absolute accuracy, with different range options. 434*> as computed by DSTEBZ. 435*> \endverbatim 436*> 437*> \param[out] WA2 438*> \verbatim 439*> WA2 is DOUBLE PRECISION array of 440*> dimension( max(NN) ) 441*> Selected eigenvalues of A, computed to high 442*> absolute accuracy, with different range options. 443*> as computed by DSTEBZ. 444*> Choose random values for IL and IU, and ask for the 445*> IL-th through IU-th eigenvalues. 446*> \endverbatim 447*> 448*> \param[out] WA3 449*> \verbatim 450*> WA3 is DOUBLE PRECISION array of 451*> dimension( max(NN) ) 452*> Selected eigenvalues of A, computed to high 453*> absolute accuracy, with different range options. 454*> as computed by DSTEBZ. 455*> Determine the values VL and VU of the IL-th and IU-th 456*> eigenvalues and ask for all eigenvalues in this range. 457*> \endverbatim 458*> 459*> \param[out] WR 460*> \verbatim 461*> WR is DOUBLE PRECISION array of 462*> dimension( max(NN) ) 463*> All eigenvalues of A, computed to high 464*> absolute accuracy, with different options. 465*> as computed by DSTEBZ. 466*> \endverbatim 467*> 468*> \param[out] U 469*> \verbatim 470*> U is COMPLEX*16 array of 471*> dimension( LDU, max(NN) ). 472*> The unitary matrix computed by ZHETRD + ZUNGTR. 473*> \endverbatim 474*> 475*> \param[in] LDU 476*> \verbatim 477*> LDU is INTEGER 478*> The leading dimension of U, Z, and V. It must be at least 1 479*> and at least max( NN ). 480*> \endverbatim 481*> 482*> \param[out] V 483*> \verbatim 484*> V is COMPLEX*16 array of 485*> dimension( LDU, max(NN) ). 486*> The Housholder vectors computed by ZHETRD in reducing A to 487*> tridiagonal form. The vectors computed with UPLO='U' are 488*> in the upper triangle, and the vectors computed with UPLO='L' 489*> are in the lower triangle. (As described in ZHETRD, the 490*> sub- and superdiagonal are not set to 1, although the 491*> true Householder vector has a 1 in that position. The 492*> routines that use V, such as ZUNGTR, set those entries to 493*> 1 before using them, and then restore them later.) 494*> \endverbatim 495*> 496*> \param[out] VP 497*> \verbatim 498*> VP is COMPLEX*16 array of 499*> dimension( max(NN)*max(NN+1)/2 ) 500*> The matrix V stored in packed format. 501*> \endverbatim 502*> 503*> \param[out] TAU 504*> \verbatim 505*> TAU is COMPLEX*16 array of 506*> dimension( max(NN) ) 507*> The Householder factors computed by ZHETRD in reducing A 508*> to tridiagonal form. 509*> \endverbatim 510*> 511*> \param[out] Z 512*> \verbatim 513*> Z is COMPLEX*16 array of 514*> dimension( LDU, max(NN) ). 515*> The unitary matrix of eigenvectors computed by ZSTEQR, 516*> ZPTEQR, and ZSTEIN. 517*> \endverbatim 518*> 519*> \param[out] WORK 520*> \verbatim 521*> WORK is COMPLEX*16 array of 522*> dimension( LWORK ) 523*> \endverbatim 524*> 525*> \param[in] LWORK 526*> \verbatim 527*> LWORK is INTEGER 528*> The number of entries in WORK. This must be at least 529*> 1 + 4 * Nmax + 2 * Nmax * lg Nmax + 3 * Nmax**2 530*> where Nmax = max( NN(j), 2 ) and lg = log base 2. 531*> \endverbatim 532*> 533*> \param[out] IWORK 534*> \verbatim 535*> IWORK is INTEGER array, 536*> Workspace. 537*> \endverbatim 538*> 539*> \param[out] LIWORK 540*> \verbatim 541*> LIWORK is INTEGER 542*> The number of entries in IWORK. This must be at least 543*> 6 + 6*Nmax + 5 * Nmax * lg Nmax 544*> where Nmax = max( NN(j), 2 ) and lg = log base 2. 545*> \endverbatim 546*> 547*> \param[out] RWORK 548*> \verbatim 549*> RWORK is DOUBLE PRECISION array 550*> \endverbatim 551*> 552*> \param[in] LRWORK 553*> \verbatim 554*> LRWORK is INTEGER 555*> The number of entries in LRWORK (dimension( ??? ) 556*> \endverbatim 557*> 558*> \param[out] RESULT 559*> \verbatim 560*> RESULT is DOUBLE PRECISION array, dimension (26) 561*> The values computed by the tests described above. 562*> The values are currently limited to 1/ulp, to avoid 563*> overflow. 564*> \endverbatim 565*> 566*> \param[out] INFO 567*> \verbatim 568*> INFO is INTEGER 569*> If 0, then everything ran OK. 570*> -1: NSIZES < 0 571*> -2: Some NN(j) < 0 572*> -3: NTYPES < 0 573*> -5: THRESH < 0 574*> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). 575*> -23: LDU < 1 or LDU < NMAX. 576*> -29: LWORK too small. 577*> If ZLATMR, CLATMS, ZHETRD, ZUNGTR, ZSTEQR, DSTERF, 578*> or ZUNMC2 returns an error code, the 579*> absolute value of it is returned. 580*> 581*>----------------------------------------------------------------------- 582*> 583*> Some Local Variables and Parameters: 584*> ---- ----- --------- --- ---------- 585*> ZERO, ONE Real 0 and 1. 586*> MAXTYP The number of types defined. 587*> NTEST The number of tests performed, or which can 588*> be performed so far, for the current matrix. 589*> NTESTT The total number of tests performed so far. 590*> NBLOCK Blocksize as returned by ENVIR. 591*> NMAX Largest value in NN. 592*> NMATS The number of matrices generated so far. 593*> NERRS The number of tests which have exceeded THRESH 594*> so far. 595*> COND, IMODE Values to be passed to the matrix generators. 596*> ANORM Norm of A; passed to matrix generators. 597*> 598*> OVFL, UNFL Overflow and underflow thresholds. 599*> ULP, ULPINV Finest relative precision and its inverse. 600*> RTOVFL, RTUNFL Square roots of the previous 2 values. 601*> The following four arrays decode JTYPE: 602*> KTYPE(j) The general type (1-10) for type "j". 603*> KMODE(j) The MODE value to be passed to the matrix 604*> generator for type "j". 605*> KMAGN(j) The order of magnitude ( O(1), 606*> O(overflow^(1/2) ), O(underflow^(1/2) ) 607*> \endverbatim 608* 609* Authors: 610* ======== 611* 612*> \author Univ. of Tennessee 613*> \author Univ. of California Berkeley 614*> \author Univ. of Colorado Denver 615*> \author NAG Ltd. 616* 617*> \ingroup complex16_eig 618* 619* ===================================================================== 620 SUBROUTINE ZCHKST2STG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 621 $ NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5, 622 $ WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK, 623 $ LWORK, RWORK, LRWORK, IWORK, LIWORK, RESULT, 624 $ INFO ) 625* 626* -- LAPACK test routine -- 627* -- LAPACK is a software package provided by Univ. of Tennessee, -- 628* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 629* 630* .. Scalar Arguments .. 631 INTEGER INFO, LDA, LDU, LIWORK, LRWORK, LWORK, NOUNIT, 632 $ NSIZES, NTYPES 633 DOUBLE PRECISION THRESH 634* .. 635* .. Array Arguments .. 636 LOGICAL DOTYPE( * ) 637 INTEGER ISEED( 4 ), IWORK( * ), NN( * ) 638 DOUBLE PRECISION D1( * ), D2( * ), D3( * ), D4( * ), D5( * ), 639 $ RESULT( * ), RWORK( * ), SD( * ), SE( * ), 640 $ WA1( * ), WA2( * ), WA3( * ), WR( * ) 641 COMPLEX*16 A( LDA, * ), AP( * ), TAU( * ), U( LDU, * ), 642 $ V( LDU, * ), VP( * ), WORK( * ), Z( LDU, * ) 643* .. 644* 645* ===================================================================== 646* 647* .. Parameters .. 648 DOUBLE PRECISION ZERO, ONE, TWO, EIGHT, TEN, HUN 649 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 650 $ EIGHT = 8.0D0, TEN = 10.0D0, HUN = 100.0D0 ) 651 COMPLEX*16 CZERO, CONE 652 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 653 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 654 DOUBLE PRECISION HALF 655 PARAMETER ( HALF = ONE / TWO ) 656 INTEGER MAXTYP 657 PARAMETER ( MAXTYP = 21 ) 658 LOGICAL CRANGE 659 PARAMETER ( CRANGE = .FALSE. ) 660 LOGICAL CREL 661 PARAMETER ( CREL = .FALSE. ) 662* .. 663* .. Local Scalars .. 664 LOGICAL BADNN, TRYRAC 665 INTEGER I, IINFO, IL, IMODE, INDE, INDRWK, ITEMP, 666 $ ITYPE, IU, J, JC, JR, JSIZE, JTYPE, LGN, 667 $ LIWEDC, LOG2UI, LRWEDC, LWEDC, M, M2, M3, 668 $ MTYPES, N, NAP, NBLOCK, NERRS, NMATS, NMAX, 669 $ NSPLIT, NTEST, NTESTT, LH, LW 670 DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL, 671 $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, 672 $ ULPINV, UNFL, VL, VU 673* .. 674* .. Local Arrays .. 675 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), 676 $ KMAGN( MAXTYP ), KMODE( MAXTYP ), 677 $ KTYPE( MAXTYP ) 678 DOUBLE PRECISION DUMMA( 1 ) 679* .. 680* .. External Functions .. 681 INTEGER ILAENV 682 DOUBLE PRECISION DLAMCH, DLARND, DSXT1 683 EXTERNAL ILAENV, DLAMCH, DLARND, DSXT1 684* .. 685* .. External Subroutines .. 686 EXTERNAL DCOPY, DLABAD, DLASUM, DSTEBZ, DSTECH, DSTERF, 687 $ XERBLA, ZCOPY, ZHET21, ZHETRD, ZHPT21, ZHPTRD, 688 $ ZLACPY, ZLASET, ZLATMR, ZLATMS, ZPTEQR, ZSTEDC, 689 $ ZSTEMR, ZSTEIN, ZSTEQR, ZSTT21, ZSTT22, ZUNGTR, 690 $ ZUPGTR, ZHETRD_2STAGE 691* .. 692* .. Intrinsic Functions .. 693 INTRINSIC ABS, DBLE, DCONJG, INT, LOG, MAX, MIN, SQRT 694* .. 695* .. Data statements .. 696 DATA KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8, 697 $ 8, 8, 9, 9, 9, 9, 9, 10 / 698 DATA KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1, 699 $ 2, 3, 1, 1, 1, 2, 3, 1 / 700 DATA KMODE / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0, 701 $ 0, 0, 4, 3, 1, 4, 4, 3 / 702* .. 703* .. Executable Statements .. 704* 705* Keep ftnchek happy 706 IDUMMA( 1 ) = 1 707* 708* Check for errors 709* 710 NTESTT = 0 711 INFO = 0 712* 713* Important constants 714* 715 BADNN = .FALSE. 716 TRYRAC = .TRUE. 717 NMAX = 1 718 DO 10 J = 1, NSIZES 719 NMAX = MAX( NMAX, NN( J ) ) 720 IF( NN( J ).LT.0 ) 721 $ BADNN = .TRUE. 722 10 CONTINUE 723* 724 NBLOCK = ILAENV( 1, 'ZHETRD', 'L', NMAX, -1, -1, -1 ) 725 NBLOCK = MIN( NMAX, MAX( 1, NBLOCK ) ) 726* 727* Check for errors 728* 729 IF( NSIZES.LT.0 ) THEN 730 INFO = -1 731 ELSE IF( BADNN ) THEN 732 INFO = -2 733 ELSE IF( NTYPES.LT.0 ) THEN 734 INFO = -3 735 ELSE IF( LDA.LT.NMAX ) THEN 736 INFO = -9 737 ELSE IF( LDU.LT.NMAX ) THEN 738 INFO = -23 739 ELSE IF( 2*MAX( 2, NMAX )**2.GT.LWORK ) THEN 740 INFO = -29 741 END IF 742* 743 IF( INFO.NE.0 ) THEN 744 CALL XERBLA( 'ZCHKST2STG', -INFO ) 745 RETURN 746 END IF 747* 748* Quick return if possible 749* 750 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 751 $ RETURN 752* 753* More Important constants 754* 755 UNFL = DLAMCH( 'Safe minimum' ) 756 OVFL = ONE / UNFL 757 CALL DLABAD( UNFL, OVFL ) 758 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 759 ULPINV = ONE / ULP 760 LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) ) 761 RTUNFL = SQRT( UNFL ) 762 RTOVFL = SQRT( OVFL ) 763* 764* Loop over sizes, types 765* 766 DO 20 I = 1, 4 767 ISEED2( I ) = ISEED( I ) 768 20 CONTINUE 769 NERRS = 0 770 NMATS = 0 771* 772 DO 310 JSIZE = 1, NSIZES 773 N = NN( JSIZE ) 774 IF( N.GT.0 ) THEN 775 LGN = INT( LOG( DBLE( N ) ) / LOG( TWO ) ) 776 IF( 2**LGN.LT.N ) 777 $ LGN = LGN + 1 778 IF( 2**LGN.LT.N ) 779 $ LGN = LGN + 1 780 LWEDC = 1 + 4*N + 2*N*LGN + 4*N**2 781 LRWEDC = 1 + 3*N + 2*N*LGN + 4*N**2 782 LIWEDC = 6 + 6*N + 5*N*LGN 783 ELSE 784 LWEDC = 8 785 LRWEDC = 7 786 LIWEDC = 12 787 END IF 788 NAP = ( N*( N+1 ) ) / 2 789 ANINV = ONE / DBLE( MAX( 1, N ) ) 790* 791 IF( NSIZES.NE.1 ) THEN 792 MTYPES = MIN( MAXTYP, NTYPES ) 793 ELSE 794 MTYPES = MIN( MAXTYP+1, NTYPES ) 795 END IF 796* 797 DO 300 JTYPE = 1, MTYPES 798 IF( .NOT.DOTYPE( JTYPE ) ) 799 $ GO TO 300 800 NMATS = NMATS + 1 801 NTEST = 0 802* 803 DO 30 J = 1, 4 804 IOLDSD( J ) = ISEED( J ) 805 30 CONTINUE 806* 807* Compute "A" 808* 809* Control parameters: 810* 811* KMAGN KMODE KTYPE 812* =1 O(1) clustered 1 zero 813* =2 large clustered 2 identity 814* =3 small exponential (none) 815* =4 arithmetic diagonal, (w/ eigenvalues) 816* =5 random log Hermitian, w/ eigenvalues 817* =6 random (none) 818* =7 random diagonal 819* =8 random Hermitian 820* =9 positive definite 821* =10 diagonally dominant tridiagonal 822* 823 IF( MTYPES.GT.MAXTYP ) 824 $ GO TO 100 825* 826 ITYPE = KTYPE( JTYPE ) 827 IMODE = KMODE( JTYPE ) 828* 829* Compute norm 830* 831 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 832* 833 40 CONTINUE 834 ANORM = ONE 835 GO TO 70 836* 837 50 CONTINUE 838 ANORM = ( RTOVFL*ULP )*ANINV 839 GO TO 70 840* 841 60 CONTINUE 842 ANORM = RTUNFL*N*ULPINV 843 GO TO 70 844* 845 70 CONTINUE 846* 847 CALL ZLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA ) 848 IINFO = 0 849 IF( JTYPE.LE.15 ) THEN 850 COND = ULPINV 851 ELSE 852 COND = ULPINV*ANINV / TEN 853 END IF 854* 855* Special Matrices -- Identity & Jordan block 856* 857* Zero 858* 859 IF( ITYPE.EQ.1 ) THEN 860 IINFO = 0 861* 862 ELSE IF( ITYPE.EQ.2 ) THEN 863* 864* Identity 865* 866 DO 80 JC = 1, N 867 A( JC, JC ) = ANORM 868 80 CONTINUE 869* 870 ELSE IF( ITYPE.EQ.4 ) THEN 871* 872* Diagonal Matrix, [Eigen]values Specified 873* 874 CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND, 875 $ ANORM, 0, 0, 'N', A, LDA, WORK, IINFO ) 876* 877* 878 ELSE IF( ITYPE.EQ.5 ) THEN 879* 880* Hermitian, eigenvalues specified 881* 882 CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND, 883 $ ANORM, N, N, 'N', A, LDA, WORK, IINFO ) 884* 885 ELSE IF( ITYPE.EQ.7 ) THEN 886* 887* Diagonal, random eigenvalues 888* 889 CALL ZLATMR( N, N, 'S', ISEED, 'H', WORK, 6, ONE, CONE, 890 $ 'T', 'N', WORK( N+1 ), 1, ONE, 891 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 892 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 893* 894 ELSE IF( ITYPE.EQ.8 ) THEN 895* 896* Hermitian, random eigenvalues 897* 898 CALL ZLATMR( N, N, 'S', ISEED, 'H', WORK, 6, ONE, CONE, 899 $ 'T', 'N', WORK( N+1 ), 1, ONE, 900 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, 901 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 902* 903 ELSE IF( ITYPE.EQ.9 ) THEN 904* 905* Positive definite, eigenvalues specified. 906* 907 CALL ZLATMS( N, N, 'S', ISEED, 'P', RWORK, IMODE, COND, 908 $ ANORM, N, N, 'N', A, LDA, WORK, IINFO ) 909* 910 ELSE IF( ITYPE.EQ.10 ) THEN 911* 912* Positive definite tridiagonal, eigenvalues specified. 913* 914 CALL ZLATMS( N, N, 'S', ISEED, 'P', RWORK, IMODE, COND, 915 $ ANORM, 1, 1, 'N', A, LDA, WORK, IINFO ) 916 DO 90 I = 2, N 917 TEMP1 = ABS( A( I-1, I ) ) 918 TEMP2 = SQRT( ABS( A( I-1, I-1 )*A( I, I ) ) ) 919 IF( TEMP1.GT.HALF*TEMP2 ) THEN 920 A( I-1, I ) = A( I-1, I )* 921 $ ( HALF*TEMP2 / ( UNFL+TEMP1 ) ) 922 A( I, I-1 ) = DCONJG( A( I-1, I ) ) 923 END IF 924 90 CONTINUE 925* 926 ELSE 927* 928 IINFO = 1 929 END IF 930* 931 IF( IINFO.NE.0 ) THEN 932 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE, 933 $ IOLDSD 934 INFO = ABS( IINFO ) 935 RETURN 936 END IF 937* 938 100 CONTINUE 939* 940* Call ZHETRD and ZUNGTR to compute S and U from 941* upper triangle. 942* 943 CALL ZLACPY( 'U', N, N, A, LDA, V, LDU ) 944* 945 NTEST = 1 946 CALL ZHETRD( 'U', N, V, LDU, SD, SE, TAU, WORK, LWORK, 947 $ IINFO ) 948* 949 IF( IINFO.NE.0 ) THEN 950 WRITE( NOUNIT, FMT = 9999 )'ZHETRD(U)', IINFO, N, JTYPE, 951 $ IOLDSD 952 INFO = ABS( IINFO ) 953 IF( IINFO.LT.0 ) THEN 954 RETURN 955 ELSE 956 RESULT( 1 ) = ULPINV 957 GO TO 280 958 END IF 959 END IF 960* 961 CALL ZLACPY( 'U', N, N, V, LDU, U, LDU ) 962* 963 NTEST = 2 964 CALL ZUNGTR( 'U', N, U, LDU, TAU, WORK, LWORK, IINFO ) 965 IF( IINFO.NE.0 ) THEN 966 WRITE( NOUNIT, FMT = 9999 )'ZUNGTR(U)', IINFO, N, JTYPE, 967 $ IOLDSD 968 INFO = ABS( IINFO ) 969 IF( IINFO.LT.0 ) THEN 970 RETURN 971 ELSE 972 RESULT( 2 ) = ULPINV 973 GO TO 280 974 END IF 975 END IF 976* 977* Do tests 1 and 2 978* 979 CALL ZHET21( 2, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V, 980 $ LDU, TAU, WORK, RWORK, RESULT( 1 ) ) 981 CALL ZHET21( 3, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V, 982 $ LDU, TAU, WORK, RWORK, RESULT( 2 ) ) 983* 984* Compute D1 the eigenvalues resulting from the tridiagonal 985* form using the standard 1-stage algorithm and use it as a 986* reference to compare with the 2-stage technique 987* 988* Compute D1 from the 1-stage and used as reference for the 989* 2-stage 990* 991 CALL DCOPY( N, SD, 1, D1, 1 ) 992 IF( N.GT.0 ) 993 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 994* 995 CALL ZSTEQR( 'N', N, D1, RWORK, WORK, LDU, RWORK( N+1 ), 996 $ IINFO ) 997 IF( IINFO.NE.0 ) THEN 998 WRITE( NOUNIT, FMT = 9999 )'ZSTEQR(N)', IINFO, N, JTYPE, 999 $ IOLDSD 1000 INFO = ABS( IINFO ) 1001 IF( IINFO.LT.0 ) THEN 1002 RETURN 1003 ELSE 1004 RESULT( 3 ) = ULPINV 1005 GO TO 280 1006 END IF 1007 END IF 1008* 1009* 2-STAGE TRD Upper case is used to compute D2. 1010* Note to set SD and SE to zero to be sure not reusing 1011* the one from above. Compare it with D1 computed 1012* using the 1-stage. 1013* 1014 CALL DLASET( 'Full', N, 1, ZERO, ZERO, SD, N ) 1015 CALL DLASET( 'Full', N, 1, ZERO, ZERO, SE, N ) 1016 CALL ZLACPY( 'U', N, N, A, LDA, V, LDU ) 1017 LH = MAX(1, 4*N) 1018 LW = LWORK - LH 1019 CALL ZHETRD_2STAGE( 'N', "U", N, V, LDU, SD, SE, TAU, 1020 $ WORK, LH, WORK( LH+1 ), LW, IINFO ) 1021* 1022* Compute D2 from the 2-stage Upper case 1023* 1024 CALL DCOPY( N, SD, 1, D2, 1 ) 1025 IF( N.GT.0 ) 1026 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1027* 1028 NTEST = 3 1029 CALL ZSTEQR( 'N', N, D2, RWORK, WORK, LDU, RWORK( N+1 ), 1030 $ IINFO ) 1031 IF( IINFO.NE.0 ) THEN 1032 WRITE( NOUNIT, FMT = 9999 )'ZSTEQR(N)', IINFO, N, JTYPE, 1033 $ IOLDSD 1034 INFO = ABS( IINFO ) 1035 IF( IINFO.LT.0 ) THEN 1036 RETURN 1037 ELSE 1038 RESULT( 3 ) = ULPINV 1039 GO TO 280 1040 END IF 1041 END IF 1042* 1043* 2-STAGE TRD Lower case is used to compute D3. 1044* Note to set SD and SE to zero to be sure not reusing 1045* the one from above. Compare it with D1 computed 1046* using the 1-stage. 1047* 1048 CALL DLASET( 'Full', N, 1, ZERO, ZERO, SD, N ) 1049 CALL DLASET( 'Full', N, 1, ZERO, ZERO, SE, N ) 1050 CALL ZLACPY( 'L', N, N, A, LDA, V, LDU ) 1051 CALL ZHETRD_2STAGE( 'N', "L", N, V, LDU, SD, SE, TAU, 1052 $ WORK, LH, WORK( LH+1 ), LW, IINFO ) 1053* 1054* Compute D3 from the 2-stage Upper case 1055* 1056 CALL DCOPY( N, SD, 1, D3, 1 ) 1057 IF( N.GT.0 ) 1058 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1059* 1060 NTEST = 4 1061 CALL ZSTEQR( 'N', N, D3, RWORK, WORK, LDU, RWORK( N+1 ), 1062 $ IINFO ) 1063 IF( IINFO.NE.0 ) THEN 1064 WRITE( NOUNIT, FMT = 9999 )'ZSTEQR(N)', IINFO, N, JTYPE, 1065 $ IOLDSD 1066 INFO = ABS( IINFO ) 1067 IF( IINFO.LT.0 ) THEN 1068 RETURN 1069 ELSE 1070 RESULT( 4 ) = ULPINV 1071 GO TO 280 1072 END IF 1073 END IF 1074* 1075* 1076* Do Tests 3 and 4 which are similar to 11 and 12 but with the 1077* D1 computed using the standard 1-stage reduction as reference 1078* 1079 NTEST = 4 1080 TEMP1 = ZERO 1081 TEMP2 = ZERO 1082 TEMP3 = ZERO 1083 TEMP4 = ZERO 1084* 1085 DO 151 J = 1, N 1086 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) ) 1087 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 1088 TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) ) 1089 TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) ) 1090 151 CONTINUE 1091* 1092 RESULT( 3 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 1093 RESULT( 4 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) ) 1094* 1095* Store the upper triangle of A in AP 1096* 1097 I = 0 1098 DO 120 JC = 1, N 1099 DO 110 JR = 1, JC 1100 I = I + 1 1101 AP( I ) = A( JR, JC ) 1102 110 CONTINUE 1103 120 CONTINUE 1104* 1105* Call ZHPTRD and ZUPGTR to compute S and U from AP 1106* 1107 CALL ZCOPY( NAP, AP, 1, VP, 1 ) 1108* 1109 NTEST = 5 1110 CALL ZHPTRD( 'U', N, VP, SD, SE, TAU, IINFO ) 1111* 1112 IF( IINFO.NE.0 ) THEN 1113 WRITE( NOUNIT, FMT = 9999 )'ZHPTRD(U)', IINFO, N, JTYPE, 1114 $ IOLDSD 1115 INFO = ABS( IINFO ) 1116 IF( IINFO.LT.0 ) THEN 1117 RETURN 1118 ELSE 1119 RESULT( 5 ) = ULPINV 1120 GO TO 280 1121 END IF 1122 END IF 1123* 1124 NTEST = 6 1125 CALL ZUPGTR( 'U', N, VP, TAU, U, LDU, WORK, IINFO ) 1126 IF( IINFO.NE.0 ) THEN 1127 WRITE( NOUNIT, FMT = 9999 )'ZUPGTR(U)', IINFO, N, JTYPE, 1128 $ IOLDSD 1129 INFO = ABS( IINFO ) 1130 IF( IINFO.LT.0 ) THEN 1131 RETURN 1132 ELSE 1133 RESULT( 6 ) = ULPINV 1134 GO TO 280 1135 END IF 1136 END IF 1137* 1138* Do tests 5 and 6 1139* 1140 CALL ZHPT21( 2, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU, 1141 $ WORK, RWORK, RESULT( 5 ) ) 1142 CALL ZHPT21( 3, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU, 1143 $ WORK, RWORK, RESULT( 6 ) ) 1144* 1145* Store the lower triangle of A in AP 1146* 1147 I = 0 1148 DO 140 JC = 1, N 1149 DO 130 JR = JC, N 1150 I = I + 1 1151 AP( I ) = A( JR, JC ) 1152 130 CONTINUE 1153 140 CONTINUE 1154* 1155* Call ZHPTRD and ZUPGTR to compute S and U from AP 1156* 1157 CALL ZCOPY( NAP, AP, 1, VP, 1 ) 1158* 1159 NTEST = 7 1160 CALL ZHPTRD( 'L', N, VP, SD, SE, TAU, IINFO ) 1161* 1162 IF( IINFO.NE.0 ) THEN 1163 WRITE( NOUNIT, FMT = 9999 )'ZHPTRD(L)', IINFO, N, JTYPE, 1164 $ IOLDSD 1165 INFO = ABS( IINFO ) 1166 IF( IINFO.LT.0 ) THEN 1167 RETURN 1168 ELSE 1169 RESULT( 7 ) = ULPINV 1170 GO TO 280 1171 END IF 1172 END IF 1173* 1174 NTEST = 8 1175 CALL ZUPGTR( 'L', N, VP, TAU, U, LDU, WORK, IINFO ) 1176 IF( IINFO.NE.0 ) THEN 1177 WRITE( NOUNIT, FMT = 9999 )'ZUPGTR(L)', IINFO, N, JTYPE, 1178 $ IOLDSD 1179 INFO = ABS( IINFO ) 1180 IF( IINFO.LT.0 ) THEN 1181 RETURN 1182 ELSE 1183 RESULT( 8 ) = ULPINV 1184 GO TO 280 1185 END IF 1186 END IF 1187* 1188 CALL ZHPT21( 2, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU, 1189 $ WORK, RWORK, RESULT( 7 ) ) 1190 CALL ZHPT21( 3, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU, 1191 $ WORK, RWORK, RESULT( 8 ) ) 1192* 1193* Call ZSTEQR to compute D1, D2, and Z, do tests. 1194* 1195* Compute D1 and Z 1196* 1197 CALL DCOPY( N, SD, 1, D1, 1 ) 1198 IF( N.GT.0 ) 1199 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1200 CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) 1201* 1202 NTEST = 9 1203 CALL ZSTEQR( 'V', N, D1, RWORK, Z, LDU, RWORK( N+1 ), 1204 $ IINFO ) 1205 IF( IINFO.NE.0 ) THEN 1206 WRITE( NOUNIT, FMT = 9999 )'ZSTEQR(V)', IINFO, N, JTYPE, 1207 $ IOLDSD 1208 INFO = ABS( IINFO ) 1209 IF( IINFO.LT.0 ) THEN 1210 RETURN 1211 ELSE 1212 RESULT( 9 ) = ULPINV 1213 GO TO 280 1214 END IF 1215 END IF 1216* 1217* Compute D2 1218* 1219 CALL DCOPY( N, SD, 1, D2, 1 ) 1220 IF( N.GT.0 ) 1221 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1222* 1223 NTEST = 11 1224 CALL ZSTEQR( 'N', N, D2, RWORK, WORK, LDU, RWORK( N+1 ), 1225 $ IINFO ) 1226 IF( IINFO.NE.0 ) THEN 1227 WRITE( NOUNIT, FMT = 9999 )'ZSTEQR(N)', IINFO, N, JTYPE, 1228 $ IOLDSD 1229 INFO = ABS( IINFO ) 1230 IF( IINFO.LT.0 ) THEN 1231 RETURN 1232 ELSE 1233 RESULT( 11 ) = ULPINV 1234 GO TO 280 1235 END IF 1236 END IF 1237* 1238* Compute D3 (using PWK method) 1239* 1240 CALL DCOPY( N, SD, 1, D3, 1 ) 1241 IF( N.GT.0 ) 1242 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1243* 1244 NTEST = 12 1245 CALL DSTERF( N, D3, RWORK, IINFO ) 1246 IF( IINFO.NE.0 ) THEN 1247 WRITE( NOUNIT, FMT = 9999 )'DSTERF', IINFO, N, JTYPE, 1248 $ IOLDSD 1249 INFO = ABS( IINFO ) 1250 IF( IINFO.LT.0 ) THEN 1251 RETURN 1252 ELSE 1253 RESULT( 12 ) = ULPINV 1254 GO TO 280 1255 END IF 1256 END IF 1257* 1258* Do Tests 9 and 10 1259* 1260 CALL ZSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK, 1261 $ RESULT( 9 ) ) 1262* 1263* Do Tests 11 and 12 1264* 1265 TEMP1 = ZERO 1266 TEMP2 = ZERO 1267 TEMP3 = ZERO 1268 TEMP4 = ZERO 1269* 1270 DO 150 J = 1, N 1271 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) ) 1272 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 1273 TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) ) 1274 TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) ) 1275 150 CONTINUE 1276* 1277 RESULT( 11 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 1278 RESULT( 12 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) ) 1279* 1280* Do Test 13 -- Sturm Sequence Test of Eigenvalues 1281* Go up by factors of two until it succeeds 1282* 1283 NTEST = 13 1284 TEMP1 = THRESH*( HALF-ULP ) 1285* 1286 DO 160 J = 0, LOG2UI 1287 CALL DSTECH( N, SD, SE, D1, TEMP1, RWORK, IINFO ) 1288 IF( IINFO.EQ.0 ) 1289 $ GO TO 170 1290 TEMP1 = TEMP1*TWO 1291 160 CONTINUE 1292* 1293 170 CONTINUE 1294 RESULT( 13 ) = TEMP1 1295* 1296* For positive definite matrices ( JTYPE.GT.15 ) call ZPTEQR 1297* and do tests 14, 15, and 16 . 1298* 1299 IF( JTYPE.GT.15 ) THEN 1300* 1301* Compute D4 and Z4 1302* 1303 CALL DCOPY( N, SD, 1, D4, 1 ) 1304 IF( N.GT.0 ) 1305 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1306 CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) 1307* 1308 NTEST = 14 1309 CALL ZPTEQR( 'V', N, D4, RWORK, Z, LDU, RWORK( N+1 ), 1310 $ IINFO ) 1311 IF( IINFO.NE.0 ) THEN 1312 WRITE( NOUNIT, FMT = 9999 )'ZPTEQR(V)', IINFO, N, 1313 $ JTYPE, IOLDSD 1314 INFO = ABS( IINFO ) 1315 IF( IINFO.LT.0 ) THEN 1316 RETURN 1317 ELSE 1318 RESULT( 14 ) = ULPINV 1319 GO TO 280 1320 END IF 1321 END IF 1322* 1323* Do Tests 14 and 15 1324* 1325 CALL ZSTT21( N, 0, SD, SE, D4, DUMMA, Z, LDU, WORK, 1326 $ RWORK, RESULT( 14 ) ) 1327* 1328* Compute D5 1329* 1330 CALL DCOPY( N, SD, 1, D5, 1 ) 1331 IF( N.GT.0 ) 1332 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1333* 1334 NTEST = 16 1335 CALL ZPTEQR( 'N', N, D5, RWORK, Z, LDU, RWORK( N+1 ), 1336 $ IINFO ) 1337 IF( IINFO.NE.0 ) THEN 1338 WRITE( NOUNIT, FMT = 9999 )'ZPTEQR(N)', IINFO, N, 1339 $ JTYPE, IOLDSD 1340 INFO = ABS( IINFO ) 1341 IF( IINFO.LT.0 ) THEN 1342 RETURN 1343 ELSE 1344 RESULT( 16 ) = ULPINV 1345 GO TO 280 1346 END IF 1347 END IF 1348* 1349* Do Test 16 1350* 1351 TEMP1 = ZERO 1352 TEMP2 = ZERO 1353 DO 180 J = 1, N 1354 TEMP1 = MAX( TEMP1, ABS( D4( J ) ), ABS( D5( J ) ) ) 1355 TEMP2 = MAX( TEMP2, ABS( D4( J )-D5( J ) ) ) 1356 180 CONTINUE 1357* 1358 RESULT( 16 ) = TEMP2 / MAX( UNFL, 1359 $ HUN*ULP*MAX( TEMP1, TEMP2 ) ) 1360 ELSE 1361 RESULT( 14 ) = ZERO 1362 RESULT( 15 ) = ZERO 1363 RESULT( 16 ) = ZERO 1364 END IF 1365* 1366* Call DSTEBZ with different options and do tests 17-18. 1367* 1368* If S is positive definite and diagonally dominant, 1369* ask for all eigenvalues with high relative accuracy. 1370* 1371 VL = ZERO 1372 VU = ZERO 1373 IL = 0 1374 IU = 0 1375 IF( JTYPE.EQ.21 ) THEN 1376 NTEST = 17 1377 ABSTOL = UNFL + UNFL 1378 CALL DSTEBZ( 'A', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, 1379 $ M, NSPLIT, WR, IWORK( 1 ), IWORK( N+1 ), 1380 $ RWORK, IWORK( 2*N+1 ), IINFO ) 1381 IF( IINFO.NE.0 ) THEN 1382 WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A,rel)', IINFO, N, 1383 $ JTYPE, IOLDSD 1384 INFO = ABS( IINFO ) 1385 IF( IINFO.LT.0 ) THEN 1386 RETURN 1387 ELSE 1388 RESULT( 17 ) = ULPINV 1389 GO TO 280 1390 END IF 1391 END IF 1392* 1393* Do test 17 1394* 1395 TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) / 1396 $ ( ONE-HALF )**4 1397* 1398 TEMP1 = ZERO 1399 DO 190 J = 1, N 1400 TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) / 1401 $ ( ABSTOL+ABS( D4( J ) ) ) ) 1402 190 CONTINUE 1403* 1404 RESULT( 17 ) = TEMP1 / TEMP2 1405 ELSE 1406 RESULT( 17 ) = ZERO 1407 END IF 1408* 1409* Now ask for all eigenvalues with high absolute accuracy. 1410* 1411 NTEST = 18 1412 ABSTOL = UNFL + UNFL 1413 CALL DSTEBZ( 'A', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, M, 1414 $ NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), RWORK, 1415 $ IWORK( 2*N+1 ), IINFO ) 1416 IF( IINFO.NE.0 ) THEN 1417 WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A)', IINFO, N, JTYPE, 1418 $ IOLDSD 1419 INFO = ABS( IINFO ) 1420 IF( IINFO.LT.0 ) THEN 1421 RETURN 1422 ELSE 1423 RESULT( 18 ) = ULPINV 1424 GO TO 280 1425 END IF 1426 END IF 1427* 1428* Do test 18 1429* 1430 TEMP1 = ZERO 1431 TEMP2 = ZERO 1432 DO 200 J = 1, N 1433 TEMP1 = MAX( TEMP1, ABS( D3( J ) ), ABS( WA1( J ) ) ) 1434 TEMP2 = MAX( TEMP2, ABS( D3( J )-WA1( J ) ) ) 1435 200 CONTINUE 1436* 1437 RESULT( 18 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 1438* 1439* Choose random values for IL and IU, and ask for the 1440* IL-th through IU-th eigenvalues. 1441* 1442 NTEST = 19 1443 IF( N.LE.1 ) THEN 1444 IL = 1 1445 IU = N 1446 ELSE 1447 IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) 1448 IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) 1449 IF( IU.LT.IL ) THEN 1450 ITEMP = IU 1451 IU = IL 1452 IL = ITEMP 1453 END IF 1454 END IF 1455* 1456 CALL DSTEBZ( 'I', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, 1457 $ M2, NSPLIT, WA2, IWORK( 1 ), IWORK( N+1 ), 1458 $ RWORK, IWORK( 2*N+1 ), IINFO ) 1459 IF( IINFO.NE.0 ) THEN 1460 WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(I)', IINFO, N, JTYPE, 1461 $ IOLDSD 1462 INFO = ABS( IINFO ) 1463 IF( IINFO.LT.0 ) THEN 1464 RETURN 1465 ELSE 1466 RESULT( 19 ) = ULPINV 1467 GO TO 280 1468 END IF 1469 END IF 1470* 1471* Determine the values VL and VU of the IL-th and IU-th 1472* eigenvalues and ask for all eigenvalues in this range. 1473* 1474 IF( N.GT.0 ) THEN 1475 IF( IL.NE.1 ) THEN 1476 VL = WA1( IL ) - MAX( HALF*( WA1( IL )-WA1( IL-1 ) ), 1477 $ ULP*ANORM, TWO*RTUNFL ) 1478 ELSE 1479 VL = WA1( 1 ) - MAX( HALF*( WA1( N )-WA1( 1 ) ), 1480 $ ULP*ANORM, TWO*RTUNFL ) 1481 END IF 1482 IF( IU.NE.N ) THEN 1483 VU = WA1( IU ) + MAX( HALF*( WA1( IU+1 )-WA1( IU ) ), 1484 $ ULP*ANORM, TWO*RTUNFL ) 1485 ELSE 1486 VU = WA1( N ) + MAX( HALF*( WA1( N )-WA1( 1 ) ), 1487 $ ULP*ANORM, TWO*RTUNFL ) 1488 END IF 1489 ELSE 1490 VL = ZERO 1491 VU = ONE 1492 END IF 1493* 1494 CALL DSTEBZ( 'V', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, 1495 $ M3, NSPLIT, WA3, IWORK( 1 ), IWORK( N+1 ), 1496 $ RWORK, IWORK( 2*N+1 ), IINFO ) 1497 IF( IINFO.NE.0 ) THEN 1498 WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(V)', IINFO, N, JTYPE, 1499 $ IOLDSD 1500 INFO = ABS( IINFO ) 1501 IF( IINFO.LT.0 ) THEN 1502 RETURN 1503 ELSE 1504 RESULT( 19 ) = ULPINV 1505 GO TO 280 1506 END IF 1507 END IF 1508* 1509 IF( M3.EQ.0 .AND. N.NE.0 ) THEN 1510 RESULT( 19 ) = ULPINV 1511 GO TO 280 1512 END IF 1513* 1514* Do test 19 1515* 1516 TEMP1 = DSXT1( 1, WA2, M2, WA3, M3, ABSTOL, ULP, UNFL ) 1517 TEMP2 = DSXT1( 1, WA3, M3, WA2, M2, ABSTOL, ULP, UNFL ) 1518 IF( N.GT.0 ) THEN 1519 TEMP3 = MAX( ABS( WA1( N ) ), ABS( WA1( 1 ) ) ) 1520 ELSE 1521 TEMP3 = ZERO 1522 END IF 1523* 1524 RESULT( 19 ) = ( TEMP1+TEMP2 ) / MAX( UNFL, TEMP3*ULP ) 1525* 1526* Call ZSTEIN to compute eigenvectors corresponding to 1527* eigenvalues in WA1. (First call DSTEBZ again, to make sure 1528* it returns these eigenvalues in the correct order.) 1529* 1530 NTEST = 21 1531 CALL DSTEBZ( 'A', 'B', N, VL, VU, IL, IU, ABSTOL, SD, SE, M, 1532 $ NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), RWORK, 1533 $ IWORK( 2*N+1 ), IINFO ) 1534 IF( IINFO.NE.0 ) THEN 1535 WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A,B)', IINFO, N, 1536 $ JTYPE, IOLDSD 1537 INFO = ABS( IINFO ) 1538 IF( IINFO.LT.0 ) THEN 1539 RETURN 1540 ELSE 1541 RESULT( 20 ) = ULPINV 1542 RESULT( 21 ) = ULPINV 1543 GO TO 280 1544 END IF 1545 END IF 1546* 1547 CALL ZSTEIN( N, SD, SE, M, WA1, IWORK( 1 ), IWORK( N+1 ), Z, 1548 $ LDU, RWORK, IWORK( 2*N+1 ), IWORK( 3*N+1 ), 1549 $ IINFO ) 1550 IF( IINFO.NE.0 ) THEN 1551 WRITE( NOUNIT, FMT = 9999 )'ZSTEIN', IINFO, N, JTYPE, 1552 $ IOLDSD 1553 INFO = ABS( IINFO ) 1554 IF( IINFO.LT.0 ) THEN 1555 RETURN 1556 ELSE 1557 RESULT( 20 ) = ULPINV 1558 RESULT( 21 ) = ULPINV 1559 GO TO 280 1560 END IF 1561 END IF 1562* 1563* Do tests 20 and 21 1564* 1565 CALL ZSTT21( N, 0, SD, SE, WA1, DUMMA, Z, LDU, WORK, RWORK, 1566 $ RESULT( 20 ) ) 1567* 1568* Call ZSTEDC(I) to compute D1 and Z, do tests. 1569* 1570* Compute D1 and Z 1571* 1572 INDE = 1 1573 INDRWK = INDE + N 1574 CALL DCOPY( N, SD, 1, D1, 1 ) 1575 IF( N.GT.0 ) 1576 $ CALL DCOPY( N-1, SE, 1, RWORK( INDE ), 1 ) 1577 CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) 1578* 1579 NTEST = 22 1580 CALL ZSTEDC( 'I', N, D1, RWORK( INDE ), Z, LDU, WORK, LWEDC, 1581 $ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO ) 1582 IF( IINFO.NE.0 ) THEN 1583 WRITE( NOUNIT, FMT = 9999 )'ZSTEDC(I)', IINFO, N, JTYPE, 1584 $ IOLDSD 1585 INFO = ABS( IINFO ) 1586 IF( IINFO.LT.0 ) THEN 1587 RETURN 1588 ELSE 1589 RESULT( 22 ) = ULPINV 1590 GO TO 280 1591 END IF 1592 END IF 1593* 1594* Do Tests 22 and 23 1595* 1596 CALL ZSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK, 1597 $ RESULT( 22 ) ) 1598* 1599* Call ZSTEDC(V) to compute D1 and Z, do tests. 1600* 1601* Compute D1 and Z 1602* 1603 CALL DCOPY( N, SD, 1, D1, 1 ) 1604 IF( N.GT.0 ) 1605 $ CALL DCOPY( N-1, SE, 1, RWORK( INDE ), 1 ) 1606 CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) 1607* 1608 NTEST = 24 1609 CALL ZSTEDC( 'V', N, D1, RWORK( INDE ), Z, LDU, WORK, LWEDC, 1610 $ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO ) 1611 IF( IINFO.NE.0 ) THEN 1612 WRITE( NOUNIT, FMT = 9999 )'ZSTEDC(V)', IINFO, N, JTYPE, 1613 $ IOLDSD 1614 INFO = ABS( IINFO ) 1615 IF( IINFO.LT.0 ) THEN 1616 RETURN 1617 ELSE 1618 RESULT( 24 ) = ULPINV 1619 GO TO 280 1620 END IF 1621 END IF 1622* 1623* Do Tests 24 and 25 1624* 1625 CALL ZSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK, 1626 $ RESULT( 24 ) ) 1627* 1628* Call ZSTEDC(N) to compute D2, do tests. 1629* 1630* Compute D2 1631* 1632 CALL DCOPY( N, SD, 1, D2, 1 ) 1633 IF( N.GT.0 ) 1634 $ CALL DCOPY( N-1, SE, 1, RWORK( INDE ), 1 ) 1635 CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) 1636* 1637 NTEST = 26 1638 CALL ZSTEDC( 'N', N, D2, RWORK( INDE ), Z, LDU, WORK, LWEDC, 1639 $ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO ) 1640 IF( IINFO.NE.0 ) THEN 1641 WRITE( NOUNIT, FMT = 9999 )'ZSTEDC(N)', IINFO, N, JTYPE, 1642 $ IOLDSD 1643 INFO = ABS( IINFO ) 1644 IF( IINFO.LT.0 ) THEN 1645 RETURN 1646 ELSE 1647 RESULT( 26 ) = ULPINV 1648 GO TO 280 1649 END IF 1650 END IF 1651* 1652* Do Test 26 1653* 1654 TEMP1 = ZERO 1655 TEMP2 = ZERO 1656* 1657 DO 210 J = 1, N 1658 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) ) 1659 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 1660 210 CONTINUE 1661* 1662 RESULT( 26 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 1663* 1664* Only test ZSTEMR if IEEE compliant 1665* 1666 IF( ILAENV( 10, 'ZSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND. 1667 $ ILAENV( 11, 'ZSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN 1668* 1669* Call ZSTEMR, do test 27 (relative eigenvalue accuracy) 1670* 1671* If S is positive definite and diagonally dominant, 1672* ask for all eigenvalues with high relative accuracy. 1673* 1674 VL = ZERO 1675 VU = ZERO 1676 IL = 0 1677 IU = 0 1678 IF( JTYPE.EQ.21 .AND. CREL ) THEN 1679 NTEST = 27 1680 ABSTOL = UNFL + UNFL 1681 CALL ZSTEMR( 'V', 'A', N, SD, SE, VL, VU, IL, IU, 1682 $ M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC, 1683 $ RWORK, LRWORK, IWORK( 2*N+1 ), LWORK-2*N, 1684 $ IINFO ) 1685 IF( IINFO.NE.0 ) THEN 1686 WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(V,A,rel)', 1687 $ IINFO, N, JTYPE, IOLDSD 1688 INFO = ABS( IINFO ) 1689 IF( IINFO.LT.0 ) THEN 1690 RETURN 1691 ELSE 1692 RESULT( 27 ) = ULPINV 1693 GO TO 270 1694 END IF 1695 END IF 1696* 1697* Do test 27 1698* 1699 TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) / 1700 $ ( ONE-HALF )**4 1701* 1702 TEMP1 = ZERO 1703 DO 220 J = 1, N 1704 TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) / 1705 $ ( ABSTOL+ABS( D4( J ) ) ) ) 1706 220 CONTINUE 1707* 1708 RESULT( 27 ) = TEMP1 / TEMP2 1709* 1710 IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) 1711 IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) 1712 IF( IU.LT.IL ) THEN 1713 ITEMP = IU 1714 IU = IL 1715 IL = ITEMP 1716 END IF 1717* 1718 IF( CRANGE ) THEN 1719 NTEST = 28 1720 ABSTOL = UNFL + UNFL 1721 CALL ZSTEMR( 'V', 'I', N, SD, SE, VL, VU, IL, IU, 1722 $ M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC, 1723 $ RWORK, LRWORK, IWORK( 2*N+1 ), 1724 $ LWORK-2*N, IINFO ) 1725* 1726 IF( IINFO.NE.0 ) THEN 1727 WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(V,I,rel)', 1728 $ IINFO, N, JTYPE, IOLDSD 1729 INFO = ABS( IINFO ) 1730 IF( IINFO.LT.0 ) THEN 1731 RETURN 1732 ELSE 1733 RESULT( 28 ) = ULPINV 1734 GO TO 270 1735 END IF 1736 END IF 1737* 1738* 1739* Do test 28 1740* 1741 TEMP2 = TWO*( TWO*N-ONE )*ULP* 1742 $ ( ONE+EIGHT*HALF**2 ) / ( ONE-HALF )**4 1743* 1744 TEMP1 = ZERO 1745 DO 230 J = IL, IU 1746 TEMP1 = MAX( TEMP1, ABS( WR( J-IL+1 )-D4( N-J+ 1747 $ 1 ) ) / ( ABSTOL+ABS( WR( J-IL+1 ) ) ) ) 1748 230 CONTINUE 1749* 1750 RESULT( 28 ) = TEMP1 / TEMP2 1751 ELSE 1752 RESULT( 28 ) = ZERO 1753 END IF 1754 ELSE 1755 RESULT( 27 ) = ZERO 1756 RESULT( 28 ) = ZERO 1757 END IF 1758* 1759* Call ZSTEMR(V,I) to compute D1 and Z, do tests. 1760* 1761* Compute D1 and Z 1762* 1763 CALL DCOPY( N, SD, 1, D5, 1 ) 1764 IF( N.GT.0 ) 1765 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1766 CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) 1767* 1768 IF( CRANGE ) THEN 1769 NTEST = 29 1770 IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) 1771 IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) 1772 IF( IU.LT.IL ) THEN 1773 ITEMP = IU 1774 IU = IL 1775 IL = ITEMP 1776 END IF 1777 CALL ZSTEMR( 'V', 'I', N, D5, RWORK, VL, VU, IL, IU, 1778 $ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC, 1779 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), 1780 $ LIWORK-2*N, IINFO ) 1781 IF( IINFO.NE.0 ) THEN 1782 WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(V,I)', IINFO, 1783 $ N, JTYPE, IOLDSD 1784 INFO = ABS( IINFO ) 1785 IF( IINFO.LT.0 ) THEN 1786 RETURN 1787 ELSE 1788 RESULT( 29 ) = ULPINV 1789 GO TO 280 1790 END IF 1791 END IF 1792* 1793* Do Tests 29 and 30 1794* 1795* 1796* Call ZSTEMR to compute D2, do tests. 1797* 1798* Compute D2 1799* 1800 CALL DCOPY( N, SD, 1, D5, 1 ) 1801 IF( N.GT.0 ) 1802 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1803* 1804 NTEST = 31 1805 CALL ZSTEMR( 'N', 'I', N, D5, RWORK, VL, VU, IL, IU, 1806 $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC, 1807 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), 1808 $ LIWORK-2*N, IINFO ) 1809 IF( IINFO.NE.0 ) THEN 1810 WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(N,I)', IINFO, 1811 $ N, JTYPE, IOLDSD 1812 INFO = ABS( IINFO ) 1813 IF( IINFO.LT.0 ) THEN 1814 RETURN 1815 ELSE 1816 RESULT( 31 ) = ULPINV 1817 GO TO 280 1818 END IF 1819 END IF 1820* 1821* Do Test 31 1822* 1823 TEMP1 = ZERO 1824 TEMP2 = ZERO 1825* 1826 DO 240 J = 1, IU - IL + 1 1827 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), 1828 $ ABS( D2( J ) ) ) 1829 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 1830 240 CONTINUE 1831* 1832 RESULT( 31 ) = TEMP2 / MAX( UNFL, 1833 $ ULP*MAX( TEMP1, TEMP2 ) ) 1834* 1835* 1836* Call ZSTEMR(V,V) to compute D1 and Z, do tests. 1837* 1838* Compute D1 and Z 1839* 1840 CALL DCOPY( N, SD, 1, D5, 1 ) 1841 IF( N.GT.0 ) 1842 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1843 CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) 1844* 1845 NTEST = 32 1846* 1847 IF( N.GT.0 ) THEN 1848 IF( IL.NE.1 ) THEN 1849 VL = D2( IL ) - MAX( HALF* 1850 $ ( D2( IL )-D2( IL-1 ) ), ULP*ANORM, 1851 $ TWO*RTUNFL ) 1852 ELSE 1853 VL = D2( 1 ) - MAX( HALF*( D2( N )-D2( 1 ) ), 1854 $ ULP*ANORM, TWO*RTUNFL ) 1855 END IF 1856 IF( IU.NE.N ) THEN 1857 VU = D2( IU ) + MAX( HALF* 1858 $ ( D2( IU+1 )-D2( IU ) ), ULP*ANORM, 1859 $ TWO*RTUNFL ) 1860 ELSE 1861 VU = D2( N ) + MAX( HALF*( D2( N )-D2( 1 ) ), 1862 $ ULP*ANORM, TWO*RTUNFL ) 1863 END IF 1864 ELSE 1865 VL = ZERO 1866 VU = ONE 1867 END IF 1868* 1869 CALL ZSTEMR( 'V', 'V', N, D5, RWORK, VL, VU, IL, IU, 1870 $ M, D1, Z, LDU, M, IWORK( 1 ), TRYRAC, 1871 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), 1872 $ LIWORK-2*N, IINFO ) 1873 IF( IINFO.NE.0 ) THEN 1874 WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(V,V)', IINFO, 1875 $ N, JTYPE, IOLDSD 1876 INFO = ABS( IINFO ) 1877 IF( IINFO.LT.0 ) THEN 1878 RETURN 1879 ELSE 1880 RESULT( 32 ) = ULPINV 1881 GO TO 280 1882 END IF 1883 END IF 1884* 1885* Do Tests 32 and 33 1886* 1887 CALL ZSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, 1888 $ M, RWORK, RESULT( 32 ) ) 1889* 1890* Call ZSTEMR to compute D2, do tests. 1891* 1892* Compute D2 1893* 1894 CALL DCOPY( N, SD, 1, D5, 1 ) 1895 IF( N.GT.0 ) 1896 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1897* 1898 NTEST = 34 1899 CALL ZSTEMR( 'N', 'V', N, D5, RWORK, VL, VU, IL, IU, 1900 $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC, 1901 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), 1902 $ LIWORK-2*N, IINFO ) 1903 IF( IINFO.NE.0 ) THEN 1904 WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(N,V)', IINFO, 1905 $ N, JTYPE, IOLDSD 1906 INFO = ABS( IINFO ) 1907 IF( IINFO.LT.0 ) THEN 1908 RETURN 1909 ELSE 1910 RESULT( 34 ) = ULPINV 1911 GO TO 280 1912 END IF 1913 END IF 1914* 1915* Do Test 34 1916* 1917 TEMP1 = ZERO 1918 TEMP2 = ZERO 1919* 1920 DO 250 J = 1, IU - IL + 1 1921 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), 1922 $ ABS( D2( J ) ) ) 1923 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 1924 250 CONTINUE 1925* 1926 RESULT( 34 ) = TEMP2 / MAX( UNFL, 1927 $ ULP*MAX( TEMP1, TEMP2 ) ) 1928 ELSE 1929 RESULT( 29 ) = ZERO 1930 RESULT( 30 ) = ZERO 1931 RESULT( 31 ) = ZERO 1932 RESULT( 32 ) = ZERO 1933 RESULT( 33 ) = ZERO 1934 RESULT( 34 ) = ZERO 1935 END IF 1936* 1937* 1938* Call ZSTEMR(V,A) to compute D1 and Z, do tests. 1939* 1940* Compute D1 and Z 1941* 1942 CALL DCOPY( N, SD, 1, D5, 1 ) 1943 IF( N.GT.0 ) 1944 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1945* 1946 NTEST = 35 1947* 1948 CALL ZSTEMR( 'V', 'A', N, D5, RWORK, VL, VU, IL, IU, 1949 $ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC, 1950 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), 1951 $ LIWORK-2*N, IINFO ) 1952 IF( IINFO.NE.0 ) THEN 1953 WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(V,A)', IINFO, N, 1954 $ JTYPE, IOLDSD 1955 INFO = ABS( IINFO ) 1956 IF( IINFO.LT.0 ) THEN 1957 RETURN 1958 ELSE 1959 RESULT( 35 ) = ULPINV 1960 GO TO 280 1961 END IF 1962 END IF 1963* 1964* Do Tests 35 and 36 1965* 1966 CALL ZSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, M, 1967 $ RWORK, RESULT( 35 ) ) 1968* 1969* Call ZSTEMR to compute D2, do tests. 1970* 1971* Compute D2 1972* 1973 CALL DCOPY( N, SD, 1, D5, 1 ) 1974 IF( N.GT.0 ) 1975 $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) 1976* 1977 NTEST = 37 1978 CALL ZSTEMR( 'N', 'A', N, D5, RWORK, VL, VU, IL, IU, 1979 $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC, 1980 $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), 1981 $ LIWORK-2*N, IINFO ) 1982 IF( IINFO.NE.0 ) THEN 1983 WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(N,A)', IINFO, N, 1984 $ JTYPE, IOLDSD 1985 INFO = ABS( IINFO ) 1986 IF( IINFO.LT.0 ) THEN 1987 RETURN 1988 ELSE 1989 RESULT( 37 ) = ULPINV 1990 GO TO 280 1991 END IF 1992 END IF 1993* 1994* Do Test 34 1995* 1996 TEMP1 = ZERO 1997 TEMP2 = ZERO 1998* 1999 DO 260 J = 1, N 2000 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) ) 2001 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 2002 260 CONTINUE 2003* 2004 RESULT( 37 ) = TEMP2 / MAX( UNFL, 2005 $ ULP*MAX( TEMP1, TEMP2 ) ) 2006 END IF 2007 270 CONTINUE 2008 280 CONTINUE 2009 NTESTT = NTESTT + NTEST 2010* 2011* End of Loop -- Check for RESULT(j) > THRESH 2012* 2013* 2014* Print out tests which fail. 2015* 2016 DO 290 JR = 1, NTEST 2017 IF( RESULT( JR ).GE.THRESH ) THEN 2018* 2019* If this is the first test to fail, 2020* print a header to the data file. 2021* 2022 IF( NERRS.EQ.0 ) THEN 2023 WRITE( NOUNIT, FMT = 9998 )'ZST' 2024 WRITE( NOUNIT, FMT = 9997 ) 2025 WRITE( NOUNIT, FMT = 9996 ) 2026 WRITE( NOUNIT, FMT = 9995 )'Hermitian' 2027 WRITE( NOUNIT, FMT = 9994 ) 2028* 2029* Tests performed 2030* 2031 WRITE( NOUNIT, FMT = 9987 ) 2032 END IF 2033 NERRS = NERRS + 1 2034 IF( RESULT( JR ).LT.10000.0D0 ) THEN 2035 WRITE( NOUNIT, FMT = 9989 )N, JTYPE, IOLDSD, JR, 2036 $ RESULT( JR ) 2037 ELSE 2038 WRITE( NOUNIT, FMT = 9988 )N, JTYPE, IOLDSD, JR, 2039 $ RESULT( JR ) 2040 END IF 2041 END IF 2042 290 CONTINUE 2043 300 CONTINUE 2044 310 CONTINUE 2045* 2046* Summary 2047* 2048 CALL DLASUM( 'ZST', NOUNIT, NERRS, NTESTT ) 2049 RETURN 2050* 2051 9999 FORMAT( ' ZCHKST2STG: ', A, ' returned INFO=', I6, '.', / 9X, 2052 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 2053* 2054 9998 FORMAT( / 1X, A3, ' -- Complex Hermitian eigenvalue problem' ) 2055 9997 FORMAT( ' Matrix types (see ZCHKST2STG for details): ' ) 2056* 2057 9996 FORMAT( / ' Special Matrices:', 2058 $ / ' 1=Zero matrix. ', 2059 $ ' 5=Diagonal: clustered entries.', 2060 $ / ' 2=Identity matrix. ', 2061 $ ' 6=Diagonal: large, evenly spaced.', 2062 $ / ' 3=Diagonal: evenly spaced entries. ', 2063 $ ' 7=Diagonal: small, evenly spaced.', 2064 $ / ' 4=Diagonal: geometr. spaced entries.' ) 2065 9995 FORMAT( ' Dense ', A, ' Matrices:', 2066 $ / ' 8=Evenly spaced eigenvals. ', 2067 $ ' 12=Small, evenly spaced eigenvals.', 2068 $ / ' 9=Geometrically spaced eigenvals. ', 2069 $ ' 13=Matrix with random O(1) entries.', 2070 $ / ' 10=Clustered eigenvalues. ', 2071 $ ' 14=Matrix with large random entries.', 2072 $ / ' 11=Large, evenly spaced eigenvals. ', 2073 $ ' 15=Matrix with small random entries.' ) 2074 9994 FORMAT( ' 16=Positive definite, evenly spaced eigenvalues', 2075 $ / ' 17=Positive definite, geometrically spaced eigenvlaues', 2076 $ / ' 18=Positive definite, clustered eigenvalues', 2077 $ / ' 19=Positive definite, small evenly spaced eigenvalues', 2078 $ / ' 20=Positive definite, large evenly spaced eigenvalues', 2079 $ / ' 21=Diagonally dominant tridiagonal, geometrically', 2080 $ ' spaced eigenvalues' ) 2081* 2082 9989 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 2083 $ 4( I4, ',' ), ' result ', I3, ' is', 0P, F8.2 ) 2084 9988 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', 2085 $ 4( I4, ',' ), ' result ', I3, ' is', 1P, D10.3 ) 2086* 2087 9987 FORMAT( / 'Test performed: see ZCHKST2STG for details.', / ) 2088* 2089* End of ZCHKST2STG 2090* 2091 END 2092