1*> \brief \b DCHKBB 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 DCHKBB( NSIZES, MVAL, NVAL, NWDTHS, KK, NTYPES, DOTYPE, 12* NRHS, ISEED, THRESH, NOUNIT, A, LDA, AB, LDAB, 13* BD, BE, Q, LDQ, P, LDP, C, LDC, CC, WORK, 14* LWORK, RESULT, INFO ) 15* 16* .. Scalar Arguments .. 17* INTEGER INFO, LDA, LDAB, LDC, LDP, LDQ, LWORK, NOUNIT, 18* $ NRHS, NSIZES, NTYPES, NWDTHS 19* DOUBLE PRECISION THRESH 20* .. 21* .. Array Arguments .. 22* LOGICAL DOTYPE( * ) 23* INTEGER ISEED( 4 ), KK( * ), MVAL( * ), NVAL( * ) 24* DOUBLE PRECISION A( LDA, * ), AB( LDAB, * ), BD( * ), BE( * ), 25* $ C( LDC, * ), CC( LDC, * ), P( LDP, * ), 26* $ Q( LDQ, * ), RESULT( * ), WORK( * ) 27* .. 28* 29* 30*> \par Purpose: 31* ============= 32*> 33*> \verbatim 34*> 35*> DCHKBB tests the reduction of a general real rectangular band 36*> matrix to bidiagonal form. 37*> 38*> DGBBRD factors a general band matrix A as Q B P* , where * means 39*> transpose, B is upper bidiagonal, and Q and P are orthogonal; 40*> DGBBRD can also overwrite a given matrix C with Q* C . 41*> 42*> For each pair of matrix dimensions (M,N) and each selected matrix 43*> type, an M by N matrix A and an M by NRHS matrix C are generated. 44*> The problem dimensions are as follows 45*> A: M x N 46*> Q: M x M 47*> P: N x N 48*> B: min(M,N) x min(M,N) 49*> C: M x NRHS 50*> 51*> For each generated matrix, 4 tests are performed: 52*> 53*> (1) | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P' 54*> 55*> (2) | I - Q' Q | / ( M ulp ) 56*> 57*> (3) | I - PT PT' | / ( N ulp ) 58*> 59*> (4) | Y - Q' C | / ( |Y| max(M,NRHS) ulp ), where Y = Q' C. 60*> 61*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); 62*> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. 63*> Currently, the list of possible types is: 64*> 65*> The possible matrix types are 66*> 67*> (1) The zero matrix. 68*> (2) The identity matrix. 69*> 70*> (3) A diagonal matrix with evenly spaced entries 71*> 1, ..., ULP and random signs. 72*> (ULP = (first number larger than 1) - 1 ) 73*> (4) A diagonal matrix with geometrically spaced entries 74*> 1, ..., ULP and random signs. 75*> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 76*> and random signs. 77*> 78*> (6) Same as (3), but multiplied by SQRT( overflow threshold ) 79*> (7) Same as (3), but multiplied by SQRT( underflow threshold ) 80*> 81*> (8) A matrix of the form U D V, where U and V are orthogonal and 82*> D has evenly spaced entries 1, ..., ULP with random signs 83*> on the diagonal. 84*> 85*> (9) A matrix of the form U D V, where U and V are orthogonal and 86*> D has geometrically spaced entries 1, ..., ULP with random 87*> signs on the diagonal. 88*> 89*> (10) A matrix of the form U D V, where U and V are orthogonal and 90*> D has "clustered" entries 1, ULP,..., ULP with random 91*> signs on the diagonal. 92*> 93*> (11) Same as (8), but multiplied by SQRT( overflow threshold ) 94*> (12) Same as (8), but multiplied by SQRT( underflow threshold ) 95*> 96*> (13) Rectangular matrix with random entries chosen from (-1,1). 97*> (14) Same as (13), but multiplied by SQRT( overflow threshold ) 98*> (15) Same as (13), but multiplied by SQRT( underflow threshold ) 99*> \endverbatim 100* 101* Arguments: 102* ========== 103* 104*> \param[in] NSIZES 105*> \verbatim 106*> NSIZES is INTEGER 107*> The number of values of M and N contained in the vectors 108*> MVAL and NVAL. The matrix sizes are used in pairs (M,N). 109*> If NSIZES is zero, DCHKBB does nothing. NSIZES must be at 110*> least zero. 111*> \endverbatim 112*> 113*> \param[in] MVAL 114*> \verbatim 115*> MVAL is INTEGER array, dimension (NSIZES) 116*> The values of the matrix row dimension M. 117*> \endverbatim 118*> 119*> \param[in] NVAL 120*> \verbatim 121*> NVAL is INTEGER array, dimension (NSIZES) 122*> The values of the matrix column dimension N. 123*> \endverbatim 124*> 125*> \param[in] NWDTHS 126*> \verbatim 127*> NWDTHS is INTEGER 128*> The number of bandwidths to use. If it is zero, 129*> DCHKBB does nothing. It must be at least zero. 130*> \endverbatim 131*> 132*> \param[in] KK 133*> \verbatim 134*> KK is INTEGER array, dimension (NWDTHS) 135*> An array containing the bandwidths to be used for the band 136*> matrices. The values must be at least zero. 137*> \endverbatim 138*> 139*> \param[in] NTYPES 140*> \verbatim 141*> NTYPES is INTEGER 142*> The number of elements in DOTYPE. If it is zero, DCHKBB 143*> does nothing. It must be at least zero. If it is MAXTYP+1 144*> and NSIZES is 1, then an additional type, MAXTYP+1 is 145*> defined, which is to use whatever matrix is in A. This 146*> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 147*> DOTYPE(MAXTYP+1) is .TRUE. . 148*> \endverbatim 149*> 150*> \param[in] DOTYPE 151*> \verbatim 152*> DOTYPE is LOGICAL array, dimension (NTYPES) 153*> If DOTYPE(j) is .TRUE., then for each size in NN a 154*> matrix of that size and of type j will be generated. 155*> If NTYPES is smaller than the maximum number of types 156*> defined (PARAMETER MAXTYP), then types NTYPES+1 through 157*> MAXTYP will not be generated. If NTYPES is larger 158*> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) 159*> will be ignored. 160*> \endverbatim 161*> 162*> \param[in] NRHS 163*> \verbatim 164*> NRHS is INTEGER 165*> The number of columns in the "right-hand side" matrix C. 166*> If NRHS = 0, then the operations on the right-hand side will 167*> not be tested. NRHS must be at least 0. 168*> \endverbatim 169*> 170*> \param[in,out] ISEED 171*> \verbatim 172*> ISEED is INTEGER array, dimension (4) 173*> On entry ISEED specifies the seed of the random number 174*> generator. The array elements should be between 0 and 4095; 175*> if not they will be reduced mod 4096. Also, ISEED(4) must 176*> be odd. The random number generator uses a linear 177*> congruential sequence limited to small integers, and so 178*> should produce machine independent random numbers. The 179*> values of ISEED are changed on exit, and can be used in the 180*> next call to DCHKBB to continue the same random number 181*> sequence. 182*> \endverbatim 183*> 184*> \param[in] THRESH 185*> \verbatim 186*> THRESH is DOUBLE PRECISION 187*> A test will count as "failed" if the "error", computed as 188*> described above, exceeds THRESH. Note that the error 189*> is scaled to be O(1), so THRESH should be a reasonably 190*> small multiple of 1, e.g., 10 or 100. In particular, 191*> it should not depend on the precision (single vs. double) 192*> or the size of the matrix. It must be at least zero. 193*> \endverbatim 194*> 195*> \param[in] NOUNIT 196*> \verbatim 197*> NOUNIT is INTEGER 198*> The FORTRAN unit number for printing out error messages 199*> (e.g., if a routine returns IINFO not equal to 0.) 200*> \endverbatim 201*> 202*> \param[in,out] A 203*> \verbatim 204*> A is DOUBLE PRECISION array, dimension 205*> (LDA, max(NN)) 206*> Used to hold the matrix A. 207*> \endverbatim 208*> 209*> \param[in] LDA 210*> \verbatim 211*> LDA is INTEGER 212*> The leading dimension of A. It must be at least 1 213*> and at least max( NN ). 214*> \endverbatim 215*> 216*> \param[out] AB 217*> \verbatim 218*> AB is DOUBLE PRECISION array, dimension (LDAB, max(NN)) 219*> Used to hold A in band storage format. 220*> \endverbatim 221*> 222*> \param[in] LDAB 223*> \verbatim 224*> LDAB is INTEGER 225*> The leading dimension of AB. It must be at least 2 (not 1!) 226*> and at least max( KK )+1. 227*> \endverbatim 228*> 229*> \param[out] BD 230*> \verbatim 231*> BD is DOUBLE PRECISION array, dimension (max(NN)) 232*> Used to hold the diagonal of the bidiagonal matrix computed 233*> by DGBBRD. 234*> \endverbatim 235*> 236*> \param[out] BE 237*> \verbatim 238*> BE is DOUBLE PRECISION array, dimension (max(NN)) 239*> Used to hold the off-diagonal of the bidiagonal matrix 240*> computed by DGBBRD. 241*> \endverbatim 242*> 243*> \param[out] Q 244*> \verbatim 245*> Q is DOUBLE PRECISION array, dimension (LDQ, max(NN)) 246*> Used to hold the orthogonal matrix Q computed by DGBBRD. 247*> \endverbatim 248*> 249*> \param[in] LDQ 250*> \verbatim 251*> LDQ is INTEGER 252*> The leading dimension of Q. It must be at least 1 253*> and at least max( NN ). 254*> \endverbatim 255*> 256*> \param[out] P 257*> \verbatim 258*> P is DOUBLE PRECISION array, dimension (LDP, max(NN)) 259*> Used to hold the orthogonal matrix P computed by DGBBRD. 260*> \endverbatim 261*> 262*> \param[in] LDP 263*> \verbatim 264*> LDP is INTEGER 265*> The leading dimension of P. It must be at least 1 266*> and at least max( NN ). 267*> \endverbatim 268*> 269*> \param[out] C 270*> \verbatim 271*> C is DOUBLE PRECISION array, dimension (LDC, max(NN)) 272*> Used to hold the matrix C updated by DGBBRD. 273*> \endverbatim 274*> 275*> \param[in] LDC 276*> \verbatim 277*> LDC is INTEGER 278*> The leading dimension of U. It must be at least 1 279*> and at least max( NN ). 280*> \endverbatim 281*> 282*> \param[out] CC 283*> \verbatim 284*> CC is DOUBLE PRECISION array, dimension (LDC, max(NN)) 285*> Used to hold a copy of the matrix C. 286*> \endverbatim 287*> 288*> \param[out] WORK 289*> \verbatim 290*> WORK is DOUBLE PRECISION array, dimension (LWORK) 291*> \endverbatim 292*> 293*> \param[in] LWORK 294*> \verbatim 295*> LWORK is INTEGER 296*> The number of entries in WORK. This must be at least 297*> max( LDA+1, max(NN)+1 )*max(NN). 298*> \endverbatim 299*> 300*> \param[out] RESULT 301*> \verbatim 302*> RESULT is DOUBLE PRECISION array, dimension (4) 303*> The values computed by the tests described above. 304*> The values are currently limited to 1/ulp, to avoid 305*> overflow. 306*> \endverbatim 307*> 308*> \param[out] INFO 309*> \verbatim 310*> INFO is INTEGER 311*> If 0, then everything ran OK. 312*> 313*>----------------------------------------------------------------------- 314*> 315*> Some Local Variables and Parameters: 316*> ---- ----- --------- --- ---------- 317*> ZERO, ONE Real 0 and 1. 318*> MAXTYP The number of types defined. 319*> NTEST The number of tests performed, or which can 320*> be performed so far, for the current matrix. 321*> NTESTT The total number of tests performed so far. 322*> NMAX Largest value in NN. 323*> NMATS The number of matrices generated so far. 324*> NERRS The number of tests which have exceeded THRESH 325*> so far. 326*> COND, IMODE Values to be passed to the matrix generators. 327*> ANORM Norm of A; passed to matrix generators. 328*> 329*> OVFL, UNFL Overflow and underflow thresholds. 330*> ULP, ULPINV Finest relative precision and its inverse. 331*> RTOVFL, RTUNFL Square roots of the previous 2 values. 332*> The following four arrays decode JTYPE: 333*> KTYPE(j) The general type (1-10) for type "j". 334*> KMODE(j) The MODE value to be passed to the matrix 335*> generator for type "j". 336*> KMAGN(j) The order of magnitude ( O(1), 337*> O(overflow^(1/2) ), O(underflow^(1/2) ) 338*> \endverbatim 339* 340* Authors: 341* ======== 342* 343*> \author Univ. of Tennessee 344*> \author Univ. of California Berkeley 345*> \author Univ. of Colorado Denver 346*> \author NAG Ltd. 347* 348*> \date November 2011 349* 350*> \ingroup double_eig 351* 352* ===================================================================== 353 SUBROUTINE DCHKBB( NSIZES, MVAL, NVAL, NWDTHS, KK, NTYPES, DOTYPE, 354 $ NRHS, ISEED, THRESH, NOUNIT, A, LDA, AB, LDAB, 355 $ BD, BE, Q, LDQ, P, LDP, C, LDC, CC, WORK, 356 $ LWORK, RESULT, INFO ) 357* 358* -- LAPACK test routine (input) -- 359* -- LAPACK is a software package provided by Univ. of Tennessee, -- 360* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 361* November 2011 362* 363* .. Scalar Arguments .. 364 INTEGER INFO, LDA, LDAB, LDC, LDP, LDQ, LWORK, NOUNIT, 365 $ NRHS, NSIZES, NTYPES, NWDTHS 366 DOUBLE PRECISION THRESH 367* .. 368* .. Array Arguments .. 369 LOGICAL DOTYPE( * ) 370 INTEGER ISEED( 4 ), KK( * ), MVAL( * ), NVAL( * ) 371 DOUBLE PRECISION A( LDA, * ), AB( LDAB, * ), BD( * ), BE( * ), 372 $ C( LDC, * ), CC( LDC, * ), P( LDP, * ), 373 $ Q( LDQ, * ), RESULT( * ), WORK( * ) 374* .. 375* 376* ===================================================================== 377* 378* .. Parameters .. 379 DOUBLE PRECISION ZERO, ONE 380 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 381 INTEGER MAXTYP 382 PARAMETER ( MAXTYP = 15 ) 383* .. 384* .. Local Scalars .. 385 LOGICAL BADMM, BADNN, BADNNB 386 INTEGER I, IINFO, IMODE, ITYPE, J, JCOL, JR, JSIZE, 387 $ JTYPE, JWIDTH, K, KL, KMAX, KU, M, MMAX, MNMAX, 388 $ MNMIN, MTYPES, N, NERRS, NMATS, NMAX, NTEST, 389 $ NTESTT 390 DOUBLE PRECISION AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL, ULP, 391 $ ULPINV, UNFL 392* .. 393* .. Local Arrays .. 394 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KMAGN( MAXTYP ), 395 $ KMODE( MAXTYP ), KTYPE( MAXTYP ) 396* .. 397* .. External Functions .. 398 DOUBLE PRECISION DLAMCH 399 EXTERNAL DLAMCH 400* .. 401* .. External Subroutines .. 402 EXTERNAL DBDT01, DBDT02, DGBBRD, DLACPY, DLAHD2, DLASET, 403 $ DLASUM, DLATMR, DLATMS, DORT01, XERBLA 404* .. 405* .. Intrinsic Functions .. 406 INTRINSIC ABS, DBLE, MAX, MIN, SQRT 407* .. 408* .. Data statements .. 409 DATA KTYPE / 1, 2, 5*4, 5*6, 3*9 / 410 DATA KMAGN / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3 / 411 DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0, 412 $ 0, 0 / 413* .. 414* .. Executable Statements .. 415* 416* Check for errors 417* 418 NTESTT = 0 419 INFO = 0 420* 421* Important constants 422* 423 BADMM = .FALSE. 424 BADNN = .FALSE. 425 MMAX = 1 426 NMAX = 1 427 MNMAX = 1 428 DO 10 J = 1, NSIZES 429 MMAX = MAX( MMAX, MVAL( J ) ) 430 IF( MVAL( J ).LT.0 ) 431 $ BADMM = .TRUE. 432 NMAX = MAX( NMAX, NVAL( J ) ) 433 IF( NVAL( J ).LT.0 ) 434 $ BADNN = .TRUE. 435 MNMAX = MAX( MNMAX, MIN( MVAL( J ), NVAL( J ) ) ) 436 10 CONTINUE 437* 438 BADNNB = .FALSE. 439 KMAX = 0 440 DO 20 J = 1, NWDTHS 441 KMAX = MAX( KMAX, KK( J ) ) 442 IF( KK( J ).LT.0 ) 443 $ BADNNB = .TRUE. 444 20 CONTINUE 445* 446* Check for errors 447* 448 IF( NSIZES.LT.0 ) THEN 449 INFO = -1 450 ELSE IF( BADMM ) THEN 451 INFO = -2 452 ELSE IF( BADNN ) THEN 453 INFO = -3 454 ELSE IF( NWDTHS.LT.0 ) THEN 455 INFO = -4 456 ELSE IF( BADNNB ) THEN 457 INFO = -5 458 ELSE IF( NTYPES.LT.0 ) THEN 459 INFO = -6 460 ELSE IF( NRHS.LT.0 ) THEN 461 INFO = -8 462 ELSE IF( LDA.LT.NMAX ) THEN 463 INFO = -13 464 ELSE IF( LDAB.LT.2*KMAX+1 ) THEN 465 INFO = -15 466 ELSE IF( LDQ.LT.NMAX ) THEN 467 INFO = -19 468 ELSE IF( LDP.LT.NMAX ) THEN 469 INFO = -21 470 ELSE IF( LDC.LT.NMAX ) THEN 471 INFO = -23 472 ELSE IF( ( MAX( LDA, NMAX )+1 )*NMAX.GT.LWORK ) THEN 473 INFO = -26 474 END IF 475* 476 IF( INFO.NE.0 ) THEN 477 CALL XERBLA( 'DCHKBB', -INFO ) 478 RETURN 479 END IF 480* 481* Quick return if possible 482* 483 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 .OR. NWDTHS.EQ.0 ) 484 $ RETURN 485* 486* More Important constants 487* 488 UNFL = DLAMCH( 'Safe minimum' ) 489 OVFL = ONE / UNFL 490 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 491 ULPINV = ONE / ULP 492 RTUNFL = SQRT( UNFL ) 493 RTOVFL = SQRT( OVFL ) 494* 495* Loop over sizes, widths, types 496* 497 NERRS = 0 498 NMATS = 0 499* 500 DO 160 JSIZE = 1, NSIZES 501 M = MVAL( JSIZE ) 502 N = NVAL( JSIZE ) 503 MNMIN = MIN( M, N ) 504 AMNINV = ONE / DBLE( MAX( 1, M, N ) ) 505* 506 DO 150 JWIDTH = 1, NWDTHS 507 K = KK( JWIDTH ) 508 IF( K.GE.M .AND. K.GE.N ) 509 $ GO TO 150 510 KL = MAX( 0, MIN( M-1, K ) ) 511 KU = MAX( 0, MIN( N-1, K ) ) 512* 513 IF( NSIZES.NE.1 ) THEN 514 MTYPES = MIN( MAXTYP, NTYPES ) 515 ELSE 516 MTYPES = MIN( MAXTYP+1, NTYPES ) 517 END IF 518* 519 DO 140 JTYPE = 1, MTYPES 520 IF( .NOT.DOTYPE( JTYPE ) ) 521 $ GO TO 140 522 NMATS = NMATS + 1 523 NTEST = 0 524* 525 DO 30 J = 1, 4 526 IOLDSD( J ) = ISEED( J ) 527 30 CONTINUE 528* 529* Compute "A". 530* 531* Control parameters: 532* 533* KMAGN KMODE KTYPE 534* =1 O(1) clustered 1 zero 535* =2 large clustered 2 identity 536* =3 small exponential (none) 537* =4 arithmetic diagonal, (w/ singular values) 538* =5 random log (none) 539* =6 random nonhermitian, w/ singular values 540* =7 (none) 541* =8 (none) 542* =9 random nonhermitian 543* 544 IF( MTYPES.GT.MAXTYP ) 545 $ GO TO 90 546* 547 ITYPE = KTYPE( JTYPE ) 548 IMODE = KMODE( JTYPE ) 549* 550* Compute norm 551* 552 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 553* 554 40 CONTINUE 555 ANORM = ONE 556 GO TO 70 557* 558 50 CONTINUE 559 ANORM = ( RTOVFL*ULP )*AMNINV 560 GO TO 70 561* 562 60 CONTINUE 563 ANORM = RTUNFL*MAX( M, N )*ULPINV 564 GO TO 70 565* 566 70 CONTINUE 567* 568 CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 569 CALL DLASET( 'Full', LDAB, N, ZERO, ZERO, AB, LDAB ) 570 IINFO = 0 571 COND = ULPINV 572* 573* Special Matrices -- Identity & Jordan block 574* 575* Zero 576* 577 IF( ITYPE.EQ.1 ) THEN 578 IINFO = 0 579* 580 ELSE IF( ITYPE.EQ.2 ) THEN 581* 582* Identity 583* 584 DO 80 JCOL = 1, N 585 A( JCOL, JCOL ) = ANORM 586 80 CONTINUE 587* 588 ELSE IF( ITYPE.EQ.4 ) THEN 589* 590* Diagonal Matrix, singular values specified 591* 592 CALL DLATMS( M, N, 'S', ISEED, 'N', WORK, IMODE, COND, 593 $ ANORM, 0, 0, 'N', A, LDA, WORK( M+1 ), 594 $ IINFO ) 595* 596 ELSE IF( ITYPE.EQ.6 ) THEN 597* 598* Nonhermitian, singular values specified 599* 600 CALL DLATMS( M, N, 'S', ISEED, 'N', WORK, IMODE, COND, 601 $ ANORM, KL, KU, 'N', A, LDA, WORK( M+1 ), 602 $ IINFO ) 603* 604 ELSE IF( ITYPE.EQ.9 ) THEN 605* 606* Nonhermitian, random entries 607* 608 CALL DLATMR( M, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 609 $ 'T', 'N', WORK( N+1 ), 1, ONE, 610 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, KL, 611 $ KU, ZERO, ANORM, 'N', A, LDA, IDUMMA, 612 $ IINFO ) 613* 614 ELSE 615* 616 IINFO = 1 617 END IF 618* 619* Generate Right-Hand Side 620* 621 CALL DLATMR( M, NRHS, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 622 $ 'T', 'N', WORK( M+1 ), 1, ONE, 623 $ WORK( 2*M+1 ), 1, ONE, 'N', IDUMMA, M, NRHS, 624 $ ZERO, ONE, 'NO', C, LDC, IDUMMA, IINFO ) 625* 626 IF( IINFO.NE.0 ) THEN 627 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, 628 $ JTYPE, IOLDSD 629 INFO = ABS( IINFO ) 630 RETURN 631 END IF 632* 633 90 CONTINUE 634* 635* Copy A to band storage. 636* 637 DO 110 J = 1, N 638 DO 100 I = MAX( 1, J-KU ), MIN( M, J+KL ) 639 AB( KU+1+I-J, J ) = A( I, J ) 640 100 CONTINUE 641 110 CONTINUE 642* 643* Copy C 644* 645 CALL DLACPY( 'Full', M, NRHS, C, LDC, CC, LDC ) 646* 647* Call DGBBRD to compute B, Q and P, and to update C. 648* 649 CALL DGBBRD( 'B', M, N, NRHS, KL, KU, AB, LDAB, BD, BE, 650 $ Q, LDQ, P, LDP, CC, LDC, WORK, IINFO ) 651* 652 IF( IINFO.NE.0 ) THEN 653 WRITE( NOUNIT, FMT = 9999 )'DGBBRD', IINFO, N, JTYPE, 654 $ IOLDSD 655 INFO = ABS( IINFO ) 656 IF( IINFO.LT.0 ) THEN 657 RETURN 658 ELSE 659 RESULT( 1 ) = ULPINV 660 GO TO 120 661 END IF 662 END IF 663* 664* Test 1: Check the decomposition A := Q * B * P' 665* 2: Check the orthogonality of Q 666* 3: Check the orthogonality of P 667* 4: Check the computation of Q' * C 668* 669 CALL DBDT01( M, N, -1, A, LDA, Q, LDQ, BD, BE, P, LDP, 670 $ WORK, RESULT( 1 ) ) 671 CALL DORT01( 'Columns', M, M, Q, LDQ, WORK, LWORK, 672 $ RESULT( 2 ) ) 673 CALL DORT01( 'Rows', N, N, P, LDP, WORK, LWORK, 674 $ RESULT( 3 ) ) 675 CALL DBDT02( M, NRHS, C, LDC, CC, LDC, Q, LDQ, WORK, 676 $ RESULT( 4 ) ) 677* 678* End of Loop -- Check for RESULT(j) > THRESH 679* 680 NTEST = 4 681 120 CONTINUE 682 NTESTT = NTESTT + NTEST 683* 684* Print out tests which fail. 685* 686 DO 130 JR = 1, NTEST 687 IF( RESULT( JR ).GE.THRESH ) THEN 688 IF( NERRS.EQ.0 ) 689 $ CALL DLAHD2( NOUNIT, 'DBB' ) 690 NERRS = NERRS + 1 691 WRITE( NOUNIT, FMT = 9998 )M, N, K, IOLDSD, JTYPE, 692 $ JR, RESULT( JR ) 693 END IF 694 130 CONTINUE 695* 696 140 CONTINUE 697 150 CONTINUE 698 160 CONTINUE 699* 700* Summary 701* 702 CALL DLASUM( 'DBB', NOUNIT, NERRS, NTESTT ) 703 RETURN 704* 705 9999 FORMAT( ' DCHKBB: ', A, ' returned INFO=', I5, '.', / 9X, 'M=', 706 $ I5, ' N=', I5, ' K=', I5, ', JTYPE=', I5, ', ISEED=(', 707 $ 3( I5, ',' ), I5, ')' ) 708 9998 FORMAT( ' M =', I4, ' N=', I4, ', K=', I3, ', seed=', 709 $ 4( I4, ',' ), ' type ', I2, ', test(', I2, ')=', G10.3 ) 710* 711* End of DCHKBB 712* 713 END 714