1*> \brief \b DCHKBD 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 DCHKBD( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS, 12* ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX, 13* Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK, 14* IWORK, NOUT, INFO ) 15* 16* .. Scalar Arguments .. 17* INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS, 18* $ NSIZES, NTYPES 19* DOUBLE PRECISION THRESH 20* .. 21* .. Array Arguments .. 22* LOGICAL DOTYPE( * ) 23* INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * ) 24* DOUBLE PRECISION A( LDA, * ), BD( * ), BE( * ), PT( LDPT, * ), 25* $ Q( LDQ, * ), S1( * ), S2( * ), U( LDPT, * ), 26* $ VT( LDPT, * ), WORK( * ), X( LDX, * ), 27* $ Y( LDX, * ), Z( LDX, * ) 28* .. 29* 30* 31*> \par Purpose: 32* ============= 33*> 34*> \verbatim 35*> 36*> DCHKBD checks the singular value decomposition (SVD) routines. 37*> 38*> DGEBRD reduces a real general m by n matrix A to upper or lower 39*> bidiagonal form B by an orthogonal transformation: Q' * A * P = B 40*> (or A = Q * B * P'). The matrix B is upper bidiagonal if m >= n 41*> and lower bidiagonal if m < n. 42*> 43*> DORGBR generates the orthogonal matrices Q and P' from DGEBRD. 44*> Note that Q and P are not necessarily square. 45*> 46*> DBDSQR computes the singular value decomposition of the bidiagonal 47*> matrix B as B = U S V'. It is called three times to compute 48*> 1) B = U S1 V', where S1 is the diagonal matrix of singular 49*> values and the columns of the matrices U and V are the left 50*> and right singular vectors, respectively, of B. 51*> 2) Same as 1), but the singular values are stored in S2 and the 52*> singular vectors are not computed. 53*> 3) A = (UQ) S (P'V'), the SVD of the original matrix A. 54*> In addition, DBDSQR has an option to apply the left orthogonal matrix 55*> U to a matrix X, useful in least squares applications. 56*> 57*> DBDSDC computes the singular value decomposition of the bidiagonal 58*> matrix B as B = U S V' using divide-and-conquer. It is called twice 59*> to compute 60*> 1) B = U S1 V', where S1 is the diagonal matrix of singular 61*> values and the columns of the matrices U and V are the left 62*> and right singular vectors, respectively, of B. 63*> 2) Same as 1), but the singular values are stored in S2 and the 64*> singular vectors are not computed. 65*> 66*> For each pair of matrix dimensions (M,N) and each selected matrix 67*> type, an M by N matrix A and an M by NRHS matrix X are generated. 68*> The problem dimensions are as follows 69*> A: M x N 70*> Q: M x min(M,N) (but M x M if NRHS > 0) 71*> P: min(M,N) x N 72*> B: min(M,N) x min(M,N) 73*> U, V: min(M,N) x min(M,N) 74*> S1, S2 diagonal, order min(M,N) 75*> X: M x NRHS 76*> 77*> For each generated matrix, 14 tests are performed: 78*> 79*> Test DGEBRD and DORGBR 80*> 81*> (1) | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P' 82*> 83*> (2) | I - Q' Q | / ( M ulp ) 84*> 85*> (3) | I - PT PT' | / ( N ulp ) 86*> 87*> Test DBDSQR on bidiagonal matrix B 88*> 89*> (4) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V' 90*> 91*> (5) | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X 92*> and Z = U' Y. 93*> (6) | I - U' U | / ( min(M,N) ulp ) 94*> 95*> (7) | I - VT VT' | / ( min(M,N) ulp ) 96*> 97*> (8) S1 contains min(M,N) nonnegative values in decreasing order. 98*> (Return 0 if true, 1/ULP if false.) 99*> 100*> (9) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without 101*> computing U and V. 102*> 103*> (10) 0 if the true singular values of B are within THRESH of 104*> those in S1. 2*THRESH if they are not. (Tested using 105*> DSVDCH) 106*> 107*> Test DBDSQR on matrix A 108*> 109*> (11) | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp ) 110*> 111*> (12) | X - (QU) Z | / ( |X| max(M,k) ulp ) 112*> 113*> (13) | I - (QU)'(QU) | / ( M ulp ) 114*> 115*> (14) | I - (VT PT) (PT'VT') | / ( N ulp ) 116*> 117*> Test DBDSDC on bidiagonal matrix B 118*> 119*> (15) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V' 120*> 121*> (16) | I - U' U | / ( min(M,N) ulp ) 122*> 123*> (17) | I - VT VT' | / ( min(M,N) ulp ) 124*> 125*> (18) S1 contains min(M,N) nonnegative values in decreasing order. 126*> (Return 0 if true, 1/ULP if false.) 127*> 128*> (19) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without 129*> computing U and V. 130*> The possible matrix types are 131*> 132*> (1) The zero matrix. 133*> (2) The identity matrix. 134*> 135*> (3) A diagonal matrix with evenly spaced entries 136*> 1, ..., ULP and random signs. 137*> (ULP = (first number larger than 1) - 1 ) 138*> (4) A diagonal matrix with geometrically spaced entries 139*> 1, ..., ULP and random signs. 140*> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 141*> and random signs. 142*> 143*> (6) Same as (3), but multiplied by SQRT( overflow threshold ) 144*> (7) Same as (3), but multiplied by SQRT( underflow threshold ) 145*> 146*> (8) A matrix of the form U D V, where U and V are orthogonal and 147*> D has evenly spaced entries 1, ..., ULP with random signs 148*> on the diagonal. 149*> 150*> (9) A matrix of the form U D V, where U and V are orthogonal and 151*> D has geometrically spaced entries 1, ..., ULP with random 152*> signs on the diagonal. 153*> 154*> (10) A matrix of the form U D V, where U and V are orthogonal and 155*> D has "clustered" entries 1, ULP,..., ULP with random 156*> signs on the diagonal. 157*> 158*> (11) Same as (8), but multiplied by SQRT( overflow threshold ) 159*> (12) Same as (8), but multiplied by SQRT( underflow threshold ) 160*> 161*> (13) Rectangular matrix with random entries chosen from (-1,1). 162*> (14) Same as (13), but multiplied by SQRT( overflow threshold ) 163*> (15) Same as (13), but multiplied by SQRT( underflow threshold ) 164*> 165*> Special case: 166*> (16) A bidiagonal matrix with random entries chosen from a 167*> logarithmic distribution on [ulp^2,ulp^(-2)] (I.e., each 168*> entry is e^x, where x is chosen uniformly on 169*> [ 2 log(ulp), -2 log(ulp) ] .) For *this* type: 170*> (a) DGEBRD is not called to reduce it to bidiagonal form. 171*> (b) the bidiagonal is min(M,N) x min(M,N); if M<N, the 172*> matrix will be lower bidiagonal, otherwise upper. 173*> (c) only tests 5--8 and 14 are performed. 174*> 175*> A subset of the full set of matrix types may be selected through 176*> the logical array DOTYPE. 177*> \endverbatim 178* 179* Arguments: 180* ========== 181* 182*> \param[in] NSIZES 183*> \verbatim 184*> NSIZES is INTEGER 185*> The number of values of M and N contained in the vectors 186*> MVAL and NVAL. The matrix sizes are used in pairs (M,N). 187*> \endverbatim 188*> 189*> \param[in] MVAL 190*> \verbatim 191*> MVAL is INTEGER array, dimension (NM) 192*> The values of the matrix row dimension M. 193*> \endverbatim 194*> 195*> \param[in] NVAL 196*> \verbatim 197*> NVAL is INTEGER array, dimension (NM) 198*> The values of the matrix column dimension N. 199*> \endverbatim 200*> 201*> \param[in] NTYPES 202*> \verbatim 203*> NTYPES is INTEGER 204*> The number of elements in DOTYPE. If it is zero, DCHKBD 205*> does nothing. It must be at least zero. If it is MAXTYP+1 206*> and NSIZES is 1, then an additional type, MAXTYP+1 is 207*> defined, which is to use whatever matrices are in A and B. 208*> This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 209*> DOTYPE(MAXTYP+1) is .TRUE. . 210*> \endverbatim 211*> 212*> \param[in] DOTYPE 213*> \verbatim 214*> DOTYPE is LOGICAL array, dimension (NTYPES) 215*> If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix 216*> of type j will be generated. If NTYPES is smaller than the 217*> maximum number of types defined (PARAMETER MAXTYP), then 218*> types NTYPES+1 through MAXTYP will not be generated. If 219*> NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through 220*> DOTYPE(NTYPES) will be ignored. 221*> \endverbatim 222*> 223*> \param[in] NRHS 224*> \verbatim 225*> NRHS is INTEGER 226*> The number of columns in the "right-hand side" matrices X, Y, 227*> and Z, used in testing DBDSQR. If NRHS = 0, then the 228*> operations on the right-hand side will not be tested. 229*> NRHS must be at least 0. 230*> \endverbatim 231*> 232*> \param[in,out] ISEED 233*> \verbatim 234*> ISEED is INTEGER array, dimension (4) 235*> On entry ISEED specifies the seed of the random number 236*> generator. The array elements should be between 0 and 4095; 237*> if not they will be reduced mod 4096. Also, ISEED(4) must 238*> be odd. The values of ISEED are changed on exit, and can be 239*> used in the next call to DCHKBD to continue the same random 240*> number sequence. 241*> \endverbatim 242*> 243*> \param[in] THRESH 244*> \verbatim 245*> THRESH is DOUBLE PRECISION 246*> The threshold value for the test ratios. A result is 247*> included in the output file if RESULT >= THRESH. To have 248*> every test ratio printed, use THRESH = 0. Note that the 249*> expected value of the test ratios is O(1), so THRESH should 250*> be a reasonably small multiple of 1, e.g., 10 or 100. 251*> \endverbatim 252*> 253*> \param[out] A 254*> \verbatim 255*> A is DOUBLE PRECISION array, dimension (LDA,NMAX) 256*> where NMAX is the maximum value of N in NVAL. 257*> \endverbatim 258*> 259*> \param[in] LDA 260*> \verbatim 261*> LDA is INTEGER 262*> The leading dimension of the array A. LDA >= max(1,MMAX), 263*> where MMAX is the maximum value of M in MVAL. 264*> \endverbatim 265*> 266*> \param[out] BD 267*> \verbatim 268*> BD is DOUBLE PRECISION array, dimension 269*> (max(min(MVAL(j),NVAL(j)))) 270*> \endverbatim 271*> 272*> \param[out] BE 273*> \verbatim 274*> BE is DOUBLE PRECISION array, dimension 275*> (max(min(MVAL(j),NVAL(j)))) 276*> \endverbatim 277*> 278*> \param[out] S1 279*> \verbatim 280*> S1 is DOUBLE PRECISION array, dimension 281*> (max(min(MVAL(j),NVAL(j)))) 282*> \endverbatim 283*> 284*> \param[out] S2 285*> \verbatim 286*> S2 is DOUBLE PRECISION array, dimension 287*> (max(min(MVAL(j),NVAL(j)))) 288*> \endverbatim 289*> 290*> \param[out] X 291*> \verbatim 292*> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 293*> \endverbatim 294*> 295*> \param[in] LDX 296*> \verbatim 297*> LDX is INTEGER 298*> The leading dimension of the arrays X, Y, and Z. 299*> LDX >= max(1,MMAX) 300*> \endverbatim 301*> 302*> \param[out] Y 303*> \verbatim 304*> Y is DOUBLE PRECISION array, dimension (LDX,NRHS) 305*> \endverbatim 306*> 307*> \param[out] Z 308*> \verbatim 309*> Z is DOUBLE PRECISION array, dimension (LDX,NRHS) 310*> \endverbatim 311*> 312*> \param[out] Q 313*> \verbatim 314*> Q is DOUBLE PRECISION array, dimension (LDQ,MMAX) 315*> \endverbatim 316*> 317*> \param[in] LDQ 318*> \verbatim 319*> LDQ is INTEGER 320*> The leading dimension of the array Q. LDQ >= max(1,MMAX). 321*> \endverbatim 322*> 323*> \param[out] PT 324*> \verbatim 325*> PT is DOUBLE PRECISION array, dimension (LDPT,NMAX) 326*> \endverbatim 327*> 328*> \param[in] LDPT 329*> \verbatim 330*> LDPT is INTEGER 331*> The leading dimension of the arrays PT, U, and V. 332*> LDPT >= max(1, max(min(MVAL(j),NVAL(j)))). 333*> \endverbatim 334*> 335*> \param[out] U 336*> \verbatim 337*> U is DOUBLE PRECISION array, dimension 338*> (LDPT,max(min(MVAL(j),NVAL(j)))) 339*> \endverbatim 340*> 341*> \param[out] VT 342*> \verbatim 343*> VT is DOUBLE PRECISION array, dimension 344*> (LDPT,max(min(MVAL(j),NVAL(j)))) 345*> \endverbatim 346*> 347*> \param[out] WORK 348*> \verbatim 349*> WORK is DOUBLE PRECISION array, dimension (LWORK) 350*> \endverbatim 351*> 352*> \param[in] LWORK 353*> \verbatim 354*> LWORK is INTEGER 355*> The number of entries in WORK. This must be at least 356*> 3(M+N) and M(M + max(M,N,k) + 1) + N*min(M,N) for all 357*> pairs (M,N)=(MM(j),NN(j)) 358*> \endverbatim 359*> 360*> \param[out] IWORK 361*> \verbatim 362*> IWORK is INTEGER array, dimension at least 8*min(M,N) 363*> \endverbatim 364*> 365*> \param[in] NOUT 366*> \verbatim 367*> NOUT is INTEGER 368*> The FORTRAN unit number for printing out error messages 369*> (e.g., if a routine returns IINFO not equal to 0.) 370*> \endverbatim 371*> 372*> \param[out] INFO 373*> \verbatim 374*> INFO is INTEGER 375*> If 0, then everything ran OK. 376*> -1: NSIZES < 0 377*> -2: Some MM(j) < 0 378*> -3: Some NN(j) < 0 379*> -4: NTYPES < 0 380*> -6: NRHS < 0 381*> -8: THRESH < 0 382*> -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ). 383*> -17: LDB < 1 or LDB < MMAX. 384*> -21: LDQ < 1 or LDQ < MMAX. 385*> -23: LDPT< 1 or LDPT< MNMAX. 386*> -27: LWORK too small. 387*> If DLATMR, SLATMS, DGEBRD, DORGBR, or DBDSQR, 388*> returns an error code, the 389*> absolute value of it is returned. 390*> 391*>----------------------------------------------------------------------- 392*> 393*> Some Local Variables and Parameters: 394*> ---- ----- --------- --- ---------- 395*> 396*> ZERO, ONE Real 0 and 1. 397*> MAXTYP The number of types defined. 398*> NTEST The number of tests performed, or which can 399*> be performed so far, for the current matrix. 400*> MMAX Largest value in NN. 401*> NMAX Largest value in NN. 402*> MNMIN min(MM(j), NN(j)) (the dimension of the bidiagonal 403*> matrix.) 404*> MNMAX The maximum value of MNMIN for j=1,...,NSIZES. 405*> NFAIL The number of tests which have exceeded THRESH 406*> COND, IMODE Values to be passed to the matrix generators. 407*> ANORM Norm of A; passed to matrix generators. 408*> 409*> OVFL, UNFL Overflow and underflow thresholds. 410*> RTOVFL, RTUNFL Square roots of the previous 2 values. 411*> ULP, ULPINV Finest relative precision and its inverse. 412*> 413*> The following four arrays decode JTYPE: 414*> KTYPE(j) The general type (1-10) for type "j". 415*> KMODE(j) The MODE value to be passed to the matrix 416*> generator for type "j". 417*> KMAGN(j) The order of magnitude ( O(1), 418*> O(overflow^(1/2) ), O(underflow^(1/2) ) 419*> \endverbatim 420* 421* Authors: 422* ======== 423* 424*> \author Univ. of Tennessee 425*> \author Univ. of California Berkeley 426*> \author Univ. of Colorado Denver 427*> \author NAG Ltd. 428* 429*> \date November 2011 430* 431*> \ingroup double_eig 432* 433* ===================================================================== 434 SUBROUTINE DCHKBD( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS, 435 $ ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX, 436 $ Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK, 437 $ IWORK, NOUT, INFO ) 438* 439* -- LAPACK test routine (version 3.4.0) -- 440* -- LAPACK is a software package provided by Univ. of Tennessee, -- 441* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 442* November 2011 443* 444* .. Scalar Arguments .. 445 INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS, 446 $ NSIZES, NTYPES 447 DOUBLE PRECISION THRESH 448* .. 449* .. Array Arguments .. 450 LOGICAL DOTYPE( * ) 451 INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * ) 452 DOUBLE PRECISION A( LDA, * ), BD( * ), BE( * ), PT( LDPT, * ), 453 $ Q( LDQ, * ), S1( * ), S2( * ), U( LDPT, * ), 454 $ VT( LDPT, * ), WORK( * ), X( LDX, * ), 455 $ Y( LDX, * ), Z( LDX, * ) 456* .. 457* 458* ====================================================================== 459* 460* .. Parameters .. 461 DOUBLE PRECISION ZERO, ONE, TWO, HALF 462 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 463 $ HALF = 0.5D0 ) 464 INTEGER MAXTYP 465 PARAMETER ( MAXTYP = 16 ) 466* .. 467* .. Local Scalars .. 468 LOGICAL BADMM, BADNN, BIDIAG 469 CHARACTER UPLO 470 CHARACTER*3 PATH 471 INTEGER I, IINFO, IMODE, ITYPE, J, JCOL, JSIZE, JTYPE, 472 $ LOG2UI, M, MINWRK, MMAX, MNMAX, MNMIN, MQ, 473 $ MTYPES, N, NFAIL, NMAX, NTEST 474 DOUBLE PRECISION AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL, 475 $ TEMP1, TEMP2, ULP, ULPINV, UNFL 476* .. 477* .. Local Arrays .. 478 INTEGER IDUM( 1 ), IOLDSD( 4 ), KMAGN( MAXTYP ), 479 $ KMODE( MAXTYP ), KTYPE( MAXTYP ) 480 DOUBLE PRECISION DUM( 1 ), DUMMA( 1 ), RESULT( 19 ) 481* .. 482* .. External Functions .. 483 DOUBLE PRECISION DLAMCH, DLARND 484 EXTERNAL DLAMCH, DLARND 485* .. 486* .. External Subroutines .. 487 EXTERNAL ALASUM, DBDSDC, DBDSQR, DBDT01, DBDT02, DBDT03, 488 $ DCOPY, DGEBRD, DGEMM, DLABAD, DLACPY, DLAHD2, 489 $ DLASET, DLATMR, DLATMS, DORGBR, DORT01, XERBLA 490* .. 491* .. Intrinsic Functions .. 492 INTRINSIC ABS, EXP, INT, LOG, MAX, MIN, SQRT 493* .. 494* .. Scalars in Common .. 495 LOGICAL LERR, OK 496 CHARACTER*32 SRNAMT 497 INTEGER INFOT, NUNIT 498* .. 499* .. Common blocks .. 500 COMMON / INFOC / INFOT, NUNIT, OK, LERR 501 COMMON / SRNAMC / SRNAMT 502* .. 503* .. Data statements .. 504 DATA KTYPE / 1, 2, 5*4, 5*6, 3*9, 10 / 505 DATA KMAGN / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 / 506 DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0, 507 $ 0, 0, 0 / 508* .. 509* .. Executable Statements .. 510* 511* Check for errors 512* 513 INFO = 0 514* 515 BADMM = .FALSE. 516 BADNN = .FALSE. 517 MMAX = 1 518 NMAX = 1 519 MNMAX = 1 520 MINWRK = 1 521 DO 10 J = 1, NSIZES 522 MMAX = MAX( MMAX, MVAL( J ) ) 523 IF( MVAL( J ).LT.0 ) 524 $ BADMM = .TRUE. 525 NMAX = MAX( NMAX, NVAL( J ) ) 526 IF( NVAL( J ).LT.0 ) 527 $ BADNN = .TRUE. 528 MNMAX = MAX( MNMAX, MIN( MVAL( J ), NVAL( J ) ) ) 529 MINWRK = MAX( MINWRK, 3*( MVAL( J )+NVAL( J ) ), 530 $ MVAL( J )*( MVAL( J )+MAX( MVAL( J ), NVAL( J ), 531 $ NRHS )+1 )+NVAL( J )*MIN( NVAL( J ), MVAL( J ) ) ) 532 10 CONTINUE 533* 534* Check for errors 535* 536 IF( NSIZES.LT.0 ) THEN 537 INFO = -1 538 ELSE IF( BADMM ) THEN 539 INFO = -2 540 ELSE IF( BADNN ) THEN 541 INFO = -3 542 ELSE IF( NTYPES.LT.0 ) THEN 543 INFO = -4 544 ELSE IF( NRHS.LT.0 ) THEN 545 INFO = -6 546 ELSE IF( LDA.LT.MMAX ) THEN 547 INFO = -11 548 ELSE IF( LDX.LT.MMAX ) THEN 549 INFO = -17 550 ELSE IF( LDQ.LT.MMAX ) THEN 551 INFO = -21 552 ELSE IF( LDPT.LT.MNMAX ) THEN 553 INFO = -23 554 ELSE IF( MINWRK.GT.LWORK ) THEN 555 INFO = -27 556 END IF 557* 558 IF( INFO.NE.0 ) THEN 559 CALL XERBLA( 'DCHKBD', -INFO ) 560 RETURN 561 END IF 562* 563* Initialize constants 564* 565 PATH( 1: 1 ) = 'Double precision' 566 PATH( 2: 3 ) = 'BD' 567 NFAIL = 0 568 NTEST = 0 569 UNFL = DLAMCH( 'Safe minimum' ) 570 OVFL = DLAMCH( 'Overflow' ) 571 CALL DLABAD( UNFL, OVFL ) 572 ULP = DLAMCH( 'Precision' ) 573 ULPINV = ONE / ULP 574 LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) ) 575 RTUNFL = SQRT( UNFL ) 576 RTOVFL = SQRT( OVFL ) 577 INFOT = 0 578* 579* Loop over sizes, types 580* 581 DO 200 JSIZE = 1, NSIZES 582 M = MVAL( JSIZE ) 583 N = NVAL( JSIZE ) 584 MNMIN = MIN( M, N ) 585 AMNINV = ONE / MAX( M, N, 1 ) 586* 587 IF( NSIZES.NE.1 ) THEN 588 MTYPES = MIN( MAXTYP, NTYPES ) 589 ELSE 590 MTYPES = MIN( MAXTYP+1, NTYPES ) 591 END IF 592* 593 DO 190 JTYPE = 1, MTYPES 594 IF( .NOT.DOTYPE( JTYPE ) ) 595 $ GO TO 190 596* 597 DO 20 J = 1, 4 598 IOLDSD( J ) = ISEED( J ) 599 20 CONTINUE 600* 601 DO 30 J = 1, 14 602 RESULT( J ) = -ONE 603 30 CONTINUE 604* 605 UPLO = ' ' 606* 607* Compute "A" 608* 609* Control parameters: 610* 611* KMAGN KMODE KTYPE 612* =1 O(1) clustered 1 zero 613* =2 large clustered 2 identity 614* =3 small exponential (none) 615* =4 arithmetic diagonal, (w/ eigenvalues) 616* =5 random symmetric, w/ eigenvalues 617* =6 nonsymmetric, w/ singular values 618* =7 random diagonal 619* =8 random symmetric 620* =9 random nonsymmetric 621* =10 random bidiagonal (log. distrib.) 622* 623 IF( MTYPES.GT.MAXTYP ) 624 $ GO TO 100 625* 626 ITYPE = KTYPE( JTYPE ) 627 IMODE = KMODE( JTYPE ) 628* 629* Compute norm 630* 631 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 632* 633 40 CONTINUE 634 ANORM = ONE 635 GO TO 70 636* 637 50 CONTINUE 638 ANORM = ( RTOVFL*ULP )*AMNINV 639 GO TO 70 640* 641 60 CONTINUE 642 ANORM = RTUNFL*MAX( M, N )*ULPINV 643 GO TO 70 644* 645 70 CONTINUE 646* 647 CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 648 IINFO = 0 649 COND = ULPINV 650* 651 BIDIAG = .FALSE. 652 IF( ITYPE.EQ.1 ) THEN 653* 654* Zero matrix 655* 656 IINFO = 0 657* 658 ELSE IF( ITYPE.EQ.2 ) THEN 659* 660* Identity 661* 662 DO 80 JCOL = 1, MNMIN 663 A( JCOL, JCOL ) = ANORM 664 80 CONTINUE 665* 666 ELSE IF( ITYPE.EQ.4 ) THEN 667* 668* Diagonal Matrix, [Eigen]values Specified 669* 670 CALL DLATMS( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, IMODE, 671 $ COND, ANORM, 0, 0, 'N', A, LDA, 672 $ WORK( MNMIN+1 ), IINFO ) 673* 674 ELSE IF( ITYPE.EQ.5 ) THEN 675* 676* Symmetric, eigenvalues specified 677* 678 CALL DLATMS( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, IMODE, 679 $ COND, ANORM, M, N, 'N', A, LDA, 680 $ WORK( MNMIN+1 ), IINFO ) 681* 682 ELSE IF( ITYPE.EQ.6 ) THEN 683* 684* Nonsymmetric, singular values specified 685* 686 CALL DLATMS( M, N, 'S', ISEED, 'N', WORK, IMODE, COND, 687 $ ANORM, M, N, 'N', A, LDA, WORK( MNMIN+1 ), 688 $ IINFO ) 689* 690 ELSE IF( ITYPE.EQ.7 ) THEN 691* 692* Diagonal, random entries 693* 694 CALL DLATMR( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, 6, ONE, 695 $ ONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE, 696 $ WORK( 2*MNMIN+1 ), 1, ONE, 'N', IWORK, 0, 0, 697 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 698* 699 ELSE IF( ITYPE.EQ.8 ) THEN 700* 701* Symmetric, random entries 702* 703 CALL DLATMR( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, 6, ONE, 704 $ ONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE, 705 $ WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N, 706 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 707* 708 ELSE IF( ITYPE.EQ.9 ) THEN 709* 710* Nonsymmetric, random entries 711* 712 CALL DLATMR( M, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 713 $ 'T', 'N', WORK( MNMIN+1 ), 1, ONE, 714 $ WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N, 715 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 716* 717 ELSE IF( ITYPE.EQ.10 ) THEN 718* 719* Bidiagonal, random entries 720* 721 TEMP1 = -TWO*LOG( ULP ) 722 DO 90 J = 1, MNMIN 723 BD( J ) = EXP( TEMP1*DLARND( 2, ISEED ) ) 724 IF( J.LT.MNMIN ) 725 $ BE( J ) = EXP( TEMP1*DLARND( 2, ISEED ) ) 726 90 CONTINUE 727* 728 IINFO = 0 729 BIDIAG = .TRUE. 730 IF( M.GE.N ) THEN 731 UPLO = 'U' 732 ELSE 733 UPLO = 'L' 734 END IF 735 ELSE 736 IINFO = 1 737 END IF 738* 739 IF( IINFO.EQ.0 ) THEN 740* 741* Generate Right-Hand Side 742* 743 IF( BIDIAG ) THEN 744 CALL DLATMR( MNMIN, NRHS, 'S', ISEED, 'N', WORK, 6, 745 $ ONE, ONE, 'T', 'N', WORK( MNMIN+1 ), 1, 746 $ ONE, WORK( 2*MNMIN+1 ), 1, ONE, 'N', 747 $ IWORK, MNMIN, NRHS, ZERO, ONE, 'NO', Y, 748 $ LDX, IWORK, IINFO ) 749 ELSE 750 CALL DLATMR( M, NRHS, 'S', ISEED, 'N', WORK, 6, ONE, 751 $ ONE, 'T', 'N', WORK( M+1 ), 1, ONE, 752 $ WORK( 2*M+1 ), 1, ONE, 'N', IWORK, M, 753 $ NRHS, ZERO, ONE, 'NO', X, LDX, IWORK, 754 $ IINFO ) 755 END IF 756 END IF 757* 758* Error Exit 759* 760 IF( IINFO.NE.0 ) THEN 761 WRITE( NOUT, FMT = 9998 )'Generator', IINFO, M, N, 762 $ JTYPE, IOLDSD 763 INFO = ABS( IINFO ) 764 RETURN 765 END IF 766* 767 100 CONTINUE 768* 769* Call DGEBRD and DORGBR to compute B, Q, and P, do tests. 770* 771 IF( .NOT.BIDIAG ) THEN 772* 773* Compute transformations to reduce A to bidiagonal form: 774* B := Q' * A * P. 775* 776 CALL DLACPY( ' ', M, N, A, LDA, Q, LDQ ) 777 CALL DGEBRD( M, N, Q, LDQ, BD, BE, WORK, WORK( MNMIN+1 ), 778 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO ) 779* 780* Check error code from DGEBRD. 781* 782 IF( IINFO.NE.0 ) THEN 783 WRITE( NOUT, FMT = 9998 )'DGEBRD', IINFO, M, N, 784 $ JTYPE, IOLDSD 785 INFO = ABS( IINFO ) 786 RETURN 787 END IF 788* 789 CALL DLACPY( ' ', M, N, Q, LDQ, PT, LDPT ) 790 IF( M.GE.N ) THEN 791 UPLO = 'U' 792 ELSE 793 UPLO = 'L' 794 END IF 795* 796* Generate Q 797* 798 MQ = M 799 IF( NRHS.LE.0 ) 800 $ MQ = MNMIN 801 CALL DORGBR( 'Q', M, MQ, N, Q, LDQ, WORK, 802 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO ) 803* 804* Check error code from DORGBR. 805* 806 IF( IINFO.NE.0 ) THEN 807 WRITE( NOUT, FMT = 9998 )'DORGBR(Q)', IINFO, M, N, 808 $ JTYPE, IOLDSD 809 INFO = ABS( IINFO ) 810 RETURN 811 END IF 812* 813* Generate P' 814* 815 CALL DORGBR( 'P', MNMIN, N, M, PT, LDPT, WORK( MNMIN+1 ), 816 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO ) 817* 818* Check error code from DORGBR. 819* 820 IF( IINFO.NE.0 ) THEN 821 WRITE( NOUT, FMT = 9998 )'DORGBR(P)', IINFO, M, N, 822 $ JTYPE, IOLDSD 823 INFO = ABS( IINFO ) 824 RETURN 825 END IF 826* 827* Apply Q' to an M by NRHS matrix X: Y := Q' * X. 828* 829 CALL DGEMM( 'Transpose', 'No transpose', M, NRHS, M, ONE, 830 $ Q, LDQ, X, LDX, ZERO, Y, LDX ) 831* 832* Test 1: Check the decomposition A := Q * B * PT 833* 2: Check the orthogonality of Q 834* 3: Check the orthogonality of PT 835* 836 CALL DBDT01( M, N, 1, A, LDA, Q, LDQ, BD, BE, PT, LDPT, 837 $ WORK, RESULT( 1 ) ) 838 CALL DORT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK, 839 $ RESULT( 2 ) ) 840 CALL DORT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK, 841 $ RESULT( 3 ) ) 842 END IF 843* 844* Use DBDSQR to form the SVD of the bidiagonal matrix B: 845* B := U * S1 * VT, and compute Z = U' * Y. 846* 847 CALL DCOPY( MNMIN, BD, 1, S1, 1 ) 848 IF( MNMIN.GT.0 ) 849 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 850 CALL DLACPY( ' ', M, NRHS, Y, LDX, Z, LDX ) 851 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, U, LDPT ) 852 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, VT, LDPT ) 853* 854 CALL DBDSQR( UPLO, MNMIN, MNMIN, MNMIN, NRHS, S1, WORK, VT, 855 $ LDPT, U, LDPT, Z, LDX, WORK( MNMIN+1 ), IINFO ) 856* 857* Check error code from DBDSQR. 858* 859 IF( IINFO.NE.0 ) THEN 860 WRITE( NOUT, FMT = 9998 )'DBDSQR(vects)', IINFO, M, N, 861 $ JTYPE, IOLDSD 862 INFO = ABS( IINFO ) 863 IF( IINFO.LT.0 ) THEN 864 RETURN 865 ELSE 866 RESULT( 4 ) = ULPINV 867 GO TO 170 868 END IF 869 END IF 870* 871* Use DBDSQR to compute only the singular values of the 872* bidiagonal matrix B; U, VT, and Z should not be modified. 873* 874 CALL DCOPY( MNMIN, BD, 1, S2, 1 ) 875 IF( MNMIN.GT.0 ) 876 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 877* 878 CALL DBDSQR( UPLO, MNMIN, 0, 0, 0, S2, WORK, VT, LDPT, U, 879 $ LDPT, Z, LDX, WORK( MNMIN+1 ), IINFO ) 880* 881* Check error code from DBDSQR. 882* 883 IF( IINFO.NE.0 ) THEN 884 WRITE( NOUT, FMT = 9998 )'DBDSQR(values)', IINFO, M, N, 885 $ JTYPE, IOLDSD 886 INFO = ABS( IINFO ) 887 IF( IINFO.LT.0 ) THEN 888 RETURN 889 ELSE 890 RESULT( 9 ) = ULPINV 891 GO TO 170 892 END IF 893 END IF 894* 895* Test 4: Check the decomposition B := U * S1 * VT 896* 5: Check the computation Z := U' * Y 897* 6: Check the orthogonality of U 898* 7: Check the orthogonality of VT 899* 900 CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT, 901 $ WORK, RESULT( 4 ) ) 902 CALL DBDT02( MNMIN, NRHS, Y, LDX, Z, LDX, U, LDPT, WORK, 903 $ RESULT( 5 ) ) 904 CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK, 905 $ RESULT( 6 ) ) 906 CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK, 907 $ RESULT( 7 ) ) 908* 909* Test 8: Check that the singular values are sorted in 910* non-increasing order and are non-negative 911* 912 RESULT( 8 ) = ZERO 913 DO 110 I = 1, MNMIN - 1 914 IF( S1( I ).LT.S1( I+1 ) ) 915 $ RESULT( 8 ) = ULPINV 916 IF( S1( I ).LT.ZERO ) 917 $ RESULT( 8 ) = ULPINV 918 110 CONTINUE 919 IF( MNMIN.GE.1 ) THEN 920 IF( S1( MNMIN ).LT.ZERO ) 921 $ RESULT( 8 ) = ULPINV 922 END IF 923* 924* Test 9: Compare DBDSQR with and without singular vectors 925* 926 TEMP2 = ZERO 927* 928 DO 120 J = 1, MNMIN 929 TEMP1 = ABS( S1( J )-S2( J ) ) / 930 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ), 931 $ ULP*MAX( ABS( S1( J ) ), ABS( S2( J ) ) ) ) 932 TEMP2 = MAX( TEMP1, TEMP2 ) 933 120 CONTINUE 934* 935 RESULT( 9 ) = TEMP2 936* 937* Test 10: Sturm sequence test of singular values 938* Go up by factors of two until it succeeds 939* 940 TEMP1 = THRESH*( HALF-ULP ) 941* 942 DO 130 J = 0, LOG2UI 943* CALL DSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO ) 944 IF( IINFO.EQ.0 ) 945 $ GO TO 140 946 TEMP1 = TEMP1*TWO 947 130 CONTINUE 948* 949 140 CONTINUE 950 RESULT( 10 ) = TEMP1 951* 952* Use DBDSQR to form the decomposition A := (QU) S (VT PT) 953* from the bidiagonal form A := Q B PT. 954* 955 IF( .NOT.BIDIAG ) THEN 956 CALL DCOPY( MNMIN, BD, 1, S2, 1 ) 957 IF( MNMIN.GT.0 ) 958 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 959* 960 CALL DBDSQR( UPLO, MNMIN, N, M, NRHS, S2, WORK, PT, LDPT, 961 $ Q, LDQ, Y, LDX, WORK( MNMIN+1 ), IINFO ) 962* 963* Test 11: Check the decomposition A := Q*U * S2 * VT*PT 964* 12: Check the computation Z := U' * Q' * X 965* 13: Check the orthogonality of Q*U 966* 14: Check the orthogonality of VT*PT 967* 968 CALL DBDT01( M, N, 0, A, LDA, Q, LDQ, S2, DUMMA, PT, 969 $ LDPT, WORK, RESULT( 11 ) ) 970 CALL DBDT02( M, NRHS, X, LDX, Y, LDX, Q, LDQ, WORK, 971 $ RESULT( 12 ) ) 972 CALL DORT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK, 973 $ RESULT( 13 ) ) 974 CALL DORT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK, 975 $ RESULT( 14 ) ) 976 END IF 977* 978* Use DBDSDC to form the SVD of the bidiagonal matrix B: 979* B := U * S1 * VT 980* 981 CALL DCOPY( MNMIN, BD, 1, S1, 1 ) 982 IF( MNMIN.GT.0 ) 983 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 984 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, U, LDPT ) 985 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, VT, LDPT ) 986* 987 CALL DBDSDC( UPLO, 'I', MNMIN, S1, WORK, U, LDPT, VT, LDPT, 988 $ DUM, IDUM, WORK( MNMIN+1 ), IWORK, IINFO ) 989* 990* Check error code from DBDSDC. 991* 992 IF( IINFO.NE.0 ) THEN 993 WRITE( NOUT, FMT = 9998 )'DBDSDC(vects)', IINFO, M, N, 994 $ JTYPE, IOLDSD 995 INFO = ABS( IINFO ) 996 IF( IINFO.LT.0 ) THEN 997 RETURN 998 ELSE 999 RESULT( 15 ) = ULPINV 1000 GO TO 170 1001 END IF 1002 END IF 1003* 1004* Use DBDSDC to compute only the singular values of the 1005* bidiagonal matrix B; U and VT should not be modified. 1006* 1007 CALL DCOPY( MNMIN, BD, 1, S2, 1 ) 1008 IF( MNMIN.GT.0 ) 1009 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 1010* 1011 CALL DBDSDC( UPLO, 'N', MNMIN, S2, WORK, DUM, 1, DUM, 1, 1012 $ DUM, IDUM, WORK( MNMIN+1 ), IWORK, IINFO ) 1013* 1014* Check error code from DBDSDC. 1015* 1016 IF( IINFO.NE.0 ) THEN 1017 WRITE( NOUT, FMT = 9998 )'DBDSDC(values)', IINFO, M, N, 1018 $ JTYPE, IOLDSD 1019 INFO = ABS( IINFO ) 1020 IF( IINFO.LT.0 ) THEN 1021 RETURN 1022 ELSE 1023 RESULT( 18 ) = ULPINV 1024 GO TO 170 1025 END IF 1026 END IF 1027* 1028* Test 15: Check the decomposition B := U * S1 * VT 1029* 16: Check the orthogonality of U 1030* 17: Check the orthogonality of VT 1031* 1032 CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT, 1033 $ WORK, RESULT( 15 ) ) 1034 CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK, 1035 $ RESULT( 16 ) ) 1036 CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK, 1037 $ RESULT( 17 ) ) 1038* 1039* Test 18: Check that the singular values are sorted in 1040* non-increasing order and are non-negative 1041* 1042 RESULT( 18 ) = ZERO 1043 DO 150 I = 1, MNMIN - 1 1044 IF( S1( I ).LT.S1( I+1 ) ) 1045 $ RESULT( 18 ) = ULPINV 1046 IF( S1( I ).LT.ZERO ) 1047 $ RESULT( 18 ) = ULPINV 1048 150 CONTINUE 1049 IF( MNMIN.GE.1 ) THEN 1050 IF( S1( MNMIN ).LT.ZERO ) 1051 $ RESULT( 18 ) = ULPINV 1052 END IF 1053* 1054* Test 19: Compare DBDSQR with and without singular vectors 1055* 1056 TEMP2 = ZERO 1057* 1058 DO 160 J = 1, MNMIN 1059 TEMP1 = ABS( S1( J )-S2( J ) ) / 1060 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ), 1061 $ ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) ) 1062 TEMP2 = MAX( TEMP1, TEMP2 ) 1063 160 CONTINUE 1064* 1065 RESULT( 19 ) = TEMP2 1066* 1067* End of Loop -- Check for RESULT(j) > THRESH 1068* 1069 170 CONTINUE 1070 DO 180 J = 1, 19 1071 IF( RESULT( J ).GE.THRESH ) THEN 1072 IF( NFAIL.EQ.0 ) 1073 $ CALL DLAHD2( NOUT, PATH ) 1074 WRITE( NOUT, FMT = 9999 )M, N, JTYPE, IOLDSD, J, 1075 $ RESULT( J ) 1076 NFAIL = NFAIL + 1 1077 END IF 1078 180 CONTINUE 1079 IF( .NOT.BIDIAG ) THEN 1080 NTEST = NTEST + 19 1081 ELSE 1082 NTEST = NTEST + 5 1083 END IF 1084* 1085 190 CONTINUE 1086 200 CONTINUE 1087* 1088* Summary 1089* 1090 CALL ALASUM( PATH, NOUT, NFAIL, NTEST, 0 ) 1091* 1092 RETURN 1093* 1094* End of DCHKBD 1095* 1096 9999 FORMAT( ' M=', I5, ', N=', I5, ', type ', I2, ', seed=', 1097 $ 4( I4, ',' ), ' test(', I2, ')=', G11.4 ) 1098 9998 FORMAT( ' DCHKBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=', 1099 $ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), 1100 $ I5, ')' ) 1101* 1102 END 1103