1*> \brief \b SCHKSBSTG 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 SCHKSB2STG( NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE, 12* ISEED, THRESH, NOUNIT, A, LDA, SD, SE, D1, 13* D2, D3, U, LDU, WORK, LWORK, RESULT, INFO ) 14* 15* .. Scalar Arguments .. 16* INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES, 17* $ NWDTHS 18* REAL THRESH 19* .. 20* .. Array Arguments .. 21* LOGICAL DOTYPE( * ) 22* INTEGER ISEED( 4 ), KK( * ), NN( * ) 23* REAL A( LDA, * ), RESULT( * ), SD( * ), SE( * ), 24* $ U( LDU, * ), WORK( * ) 25* .. 26* 27* 28*> \par Purpose: 29* ============= 30*> 31*> \verbatim 32*> 33*> SCHKSBSTG tests the reduction of a symmetric band matrix to tridiagonal 34*> form, used with the symmetric eigenvalue problem. 35*> 36*> SSBTRD factors a symmetric band matrix A as U S U' , where ' means 37*> transpose, S is symmetric tridiagonal, and U is orthogonal. 38*> SSBTRD can use either just the lower or just the upper triangle 39*> of A; SCHKSBSTG checks both cases. 40*> 41*> SSYTRD_SB2ST factors a symmetric band matrix A as U S U' , 42*> where ' means transpose, S is symmetric tridiagonal, and U is 43*> orthogonal. SSYTRD_SB2ST can use either just the lower or just 44*> the upper triangle of A; SCHKSBSTG checks both cases. 45*> 46*> SSTEQR factors S as Z D1 Z'. 47*> D1 is the matrix of eigenvalues computed when Z is not computed 48*> and from the S resulting of SSBTRD "U" (used as reference for SSYTRD_SB2ST) 49*> D2 is the matrix of eigenvalues computed when Z is not computed 50*> and from the S resulting of SSYTRD_SB2ST "U". 51*> D3 is the matrix of eigenvalues computed when Z is not computed 52*> and from the S resulting of SSYTRD_SB2ST "L". 53*> 54*> When SCHKSBSTG is called, a number of matrix "sizes" ("n's"), a number 55*> of bandwidths ("k's"), and a number of matrix "types" are 56*> specified. For each size ("n"), each bandwidth ("k") less than or 57*> equal to "n", and each type of matrix, one matrix will be generated 58*> and used to test the symmetric banded reduction routine. For each 59*> matrix, a number of tests will be performed: 60*> 61*> (1) | A - V S V' | / ( |A| n ulp ) computed by SSBTRD with 62*> UPLO='U' 63*> 64*> (2) | I - UU' | / ( n ulp ) 65*> 66*> (3) | A - V S V' | / ( |A| n ulp ) computed by SSBTRD with 67*> UPLO='L' 68*> 69*> (4) | I - UU' | / ( n ulp ) 70*> 71*> (5) | D1 - D2 | / ( |D1| ulp ) where D1 is computed by 72*> SSBTRD with UPLO='U' and 73*> D2 is computed by 74*> SSYTRD_SB2ST with UPLO='U' 75*> 76*> (6) | D1 - D3 | / ( |D1| ulp ) where D1 is computed by 77*> SSBTRD with UPLO='U' and 78*> D3 is computed by 79*> SSYTRD_SB2ST with UPLO='L' 80*> 81*> The "sizes" are specified by an array NN(1:NSIZES); the value of 82*> each element NN(j) specifies one size. 83*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); 84*> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 85*> Currently, the list of possible types is: 86*> 87*> (1) The zero matrix. 88*> (2) The identity matrix. 89*> 90*> (3) A diagonal matrix with evenly spaced entries 91*> 1, ..., ULP and random signs. 92*> (ULP = (first number larger than 1) - 1 ) 93*> (4) A diagonal matrix with geometrically spaced entries 94*> 1, ..., ULP and random signs. 95*> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 96*> and random signs. 97*> 98*> (6) Same as (4), but multiplied by SQRT( overflow threshold ) 99*> (7) Same as (4), but multiplied by SQRT( underflow threshold ) 100*> 101*> (8) A matrix of the form U' D U, where U is orthogonal and 102*> D has evenly spaced entries 1, ..., ULP with random signs 103*> on the diagonal. 104*> 105*> (9) A matrix of the form U' D U, where U is orthogonal and 106*> D has geometrically spaced entries 1, ..., ULP with random 107*> signs on the diagonal. 108*> 109*> (10) A matrix of the form U' D U, where U is orthogonal and 110*> D has "clustered" entries 1, ULP,..., ULP with random 111*> signs on the diagonal. 112*> 113*> (11) Same as (8), but multiplied by SQRT( overflow threshold ) 114*> (12) Same as (8), but multiplied by SQRT( underflow threshold ) 115*> 116*> (13) Symmetric matrix with random entries chosen from (-1,1). 117*> (14) Same as (13), but multiplied by SQRT( overflow threshold ) 118*> (15) Same as (13), but multiplied by SQRT( underflow threshold ) 119*> \endverbatim 120* 121* Arguments: 122* ========== 123* 124*> \param[in] NSIZES 125*> \verbatim 126*> NSIZES is INTEGER 127*> The number of sizes of matrices to use. If it is zero, 128*> SCHKSBSTG does nothing. It must be at least zero. 129*> \endverbatim 130*> 131*> \param[in] NN 132*> \verbatim 133*> NN is INTEGER array, dimension (NSIZES) 134*> An array containing the sizes to be used for the matrices. 135*> Zero values will be skipped. The values must be at least 136*> zero. 137*> \endverbatim 138*> 139*> \param[in] NWDTHS 140*> \verbatim 141*> NWDTHS is INTEGER 142*> The number of bandwidths to use. If it is zero, 143*> SCHKSBSTG does nothing. It must be at least zero. 144*> \endverbatim 145*> 146*> \param[in] KK 147*> \verbatim 148*> KK is INTEGER array, dimension (NWDTHS) 149*> An array containing the bandwidths to be used for the band 150*> matrices. The values must be at least zero. 151*> \endverbatim 152*> 153*> \param[in] NTYPES 154*> \verbatim 155*> NTYPES is INTEGER 156*> The number of elements in DOTYPE. If it is zero, SCHKSBSTG 157*> does nothing. It must be at least zero. If it is MAXTYP+1 158*> and NSIZES is 1, then an additional type, MAXTYP+1 is 159*> defined, which is to use whatever matrix is in A. This 160*> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 161*> DOTYPE(MAXTYP+1) is .TRUE. . 162*> \endverbatim 163*> 164*> \param[in] DOTYPE 165*> \verbatim 166*> DOTYPE is LOGICAL array, dimension (NTYPES) 167*> If DOTYPE(j) is .TRUE., then for each size in NN a 168*> matrix of that size and of type j will be generated. 169*> If NTYPES is smaller than the maximum number of types 170*> defined (PARAMETER MAXTYP), then types NTYPES+1 through 171*> MAXTYP will not be generated. If NTYPES is larger 172*> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 173*> will be ignored. 174*> \endverbatim 175*> 176*> \param[in,out] ISEED 177*> \verbatim 178*> ISEED is INTEGER array, dimension (4) 179*> On entry ISEED specifies the seed of the random number 180*> generator. The array elements should be between 0 and 4095; 181*> if not they will be reduced mod 4096. Also, ISEED(4) must 182*> be odd. The random number generator uses a linear 183*> congruential sequence limited to small integers, and so 184*> should produce machine independent random numbers. The 185*> values of ISEED are changed on exit, and can be used in the 186*> next call to SCHKSBSTG to continue the same random number 187*> sequence. 188*> \endverbatim 189*> 190*> \param[in] THRESH 191*> \verbatim 192*> THRESH is REAL 193*> A test will count as "failed" if the "error", computed as 194*> described above, exceeds THRESH. Note that the error 195*> is scaled to be O(1), so THRESH should be a reasonably 196*> small multiple of 1, e.g., 10 or 100. In particular, 197*> it should not depend on the precision (single vs. double) 198*> or the size of the matrix. It must be at least zero. 199*> \endverbatim 200*> 201*> \param[in] NOUNIT 202*> \verbatim 203*> NOUNIT is INTEGER 204*> The FORTRAN unit number for printing out error messages 205*> (e.g., if a routine returns IINFO not equal to 0.) 206*> \endverbatim 207*> 208*> \param[in,out] A 209*> \verbatim 210*> A is REAL array, dimension 211*> (LDA, max(NN)) 212*> Used to hold the matrix whose eigenvalues are to be 213*> computed. 214*> \endverbatim 215*> 216*> \param[in] LDA 217*> \verbatim 218*> LDA is INTEGER 219*> The leading dimension of A. It must be at least 2 (not 1!) 220*> and at least max( KK )+1. 221*> \endverbatim 222*> 223*> \param[out] SD 224*> \verbatim 225*> SD is REAL array, dimension (max(NN)) 226*> Used to hold the diagonal of the tridiagonal matrix computed 227*> by SSBTRD. 228*> \endverbatim 229*> 230*> \param[out] SE 231*> \verbatim 232*> SE is REAL array, dimension (max(NN)) 233*> Used to hold the off-diagonal of the tridiagonal matrix 234*> computed by SSBTRD. 235*> \endverbatim 236*> 237*> \param[out] U 238*> \verbatim 239*> U is REAL array, dimension (LDU, max(NN)) 240*> Used to hold the orthogonal matrix computed by SSBTRD. 241*> \endverbatim 242*> 243*> \param[in] LDU 244*> \verbatim 245*> LDU is INTEGER 246*> The leading dimension of U. It must be at least 1 247*> and at least max( NN ). 248*> \endverbatim 249*> 250*> \param[out] WORK 251*> \verbatim 252*> WORK is REAL array, dimension (LWORK) 253*> \endverbatim 254*> 255*> \param[in] LWORK 256*> \verbatim 257*> LWORK is INTEGER 258*> The number of entries in WORK. This must be at least 259*> max( LDA+1, max(NN)+1 )*max(NN). 260*> \endverbatim 261*> 262*> \param[out] RESULT 263*> \verbatim 264*> RESULT is REAL array, dimension (4) 265*> The values computed by the tests described above. 266*> The values are currently limited to 1/ulp, to avoid 267*> overflow. 268*> \endverbatim 269*> 270*> \param[out] INFO 271*> \verbatim 272*> INFO is INTEGER 273*> If 0, then everything ran OK. 274*> 275*>----------------------------------------------------------------------- 276*> 277*> Some Local Variables and Parameters: 278*> ---- ----- --------- --- ---------- 279*> ZERO, ONE Real 0 and 1. 280*> MAXTYP The number of types defined. 281*> NTEST The number of tests performed, or which can 282*> be performed so far, for the current matrix. 283*> NTESTT The total number of tests performed so far. 284*> NMAX Largest value in NN. 285*> NMATS The number of matrices generated so far. 286*> NERRS The number of tests which have exceeded THRESH 287*> so far. 288*> COND, IMODE Values to be passed to the matrix generators. 289*> ANORM Norm of A; passed to matrix generators. 290*> 291*> OVFL, UNFL Overflow and underflow thresholds. 292*> ULP, ULPINV Finest relative precision and its inverse. 293*> RTOVFL, RTUNFL Square roots of the previous 2 values. 294*> The following four arrays decode JTYPE: 295*> KTYPE(j) The general type (1-10) for type "j". 296*> KMODE(j) The MODE value to be passed to the matrix 297*> generator for type "j". 298*> KMAGN(j) The order of magnitude ( O(1), 299*> O(overflow^(1/2) ), O(underflow^(1/2) ) 300*> \endverbatim 301* 302* Authors: 303* ======== 304* 305*> \author Univ. of Tennessee 306*> \author Univ. of California Berkeley 307*> \author Univ. of Colorado Denver 308*> \author NAG Ltd. 309* 310*> \date June 2017 311* 312*> \ingroup single_eig 313* 314* ===================================================================== 315 SUBROUTINE SCHKSB2STG( NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE, 316 $ ISEED, THRESH, NOUNIT, A, LDA, SD, SE, D1, 317 $ D2, D3, U, LDU, WORK, LWORK, RESULT, INFO ) 318* 319* -- LAPACK test routine (version 3.7.1) -- 320* -- LAPACK is a software package provided by Univ. of Tennessee, -- 321* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 322* June 2017 323* 324* .. Scalar Arguments .. 325 INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES, 326 $ NWDTHS 327 REAL THRESH 328* .. 329* .. Array Arguments .. 330 LOGICAL DOTYPE( * ) 331 INTEGER ISEED( 4 ), KK( * ), NN( * ) 332 REAL A( LDA, * ), RESULT( * ), SD( * ), SE( * ), 333 $ D1( * ), D2( * ), D3( * ), 334 $ U( LDU, * ), WORK( * ) 335* .. 336* 337* ===================================================================== 338* 339* .. Parameters .. 340 REAL ZERO, ONE, TWO, TEN 341 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 342 $ TEN = 10.0E0 ) 343 REAL HALF 344 PARAMETER ( HALF = ONE / TWO ) 345 INTEGER MAXTYP 346 PARAMETER ( MAXTYP = 15 ) 347* .. 348* .. Local Scalars .. 349 LOGICAL BADNN, BADNNB 350 INTEGER I, IINFO, IMODE, ITYPE, J, JC, JCOL, JR, JSIZE, 351 $ JTYPE, JWIDTH, K, KMAX, LH, LW, MTYPES, N, 352 $ NERRS, NMATS, NMAX, NTEST, NTESTT 353 REAL ANINV, ANORM, COND, OVFL, RTOVFL, RTUNFL, 354 $ TEMP1, TEMP2, TEMP3, TEMP4, ULP, ULPINV, UNFL 355* .. 356* .. Local Arrays .. 357 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KMAGN( MAXTYP ), 358 $ KMODE( MAXTYP ), KTYPE( MAXTYP ) 359* .. 360* .. External Functions .. 361 REAL SLAMCH 362 EXTERNAL SLAMCH 363* .. 364* .. External Subroutines .. 365 EXTERNAL SLACPY, SLASET, SLASUM, SLATMR, SLATMS, SSBT21, 366 $ SSBTRD, XERBLA, SSYTRD_SB2ST, SSTEQR 367* .. 368* .. Intrinsic Functions .. 369 INTRINSIC ABS, REAL, MAX, MIN, SQRT 370* .. 371* .. Data statements .. 372 DATA KTYPE / 1, 2, 5*4, 5*5, 3*8 / 373 DATA KMAGN / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1, 374 $ 2, 3 / 375 DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0, 376 $ 0, 0 / 377* .. 378* .. Executable Statements .. 379* 380* Check for errors 381* 382 NTESTT = 0 383 INFO = 0 384* 385* Important constants 386* 387 BADNN = .FALSE. 388 NMAX = 1 389 DO 10 J = 1, NSIZES 390 NMAX = MAX( NMAX, NN( J ) ) 391 IF( NN( J ).LT.0 ) 392 $ BADNN = .TRUE. 393 10 CONTINUE 394* 395 BADNNB = .FALSE. 396 KMAX = 0 397 DO 20 J = 1, NSIZES 398 KMAX = MAX( KMAX, KK( J ) ) 399 IF( KK( J ).LT.0 ) 400 $ BADNNB = .TRUE. 401 20 CONTINUE 402 KMAX = MIN( NMAX-1, KMAX ) 403* 404* Check for errors 405* 406 IF( NSIZES.LT.0 ) THEN 407 INFO = -1 408 ELSE IF( BADNN ) THEN 409 INFO = -2 410 ELSE IF( NWDTHS.LT.0 ) THEN 411 INFO = -3 412 ELSE IF( BADNNB ) THEN 413 INFO = -4 414 ELSE IF( NTYPES.LT.0 ) THEN 415 INFO = -5 416 ELSE IF( LDA.LT.KMAX+1 ) THEN 417 INFO = -11 418 ELSE IF( LDU.LT.NMAX ) THEN 419 INFO = -15 420 ELSE IF( ( MAX( LDA, NMAX )+1 )*NMAX.GT.LWORK ) THEN 421 INFO = -17 422 END IF 423* 424 IF( INFO.NE.0 ) THEN 425 CALL XERBLA( 'SCHKSBSTG', -INFO ) 426 RETURN 427 END IF 428* 429* Quick return if possible 430* 431 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 .OR. NWDTHS.EQ.0 ) 432 $ RETURN 433* 434* More Important constants 435* 436 UNFL = SLAMCH( 'Safe minimum' ) 437 OVFL = ONE / UNFL 438 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 439 ULPINV = ONE / ULP 440 RTUNFL = SQRT( UNFL ) 441 RTOVFL = SQRT( OVFL ) 442* 443* Loop over sizes, types 444* 445 NERRS = 0 446 NMATS = 0 447* 448 DO 190 JSIZE = 1, NSIZES 449 N = NN( JSIZE ) 450 ANINV = ONE / REAL( MAX( 1, N ) ) 451* 452 DO 180 JWIDTH = 1, NWDTHS 453 K = KK( JWIDTH ) 454 IF( K.GT.N ) 455 $ GO TO 180 456 K = MAX( 0, MIN( N-1, K ) ) 457* 458 IF( NSIZES.NE.1 ) THEN 459 MTYPES = MIN( MAXTYP, NTYPES ) 460 ELSE 461 MTYPES = MIN( MAXTYP+1, NTYPES ) 462 END IF 463* 464 DO 170 JTYPE = 1, MTYPES 465 IF( .NOT.DOTYPE( JTYPE ) ) 466 $ GO TO 170 467 NMATS = NMATS + 1 468 NTEST = 0 469* 470 DO 30 J = 1, 4 471 IOLDSD( J ) = ISEED( J ) 472 30 CONTINUE 473* 474* Compute "A". 475* Store as "Upper"; later, we will copy to other format. 476* 477* Control parameters: 478* 479* KMAGN KMODE KTYPE 480* =1 O(1) clustered 1 zero 481* =2 large clustered 2 identity 482* =3 small exponential (none) 483* =4 arithmetic diagonal, (w/ eigenvalues) 484* =5 random log symmetric, w/ eigenvalues 485* =6 random (none) 486* =7 random diagonal 487* =8 random symmetric 488* =9 positive definite 489* =10 diagonally dominant tridiagonal 490* 491 IF( MTYPES.GT.MAXTYP ) 492 $ GO TO 100 493* 494 ITYPE = KTYPE( JTYPE ) 495 IMODE = KMODE( JTYPE ) 496* 497* Compute norm 498* 499 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 500* 501 40 CONTINUE 502 ANORM = ONE 503 GO TO 70 504* 505 50 CONTINUE 506 ANORM = ( RTOVFL*ULP )*ANINV 507 GO TO 70 508* 509 60 CONTINUE 510 ANORM = RTUNFL*N*ULPINV 511 GO TO 70 512* 513 70 CONTINUE 514* 515 CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 516 IINFO = 0 517 IF( JTYPE.LE.15 ) THEN 518 COND = ULPINV 519 ELSE 520 COND = ULPINV*ANINV / TEN 521 END IF 522* 523* Special Matrices -- Identity & Jordan block 524* 525* Zero 526* 527 IF( ITYPE.EQ.1 ) THEN 528 IINFO = 0 529* 530 ELSE IF( ITYPE.EQ.2 ) THEN 531* 532* Identity 533* 534 DO 80 JCOL = 1, N 535 A( K+1, JCOL ) = ANORM 536 80 CONTINUE 537* 538 ELSE IF( ITYPE.EQ.4 ) THEN 539* 540* Diagonal Matrix, [Eigen]values Specified 541* 542 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 543 $ ANORM, 0, 0, 'Q', A( K+1, 1 ), LDA, 544 $ WORK( N+1 ), IINFO ) 545* 546 ELSE IF( ITYPE.EQ.5 ) THEN 547* 548* Symmetric, eigenvalues specified 549* 550 CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND, 551 $ ANORM, K, K, 'Q', A, LDA, WORK( N+1 ), 552 $ IINFO ) 553* 554 ELSE IF( ITYPE.EQ.7 ) THEN 555* 556* Diagonal, random eigenvalues 557* 558 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 559 $ 'T', 'N', WORK( N+1 ), 1, ONE, 560 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, 561 $ ZERO, ANORM, 'Q', A( K+1, 1 ), LDA, 562 $ IDUMMA, IINFO ) 563* 564 ELSE IF( ITYPE.EQ.8 ) THEN 565* 566* Symmetric, random eigenvalues 567* 568 CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE, 569 $ 'T', 'N', WORK( N+1 ), 1, ONE, 570 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, K, K, 571 $ ZERO, ANORM, 'Q', A, LDA, IDUMMA, IINFO ) 572* 573 ELSE IF( ITYPE.EQ.9 ) THEN 574* 575* Positive definite, eigenvalues specified. 576* 577 CALL SLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND, 578 $ ANORM, K, K, 'Q', A, LDA, WORK( N+1 ), 579 $ IINFO ) 580* 581 ELSE IF( ITYPE.EQ.10 ) THEN 582* 583* Positive definite tridiagonal, eigenvalues specified. 584* 585 IF( N.GT.1 ) 586 $ K = MAX( 1, K ) 587 CALL SLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND, 588 $ ANORM, 1, 1, 'Q', A( K, 1 ), LDA, 589 $ WORK( N+1 ), IINFO ) 590 DO 90 I = 2, N 591 TEMP1 = ABS( A( K, I ) ) / 592 $ SQRT( ABS( A( K+1, I-1 )*A( K+1, I ) ) ) 593 IF( TEMP1.GT.HALF ) THEN 594 A( K, I ) = HALF*SQRT( ABS( A( K+1, 595 $ I-1 )*A( K+1, I ) ) ) 596 END IF 597 90 CONTINUE 598* 599 ELSE 600* 601 IINFO = 1 602 END IF 603* 604 IF( IINFO.NE.0 ) THEN 605 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, 606 $ JTYPE, IOLDSD 607 INFO = ABS( IINFO ) 608 RETURN 609 END IF 610* 611 100 CONTINUE 612* 613* Call SSBTRD to compute S and U from upper triangle. 614* 615 CALL SLACPY( ' ', K+1, N, A, LDA, WORK, LDA ) 616* 617 NTEST = 1 618 CALL SSBTRD( 'V', 'U', N, K, WORK, LDA, SD, SE, U, LDU, 619 $ WORK( LDA*N+1 ), IINFO ) 620* 621 IF( IINFO.NE.0 ) THEN 622 WRITE( NOUNIT, FMT = 9999 )'SSBTRD(U)', IINFO, N, 623 $ JTYPE, IOLDSD 624 INFO = ABS( IINFO ) 625 IF( IINFO.LT.0 ) THEN 626 RETURN 627 ELSE 628 RESULT( 1 ) = ULPINV 629 GO TO 150 630 END IF 631 END IF 632* 633* Do tests 1 and 2 634* 635 CALL SSBT21( 'Upper', N, K, 1, A, LDA, SD, SE, U, LDU, 636 $ WORK, RESULT( 1 ) ) 637* 638* Before converting A into lower for SSBTRD, run SSYTRD_SB2ST 639* otherwise matrix A will be converted to lower and then need 640* to be converted back to upper in order to run the upper case 641* ofSSYTRD_SB2ST 642* 643* Compute D1 the eigenvalues resulting from the tridiagonal 644* form using the SSBTRD and used as reference to compare 645* with the SSYTRD_SB2ST routine 646* 647* Compute D1 from the SSBTRD and used as reference for the 648* SSYTRD_SB2ST 649* 650 CALL SCOPY( N, SD, 1, D1, 1 ) 651 IF( N.GT.0 ) 652 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 653* 654 CALL SSTEQR( 'N', N, D1, WORK, WORK( N+1 ), LDU, 655 $ WORK( N+1 ), IINFO ) 656 IF( IINFO.NE.0 ) THEN 657 WRITE( NOUNIT, FMT = 9999 )'SSTEQR(N)', IINFO, N, 658 $ JTYPE, IOLDSD 659 INFO = ABS( IINFO ) 660 IF( IINFO.LT.0 ) THEN 661 RETURN 662 ELSE 663 RESULT( 5 ) = ULPINV 664 GO TO 150 665 END IF 666 END IF 667* 668* SSYTRD_SB2ST Upper case is used to compute D2. 669* Note to set SD and SE to zero to be sure not reusing 670* the one from above. Compare it with D1 computed 671* using the SSBTRD. 672* 673 CALL SLASET( 'Full', N, 1, ZERO, ZERO, SD, N ) 674 CALL SLASET( 'Full', N, 1, ZERO, ZERO, SE, N ) 675 CALL SLACPY( ' ', K+1, N, A, LDA, U, LDU ) 676 LH = MAX(1, 4*N) 677 LW = LWORK - LH 678 CALL SSYTRD_SB2ST( 'N', 'N', "U", N, K, U, LDU, SD, SE, 679 $ WORK, LH, WORK( LH+1 ), LW, IINFO ) 680* 681* Compute D2 from the SSYTRD_SB2ST Upper case 682* 683 CALL SCOPY( N, SD, 1, D2, 1 ) 684 IF( N.GT.0 ) 685 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 686* 687 CALL SSTEQR( 'N', N, D2, WORK, WORK( N+1 ), LDU, 688 $ WORK( N+1 ), IINFO ) 689 IF( IINFO.NE.0 ) THEN 690 WRITE( NOUNIT, FMT = 9999 )'SSTEQR(N)', IINFO, N, 691 $ JTYPE, IOLDSD 692 INFO = ABS( IINFO ) 693 IF( IINFO.LT.0 ) THEN 694 RETURN 695 ELSE 696 RESULT( 5 ) = ULPINV 697 GO TO 150 698 END IF 699 END IF 700* 701* Convert A from Upper-Triangle-Only storage to 702* Lower-Triangle-Only storage. 703* 704 DO 120 JC = 1, N 705 DO 110 JR = 0, MIN( K, N-JC ) 706 A( JR+1, JC ) = A( K+1-JR, JC+JR ) 707 110 CONTINUE 708 120 CONTINUE 709 DO 140 JC = N + 1 - K, N 710 DO 130 JR = MIN( K, N-JC ) + 1, K 711 A( JR+1, JC ) = ZERO 712 130 CONTINUE 713 140 CONTINUE 714* 715* Call SSBTRD to compute S and U from lower triangle 716* 717 CALL SLACPY( ' ', K+1, N, A, LDA, WORK, LDA ) 718* 719 NTEST = 3 720 CALL SSBTRD( 'V', 'L', N, K, WORK, LDA, SD, SE, U, LDU, 721 $ WORK( LDA*N+1 ), IINFO ) 722* 723 IF( IINFO.NE.0 ) THEN 724 WRITE( NOUNIT, FMT = 9999 )'SSBTRD(L)', IINFO, N, 725 $ JTYPE, IOLDSD 726 INFO = ABS( IINFO ) 727 IF( IINFO.LT.0 ) THEN 728 RETURN 729 ELSE 730 RESULT( 3 ) = ULPINV 731 GO TO 150 732 END IF 733 END IF 734 NTEST = 4 735* 736* Do tests 3 and 4 737* 738 CALL SSBT21( 'Lower', N, K, 1, A, LDA, SD, SE, U, LDU, 739 $ WORK, RESULT( 3 ) ) 740* 741* SSYTRD_SB2ST Lower case is used to compute D3. 742* Note to set SD and SE to zero to be sure not reusing 743* the one from above. Compare it with D1 computed 744* using the SSBTRD. 745* 746 CALL SLASET( 'Full', N, 1, ZERO, ZERO, SD, 1 ) 747 CALL SLASET( 'Full', N, 1, ZERO, ZERO, SE, 1 ) 748 CALL SLACPY( ' ', K+1, N, A, LDA, U, LDU ) 749 LH = MAX(1, 4*N) 750 LW = LWORK - LH 751 CALL SSYTRD_SB2ST( 'N', 'N', "L", N, K, U, LDU, SD, SE, 752 $ WORK, LH, WORK( LH+1 ), LW, IINFO ) 753* 754* Compute D3 from the 2-stage Upper case 755* 756 CALL SCOPY( N, SD, 1, D3, 1 ) 757 IF( N.GT.0 ) 758 $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) 759* 760 CALL SSTEQR( 'N', N, D3, WORK, WORK( N+1 ), LDU, 761 $ WORK( N+1 ), IINFO ) 762 IF( IINFO.NE.0 ) THEN 763 WRITE( NOUNIT, FMT = 9999 )'SSTEQR(N)', IINFO, N, 764 $ JTYPE, IOLDSD 765 INFO = ABS( IINFO ) 766 IF( IINFO.LT.0 ) THEN 767 RETURN 768 ELSE 769 RESULT( 6 ) = ULPINV 770 GO TO 150 771 END IF 772 END IF 773* 774* 775* Do Tests 3 and 4 which are similar to 11 and 12 but with the 776* D1 computed using the standard 1-stage reduction as reference 777* 778 NTEST = 6 779 TEMP1 = ZERO 780 TEMP2 = ZERO 781 TEMP3 = ZERO 782 TEMP4 = ZERO 783* 784 DO 151 J = 1, N 785 TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) ) 786 TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) ) 787 TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) ) 788 TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) ) 789 151 CONTINUE 790* 791 RESULT(5) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) 792 RESULT(6) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) ) 793* 794* End of Loop -- Check for RESULT(j) > THRESH 795* 796 150 CONTINUE 797 NTESTT = NTESTT + NTEST 798* 799* Print out tests which fail. 800* 801 DO 160 JR = 1, NTEST 802 IF( RESULT( JR ).GE.THRESH ) THEN 803* 804* If this is the first test to fail, 805* print a header to the data file. 806* 807 IF( NERRS.EQ.0 ) THEN 808 WRITE( NOUNIT, FMT = 9998 )'SSB' 809 WRITE( NOUNIT, FMT = 9997 ) 810 WRITE( NOUNIT, FMT = 9996 ) 811 WRITE( NOUNIT, FMT = 9995 )'Symmetric' 812 WRITE( NOUNIT, FMT = 9994 )'orthogonal', '''', 813 $ 'transpose', ( '''', J = 1, 6 ) 814 END IF 815 NERRS = NERRS + 1 816 WRITE( NOUNIT, FMT = 9993 )N, K, IOLDSD, JTYPE, 817 $ JR, RESULT( JR ) 818 END IF 819 160 CONTINUE 820* 821 170 CONTINUE 822 180 CONTINUE 823 190 CONTINUE 824* 825* Summary 826* 827 CALL SLASUM( 'SSB', NOUNIT, NERRS, NTESTT ) 828 RETURN 829* 830 9999 FORMAT( ' SCHKSBSTG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', 831 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) 832* 833 9998 FORMAT( / 1X, A3, 834 $ ' -- Real Symmetric Banded Tridiagonal Reduction Routines' ) 835 9997 FORMAT( ' Matrix types (see SCHKSBSTG for details): ' ) 836* 837 9996 FORMAT( / ' Special Matrices:', 838 $ / ' 1=Zero matrix. ', 839 $ ' 5=Diagonal: clustered entries.', 840 $ / ' 2=Identity matrix. ', 841 $ ' 6=Diagonal: large, evenly spaced.', 842 $ / ' 3=Diagonal: evenly spaced entries. ', 843 $ ' 7=Diagonal: small, evenly spaced.', 844 $ / ' 4=Diagonal: geometr. spaced entries.' ) 845 9995 FORMAT( ' Dense ', A, ' Banded Matrices:', 846 $ / ' 8=Evenly spaced eigenvals. ', 847 $ ' 12=Small, evenly spaced eigenvals.', 848 $ / ' 9=Geometrically spaced eigenvals. ', 849 $ ' 13=Matrix with random O(1) entries.', 850 $ / ' 10=Clustered eigenvalues. ', 851 $ ' 14=Matrix with large random entries.', 852 $ / ' 11=Large, evenly spaced eigenvals. ', 853 $ ' 15=Matrix with small random entries.' ) 854* 855 9994 FORMAT( / ' Tests performed: (S is Tridiag, U is ', A, ',', 856 $ / 20X, A, ' means ', A, '.', / ' UPLO=''U'':', 857 $ / ' 1= | A - U S U', A1, ' | / ( |A| n ulp ) ', 858 $ ' 2= | I - U U', A1, ' | / ( n ulp )', / ' UPLO=''L'':', 859 $ / ' 3= | A - U S U', A1, ' | / ( |A| n ulp ) ', 860 $ ' 4= | I - U U', A1, ' | / ( n ulp )' / ' Eig check:', 861 $ /' 5= | D1 - D2', '', ' | / ( |D1| ulp ) ', 862 $ ' 6= | D1 - D3', '', ' | / ( |D1| ulp ) ' ) 863 9993 FORMAT( ' N=', I5, ', K=', I4, ', seed=', 4( I4, ',' ), ' type ', 864 $ I2, ', test(', I2, ')=', G10.3 ) 865* 866* End of SCHKSBSTG 867* 868 END 869