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*> DBDSVDX computes the singular value decomposition of the bidiagonal 67*> matrix B as B = U S V' using bisection and inverse iteration. It is 68*> called six times to compute 69*> 1) B = U S1 V', RANGE='A', where S1 is the diagonal matrix of singular 70*> values and the columns of the matrices U and V are the left 71*> and right singular vectors, respectively, of B. 72*> 2) Same as 1), but the singular values are stored in S2 and the 73*> singular vectors are not computed. 74*> 3) B = U S1 V', RANGE='I', with where S1 is the diagonal matrix of singular 75*> values and the columns of the matrices U and V are the left 76*> and right singular vectors, respectively, of B 77*> 4) Same as 3), but the singular values are stored in S2 and the 78*> singular vectors are not computed. 79*> 5) B = U S1 V', RANGE='V', with where S1 is the diagonal matrix of singular 80*> values and the columns of the matrices U and V are the left 81*> and right singular vectors, respectively, of B 82*> 6) Same as 5), but the singular values are stored in S2 and the 83*> singular vectors are not computed. 84*> 85*> For each pair of matrix dimensions (M,N) and each selected matrix 86*> type, an M by N matrix A and an M by NRHS matrix X are generated. 87*> The problem dimensions are as follows 88*> A: M x N 89*> Q: M x min(M,N) (but M x M if NRHS > 0) 90*> P: min(M,N) x N 91*> B: min(M,N) x min(M,N) 92*> U, V: min(M,N) x min(M,N) 93*> S1, S2 diagonal, order min(M,N) 94*> X: M x NRHS 95*> 96*> For each generated matrix, 14 tests are performed: 97*> 98*> Test DGEBRD and DORGBR 99*> 100*> (1) | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P' 101*> 102*> (2) | I - Q' Q | / ( M ulp ) 103*> 104*> (3) | I - PT PT' | / ( N ulp ) 105*> 106*> Test DBDSQR on bidiagonal matrix B 107*> 108*> (4) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V' 109*> 110*> (5) | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X 111*> and Z = U' Y. 112*> (6) | I - U' U | / ( min(M,N) ulp ) 113*> 114*> (7) | I - VT VT' | / ( min(M,N) ulp ) 115*> 116*> (8) S1 contains min(M,N) nonnegative values in decreasing order. 117*> (Return 0 if true, 1/ULP if false.) 118*> 119*> (9) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without 120*> computing U and V. 121*> 122*> (10) 0 if the true singular values of B are within THRESH of 123*> those in S1. 2*THRESH if they are not. (Tested using 124*> DSVDCH) 125*> 126*> Test DBDSQR on matrix A 127*> 128*> (11) | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp ) 129*> 130*> (12) | X - (QU) Z | / ( |X| max(M,k) ulp ) 131*> 132*> (13) | I - (QU)'(QU) | / ( M ulp ) 133*> 134*> (14) | I - (VT PT) (PT'VT') | / ( N ulp ) 135*> 136*> Test DBDSDC on bidiagonal matrix B 137*> 138*> (15) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V' 139*> 140*> (16) | I - U' U | / ( min(M,N) ulp ) 141*> 142*> (17) | I - VT VT' | / ( min(M,N) ulp ) 143*> 144*> (18) S1 contains min(M,N) nonnegative values in decreasing order. 145*> (Return 0 if true, 1/ULP if false.) 146*> 147*> (19) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without 148*> computing U and V. 149*> Test DBDSVDX on bidiagonal matrix B 150*> 151*> (20) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V' 152*> 153*> (21) | I - U' U | / ( min(M,N) ulp ) 154*> 155*> (22) | I - VT VT' | / ( min(M,N) ulp ) 156*> 157*> (23) S1 contains min(M,N) nonnegative values in decreasing order. 158*> (Return 0 if true, 1/ULP if false.) 159*> 160*> (24) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without 161*> computing U and V. 162*> 163*> (25) | S1 - U' B VT' | / ( |S| n ulp ) DBDSVDX('V', 'I') 164*> 165*> (26) | I - U' U | / ( min(M,N) ulp ) 166*> 167*> (27) | I - VT VT' | / ( min(M,N) ulp ) 168*> 169*> (28) S1 contains min(M,N) nonnegative values in decreasing order. 170*> (Return 0 if true, 1/ULP if false.) 171*> 172*> (29) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without 173*> computing U and V. 174*> 175*> (30) | S1 - U' B VT' | / ( |S1| n ulp ) DBDSVDX('V', 'V') 176*> 177*> (31) | I - U' U | / ( min(M,N) ulp ) 178*> 179*> (32) | I - VT VT' | / ( min(M,N) ulp ) 180*> 181*> (33) S1 contains min(M,N) nonnegative values in decreasing order. 182*> (Return 0 if true, 1/ULP if false.) 183*> 184*> (34) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without 185*> computing U and V. 186*> 187*> The possible matrix types are 188*> 189*> (1) The zero matrix. 190*> (2) The identity matrix. 191*> 192*> (3) A diagonal matrix with evenly spaced entries 193*> 1, ..., ULP and random signs. 194*> (ULP = (first number larger than 1) - 1 ) 195*> (4) A diagonal matrix with geometrically spaced entries 196*> 1, ..., ULP and random signs. 197*> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 198*> and random signs. 199*> 200*> (6) Same as (3), but multiplied by SQRT( overflow threshold ) 201*> (7) Same as (3), but multiplied by SQRT( underflow threshold ) 202*> 203*> (8) A matrix of the form U D V, where U and V are orthogonal and 204*> D has evenly spaced entries 1, ..., ULP with random signs 205*> on the diagonal. 206*> 207*> (9) A matrix of the form U D V, where U and V are orthogonal and 208*> D has geometrically spaced entries 1, ..., ULP with random 209*> signs on the diagonal. 210*> 211*> (10) A matrix of the form U D V, where U and V are orthogonal and 212*> D has "clustered" entries 1, ULP,..., ULP with random 213*> signs on the diagonal. 214*> 215*> (11) Same as (8), but multiplied by SQRT( overflow threshold ) 216*> (12) Same as (8), but multiplied by SQRT( underflow threshold ) 217*> 218*> (13) Rectangular matrix with random entries chosen from (-1,1). 219*> (14) Same as (13), but multiplied by SQRT( overflow threshold ) 220*> (15) Same as (13), but multiplied by SQRT( underflow threshold ) 221*> 222*> Special case: 223*> (16) A bidiagonal matrix with random entries chosen from a 224*> logarithmic distribution on [ulp^2,ulp^(-2)] (I.e., each 225*> entry is e^x, where x is chosen uniformly on 226*> [ 2 log(ulp), -2 log(ulp) ] .) For *this* type: 227*> (a) DGEBRD is not called to reduce it to bidiagonal form. 228*> (b) the bidiagonal is min(M,N) x min(M,N); if M<N, the 229*> matrix will be lower bidiagonal, otherwise upper. 230*> (c) only tests 5--8 and 14 are performed. 231*> 232*> A subset of the full set of matrix types may be selected through 233*> the logical array DOTYPE. 234*> \endverbatim 235* 236* Arguments: 237* ========== 238* 239*> \param[in] NSIZES 240*> \verbatim 241*> NSIZES is INTEGER 242*> The number of values of M and N contained in the vectors 243*> MVAL and NVAL. The matrix sizes are used in pairs (M,N). 244*> \endverbatim 245*> 246*> \param[in] MVAL 247*> \verbatim 248*> MVAL is INTEGER array, dimension (NM) 249*> The values of the matrix row dimension M. 250*> \endverbatim 251*> 252*> \param[in] NVAL 253*> \verbatim 254*> NVAL is INTEGER array, dimension (NM) 255*> The values of the matrix column dimension N. 256*> \endverbatim 257*> 258*> \param[in] NTYPES 259*> \verbatim 260*> NTYPES is INTEGER 261*> The number of elements in DOTYPE. If it is zero, DCHKBD 262*> does nothing. It must be at least zero. If it is MAXTYP+1 263*> and NSIZES is 1, then an additional type, MAXTYP+1 is 264*> defined, which is to use whatever matrices are in A and B. 265*> This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 266*> DOTYPE(MAXTYP+1) is .TRUE. . 267*> \endverbatim 268*> 269*> \param[in] DOTYPE 270*> \verbatim 271*> DOTYPE is LOGICAL array, dimension (NTYPES) 272*> If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix 273*> of type j will be generated. If NTYPES is smaller than the 274*> maximum number of types defined (PARAMETER MAXTYP), then 275*> types NTYPES+1 through MAXTYP will not be generated. If 276*> NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through 277*> DOTYPE(NTYPES) will be ignored. 278*> \endverbatim 279*> 280*> \param[in] NRHS 281*> \verbatim 282*> NRHS is INTEGER 283*> The number of columns in the "right-hand side" matrices X, Y, 284*> and Z, used in testing DBDSQR. If NRHS = 0, then the 285*> operations on the right-hand side will not be tested. 286*> NRHS must be at least 0. 287*> \endverbatim 288*> 289*> \param[in,out] ISEED 290*> \verbatim 291*> ISEED is INTEGER array, dimension (4) 292*> On entry ISEED specifies the seed of the random number 293*> generator. The array elements should be between 0 and 4095; 294*> if not they will be reduced mod 4096. Also, ISEED(4) must 295*> be odd. The values of ISEED are changed on exit, and can be 296*> used in the next call to DCHKBD to continue the same random 297*> number sequence. 298*> \endverbatim 299*> 300*> \param[in] THRESH 301*> \verbatim 302*> THRESH is DOUBLE PRECISION 303*> The threshold value for the test ratios. A result is 304*> included in the output file if RESULT >= THRESH. To have 305*> every test ratio printed, use THRESH = 0. Note that the 306*> expected value of the test ratios is O(1), so THRESH should 307*> be a reasonably small multiple of 1, e.g., 10 or 100. 308*> \endverbatim 309*> 310*> \param[out] A 311*> \verbatim 312*> A is DOUBLE PRECISION array, dimension (LDA,NMAX) 313*> where NMAX is the maximum value of N in NVAL. 314*> \endverbatim 315*> 316*> \param[in] LDA 317*> \verbatim 318*> LDA is INTEGER 319*> The leading dimension of the array A. LDA >= max(1,MMAX), 320*> where MMAX is the maximum value of M in MVAL. 321*> \endverbatim 322*> 323*> \param[out] BD 324*> \verbatim 325*> BD is DOUBLE PRECISION array, dimension 326*> (max(min(MVAL(j),NVAL(j)))) 327*> \endverbatim 328*> 329*> \param[out] BE 330*> \verbatim 331*> BE is DOUBLE PRECISION array, dimension 332*> (max(min(MVAL(j),NVAL(j)))) 333*> \endverbatim 334*> 335*> \param[out] S1 336*> \verbatim 337*> S1 is DOUBLE PRECISION array, dimension 338*> (max(min(MVAL(j),NVAL(j)))) 339*> \endverbatim 340*> 341*> \param[out] S2 342*> \verbatim 343*> S2 is DOUBLE PRECISION array, dimension 344*> (max(min(MVAL(j),NVAL(j)))) 345*> \endverbatim 346*> 347*> \param[out] X 348*> \verbatim 349*> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 350*> \endverbatim 351*> 352*> \param[in] LDX 353*> \verbatim 354*> LDX is INTEGER 355*> The leading dimension of the arrays X, Y, and Z. 356*> LDX >= max(1,MMAX) 357*> \endverbatim 358*> 359*> \param[out] Y 360*> \verbatim 361*> Y is DOUBLE PRECISION array, dimension (LDX,NRHS) 362*> \endverbatim 363*> 364*> \param[out] Z 365*> \verbatim 366*> Z is DOUBLE PRECISION array, dimension (LDX,NRHS) 367*> \endverbatim 368*> 369*> \param[out] Q 370*> \verbatim 371*> Q is DOUBLE PRECISION array, dimension (LDQ,MMAX) 372*> \endverbatim 373*> 374*> \param[in] LDQ 375*> \verbatim 376*> LDQ is INTEGER 377*> The leading dimension of the array Q. LDQ >= max(1,MMAX). 378*> \endverbatim 379*> 380*> \param[out] PT 381*> \verbatim 382*> PT is DOUBLE PRECISION array, dimension (LDPT,NMAX) 383*> \endverbatim 384*> 385*> \param[in] LDPT 386*> \verbatim 387*> LDPT is INTEGER 388*> The leading dimension of the arrays PT, U, and V. 389*> LDPT >= max(1, max(min(MVAL(j),NVAL(j)))). 390*> \endverbatim 391*> 392*> \param[out] U 393*> \verbatim 394*> U is DOUBLE PRECISION array, dimension 395*> (LDPT,max(min(MVAL(j),NVAL(j)))) 396*> \endverbatim 397*> 398*> \param[out] VT 399*> \verbatim 400*> VT is DOUBLE PRECISION array, dimension 401*> (LDPT,max(min(MVAL(j),NVAL(j)))) 402*> \endverbatim 403*> 404*> \param[out] WORK 405*> \verbatim 406*> WORK is DOUBLE PRECISION array, dimension (LWORK) 407*> \endverbatim 408*> 409*> \param[in] LWORK 410*> \verbatim 411*> LWORK is INTEGER 412*> The number of entries in WORK. This must be at least 413*> 3(M+N) and M(M + max(M,N,k) + 1) + N*min(M,N) for all 414*> pairs (M,N)=(MM(j),NN(j)) 415*> \endverbatim 416*> 417*> \param[out] IWORK 418*> \verbatim 419*> IWORK is INTEGER array, dimension at least 8*min(M,N) 420*> \endverbatim 421*> 422*> \param[in] NOUT 423*> \verbatim 424*> NOUT is INTEGER 425*> The FORTRAN unit number for printing out error messages 426*> (e.g., if a routine returns IINFO not equal to 0.) 427*> \endverbatim 428*> 429*> \param[out] INFO 430*> \verbatim 431*> INFO is INTEGER 432*> If 0, then everything ran OK. 433*> -1: NSIZES < 0 434*> -2: Some MM(j) < 0 435*> -3: Some NN(j) < 0 436*> -4: NTYPES < 0 437*> -6: NRHS < 0 438*> -8: THRESH < 0 439*> -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ). 440*> -17: LDB < 1 or LDB < MMAX. 441*> -21: LDQ < 1 or LDQ < MMAX. 442*> -23: LDPT< 1 or LDPT< MNMAX. 443*> -27: LWORK too small. 444*> If DLATMR, SLATMS, DGEBRD, DORGBR, or DBDSQR, 445*> returns an error code, the 446*> absolute value of it is returned. 447*> 448*>----------------------------------------------------------------------- 449*> 450*> Some Local Variables and Parameters: 451*> ---- ----- --------- --- ---------- 452*> 453*> ZERO, ONE Real 0 and 1. 454*> MAXTYP The number of types defined. 455*> NTEST The number of tests performed, or which can 456*> be performed so far, for the current matrix. 457*> MMAX Largest value in NN. 458*> NMAX Largest value in NN. 459*> MNMIN min(MM(j), NN(j)) (the dimension of the bidiagonal 460*> matrix.) 461*> MNMAX The maximum value of MNMIN for j=1,...,NSIZES. 462*> NFAIL The number of tests which have exceeded THRESH 463*> COND, IMODE Values to be passed to the matrix generators. 464*> ANORM Norm of A; passed to matrix generators. 465*> 466*> OVFL, UNFL Overflow and underflow thresholds. 467*> RTOVFL, RTUNFL Square roots of the previous 2 values. 468*> ULP, ULPINV Finest relative precision and its inverse. 469*> 470*> The following four arrays decode JTYPE: 471*> KTYPE(j) The general type (1-10) for type "j". 472*> KMODE(j) The MODE value to be passed to the matrix 473*> generator for type "j". 474*> KMAGN(j) The order of magnitude ( O(1), 475*> O(overflow^(1/2) ), O(underflow^(1/2) ) 476*> \endverbatim 477* 478* Authors: 479* ======== 480* 481*> \author Univ. of Tennessee 482*> \author Univ. of California Berkeley 483*> \author Univ. of Colorado Denver 484*> \author NAG Ltd. 485* 486*> \date June 2016 487* 488*> \ingroup double_eig 489* 490* ===================================================================== 491 SUBROUTINE DCHKBD( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS, 492 $ ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX, 493 $ Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK, 494 $ IWORK, NOUT, INFO ) 495* 496* -- LAPACK test routine (version 3.7.0) -- 497* -- LAPACK is a software package provided by Univ. of Tennessee, -- 498* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 499* June 2016 500* 501* .. Scalar Arguments .. 502 INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS, 503 $ NSIZES, NTYPES 504 DOUBLE PRECISION THRESH 505* .. 506* .. Array Arguments .. 507 LOGICAL DOTYPE( * ) 508 INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * ) 509 DOUBLE PRECISION A( LDA, * ), BD( * ), BE( * ), PT( LDPT, * ), 510 $ Q( LDQ, * ), S1( * ), S2( * ), U( LDPT, * ), 511 $ VT( LDPT, * ), WORK( * ), X( LDX, * ), 512 $ Y( LDX, * ), Z( LDX, * ) 513* .. 514* 515* ====================================================================== 516* 517* .. Parameters .. 518 DOUBLE PRECISION ZERO, ONE, TWO, HALF 519 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 520 $ HALF = 0.5D0 ) 521 INTEGER MAXTYP 522 PARAMETER ( MAXTYP = 16 ) 523* .. 524* .. Local Scalars .. 525 LOGICAL BADMM, BADNN, BIDIAG 526 CHARACTER UPLO 527 CHARACTER*3 PATH 528 INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, IWBD, 529 $ IWBE, IWBS, IWBZ, IWWORK, J, JCOL, JSIZE, 530 $ JTYPE, LOG2UI, M, MINWRK, MMAX, MNMAX, MNMIN, 531 $ MNMIN2, MQ, MTYPES, N, NFAIL, NMAX, 532 $ NS1, NS2, NTEST 533 DOUBLE PRECISION ABSTOL, AMNINV, ANORM, COND, OVFL, RTOVFL, 534 $ RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL, 535 $ VL, VU 536* .. 537* .. Local Arrays .. 538 INTEGER IDUM( 1 ), IOLDSD( 4 ), ISEED2( 4 ), 539 $ KMAGN( MAXTYP ), KMODE( MAXTYP ), 540 $ KTYPE( MAXTYP ) 541 DOUBLE PRECISION DUM( 1 ), DUMMA( 1 ), RESULT( 40 ) 542* .. 543* .. External Functions .. 544 DOUBLE PRECISION DLAMCH, DLARND, DSXT1 545 EXTERNAL DLAMCH, DLARND, DSXT1 546* .. 547* .. External Subroutines .. 548 EXTERNAL ALASUM, DBDSDC, DBDSQR, DBDSVDX, DBDT01, 549 $ DBDT02, DBDT03, DBDT04, DCOPY, DGEBRD, 550 $ DGEMM, DLABAD, DLACPY, DLAHD2, DLASET, 551 $ DLATMR, DLATMS, DORGBR, DORT01, XERBLA 552* .. 553* .. Intrinsic Functions .. 554 INTRINSIC ABS, EXP, INT, LOG, MAX, MIN, SQRT 555* .. 556* .. Scalars in Common .. 557 LOGICAL LERR, OK 558 CHARACTER*32 SRNAMT 559 INTEGER INFOT, NUNIT 560* .. 561* .. Common blocks .. 562 COMMON / INFOC / INFOT, NUNIT, OK, LERR 563 COMMON / SRNAMC / SRNAMT 564* .. 565* .. Data statements .. 566 DATA KTYPE / 1, 2, 5*4, 5*6, 3*9, 10 / 567 DATA KMAGN / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 / 568 DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0, 569 $ 0, 0, 0 / 570* .. 571* .. Executable Statements .. 572* 573* Check for errors 574* 575 INFO = 0 576* 577 BADMM = .FALSE. 578 BADNN = .FALSE. 579 MMAX = 1 580 NMAX = 1 581 MNMAX = 1 582 MINWRK = 1 583 DO 10 J = 1, NSIZES 584 MMAX = MAX( MMAX, MVAL( J ) ) 585 IF( MVAL( J ).LT.0 ) 586 $ BADMM = .TRUE. 587 NMAX = MAX( NMAX, NVAL( J ) ) 588 IF( NVAL( J ).LT.0 ) 589 $ BADNN = .TRUE. 590 MNMAX = MAX( MNMAX, MIN( MVAL( J ), NVAL( J ) ) ) 591 MINWRK = MAX( MINWRK, 3*( MVAL( J )+NVAL( J ) ), 592 $ MVAL( J )*( MVAL( J )+MAX( MVAL( J ), NVAL( J ), 593 $ NRHS )+1 )+NVAL( J )*MIN( NVAL( J ), MVAL( J ) ) ) 594 10 CONTINUE 595* 596* Check for errors 597* 598 IF( NSIZES.LT.0 ) THEN 599 INFO = -1 600 ELSE IF( BADMM ) THEN 601 INFO = -2 602 ELSE IF( BADNN ) THEN 603 INFO = -3 604 ELSE IF( NTYPES.LT.0 ) THEN 605 INFO = -4 606 ELSE IF( NRHS.LT.0 ) THEN 607 INFO = -6 608 ELSE IF( LDA.LT.MMAX ) THEN 609 INFO = -11 610 ELSE IF( LDX.LT.MMAX ) THEN 611 INFO = -17 612 ELSE IF( LDQ.LT.MMAX ) THEN 613 INFO = -21 614 ELSE IF( LDPT.LT.MNMAX ) THEN 615 INFO = -23 616 ELSE IF( MINWRK.GT.LWORK ) THEN 617 INFO = -27 618 END IF 619* 620 IF( INFO.NE.0 ) THEN 621 CALL XERBLA( 'DCHKBD', -INFO ) 622 RETURN 623 END IF 624* 625* Initialize constants 626* 627 PATH( 1: 1 ) = 'Double precision' 628 PATH( 2: 3 ) = 'BD' 629 NFAIL = 0 630 NTEST = 0 631 UNFL = DLAMCH( 'Safe minimum' ) 632 OVFL = DLAMCH( 'Overflow' ) 633 CALL DLABAD( UNFL, OVFL ) 634 ULP = DLAMCH( 'Precision' ) 635 ULPINV = ONE / ULP 636 LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) ) 637 RTUNFL = SQRT( UNFL ) 638 RTOVFL = SQRT( OVFL ) 639 INFOT = 0 640 ABSTOL = 2*UNFL 641* 642* Loop over sizes, types 643* 644 DO 300 JSIZE = 1, NSIZES 645 M = MVAL( JSIZE ) 646 N = NVAL( JSIZE ) 647 MNMIN = MIN( M, N ) 648 AMNINV = ONE / MAX( M, N, 1 ) 649* 650 IF( NSIZES.NE.1 ) THEN 651 MTYPES = MIN( MAXTYP, NTYPES ) 652 ELSE 653 MTYPES = MIN( MAXTYP+1, NTYPES ) 654 END IF 655* 656 DO 290 JTYPE = 1, MTYPES 657 IF( .NOT.DOTYPE( JTYPE ) ) 658 $ GO TO 290 659* 660 DO 20 J = 1, 4 661 IOLDSD( J ) = ISEED( J ) 662 20 CONTINUE 663* 664 DO 30 J = 1, 34 665 RESULT( J ) = -ONE 666 30 CONTINUE 667* 668 UPLO = ' ' 669* 670* Compute "A" 671* 672* Control parameters: 673* 674* KMAGN KMODE KTYPE 675* =1 O(1) clustered 1 zero 676* =2 large clustered 2 identity 677* =3 small exponential (none) 678* =4 arithmetic diagonal, (w/ eigenvalues) 679* =5 random symmetric, w/ eigenvalues 680* =6 nonsymmetric, w/ singular values 681* =7 random diagonal 682* =8 random symmetric 683* =9 random nonsymmetric 684* =10 random bidiagonal (log. distrib.) 685* 686 IF( MTYPES.GT.MAXTYP ) 687 $ GO TO 100 688* 689 ITYPE = KTYPE( JTYPE ) 690 IMODE = KMODE( JTYPE ) 691* 692* Compute norm 693* 694 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 695* 696 40 CONTINUE 697 ANORM = ONE 698 GO TO 70 699* 700 50 CONTINUE 701 ANORM = ( RTOVFL*ULP )*AMNINV 702 GO TO 70 703* 704 60 CONTINUE 705 ANORM = RTUNFL*MAX( M, N )*ULPINV 706 GO TO 70 707* 708 70 CONTINUE 709* 710 CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 711 IINFO = 0 712 COND = ULPINV 713* 714 BIDIAG = .FALSE. 715 IF( ITYPE.EQ.1 ) THEN 716* 717* Zero matrix 718* 719 IINFO = 0 720* 721 ELSE IF( ITYPE.EQ.2 ) THEN 722* 723* Identity 724* 725 DO 80 JCOL = 1, MNMIN 726 A( JCOL, JCOL ) = ANORM 727 80 CONTINUE 728* 729 ELSE IF( ITYPE.EQ.4 ) THEN 730* 731* Diagonal Matrix, [Eigen]values Specified 732* 733 CALL DLATMS( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, IMODE, 734 $ COND, ANORM, 0, 0, 'N', A, LDA, 735 $ WORK( MNMIN+1 ), IINFO ) 736* 737 ELSE IF( ITYPE.EQ.5 ) THEN 738* 739* Symmetric, eigenvalues specified 740* 741 CALL DLATMS( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, IMODE, 742 $ COND, ANORM, M, N, 'N', A, LDA, 743 $ WORK( MNMIN+1 ), IINFO ) 744* 745 ELSE IF( ITYPE.EQ.6 ) THEN 746* 747* Nonsymmetric, singular values specified 748* 749 CALL DLATMS( M, N, 'S', ISEED, 'N', WORK, IMODE, COND, 750 $ ANORM, M, N, 'N', A, LDA, WORK( MNMIN+1 ), 751 $ IINFO ) 752* 753 ELSE IF( ITYPE.EQ.7 ) THEN 754* 755* Diagonal, random entries 756* 757 CALL DLATMR( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, 6, ONE, 758 $ ONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE, 759 $ WORK( 2*MNMIN+1 ), 1, ONE, 'N', IWORK, 0, 0, 760 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 761* 762 ELSE IF( ITYPE.EQ.8 ) THEN 763* 764* Symmetric, random entries 765* 766 CALL DLATMR( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, 6, ONE, 767 $ ONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE, 768 $ WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N, 769 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 770* 771 ELSE IF( ITYPE.EQ.9 ) THEN 772* 773* Nonsymmetric, random entries 774* 775 CALL DLATMR( M, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 776 $ 'T', 'N', WORK( MNMIN+1 ), 1, ONE, 777 $ WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N, 778 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 779* 780 ELSE IF( ITYPE.EQ.10 ) THEN 781* 782* Bidiagonal, random entries 783* 784 TEMP1 = -TWO*LOG( ULP ) 785 DO 90 J = 1, MNMIN 786 BD( J ) = EXP( TEMP1*DLARND( 2, ISEED ) ) 787 IF( J.LT.MNMIN ) 788 $ BE( J ) = EXP( TEMP1*DLARND( 2, ISEED ) ) 789 90 CONTINUE 790* 791 IINFO = 0 792 BIDIAG = .TRUE. 793 IF( M.GE.N ) THEN 794 UPLO = 'U' 795 ELSE 796 UPLO = 'L' 797 END IF 798 ELSE 799 IINFO = 1 800 END IF 801* 802 IF( IINFO.EQ.0 ) THEN 803* 804* Generate Right-Hand Side 805* 806 IF( BIDIAG ) THEN 807 CALL DLATMR( MNMIN, NRHS, 'S', ISEED, 'N', WORK, 6, 808 $ ONE, ONE, 'T', 'N', WORK( MNMIN+1 ), 1, 809 $ ONE, WORK( 2*MNMIN+1 ), 1, ONE, 'N', 810 $ IWORK, MNMIN, NRHS, ZERO, ONE, 'NO', Y, 811 $ LDX, IWORK, IINFO ) 812 ELSE 813 CALL DLATMR( M, NRHS, 'S', ISEED, 'N', WORK, 6, ONE, 814 $ ONE, 'T', 'N', WORK( M+1 ), 1, ONE, 815 $ WORK( 2*M+1 ), 1, ONE, 'N', IWORK, M, 816 $ NRHS, ZERO, ONE, 'NO', X, LDX, IWORK, 817 $ IINFO ) 818 END IF 819 END IF 820* 821* Error Exit 822* 823 IF( IINFO.NE.0 ) THEN 824 WRITE( NOUT, FMT = 9998 )'Generator', IINFO, M, N, 825 $ JTYPE, IOLDSD 826 INFO = ABS( IINFO ) 827 RETURN 828 END IF 829* 830 100 CONTINUE 831* 832* Call DGEBRD and DORGBR to compute B, Q, and P, do tests. 833* 834 IF( .NOT.BIDIAG ) THEN 835* 836* Compute transformations to reduce A to bidiagonal form: 837* B := Q' * A * P. 838* 839 CALL DLACPY( ' ', M, N, A, LDA, Q, LDQ ) 840 CALL DGEBRD( M, N, Q, LDQ, BD, BE, WORK, WORK( MNMIN+1 ), 841 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO ) 842* 843* Check error code from DGEBRD. 844* 845 IF( IINFO.NE.0 ) THEN 846 WRITE( NOUT, FMT = 9998 )'DGEBRD', IINFO, M, N, 847 $ JTYPE, IOLDSD 848 INFO = ABS( IINFO ) 849 RETURN 850 END IF 851* 852 CALL DLACPY( ' ', M, N, Q, LDQ, PT, LDPT ) 853 IF( M.GE.N ) THEN 854 UPLO = 'U' 855 ELSE 856 UPLO = 'L' 857 END IF 858* 859* Generate Q 860* 861 MQ = M 862 IF( NRHS.LE.0 ) 863 $ MQ = MNMIN 864 CALL DORGBR( 'Q', M, MQ, N, Q, LDQ, WORK, 865 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO ) 866* 867* Check error code from DORGBR. 868* 869 IF( IINFO.NE.0 ) THEN 870 WRITE( NOUT, FMT = 9998 )'DORGBR(Q)', IINFO, M, N, 871 $ JTYPE, IOLDSD 872 INFO = ABS( IINFO ) 873 RETURN 874 END IF 875* 876* Generate P' 877* 878 CALL DORGBR( 'P', MNMIN, N, M, PT, LDPT, WORK( MNMIN+1 ), 879 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO ) 880* 881* Check error code from DORGBR. 882* 883 IF( IINFO.NE.0 ) THEN 884 WRITE( NOUT, FMT = 9998 )'DORGBR(P)', IINFO, M, N, 885 $ JTYPE, IOLDSD 886 INFO = ABS( IINFO ) 887 RETURN 888 END IF 889* 890* Apply Q' to an M by NRHS matrix X: Y := Q' * X. 891* 892 CALL DGEMM( 'Transpose', 'No transpose', M, NRHS, M, ONE, 893 $ Q, LDQ, X, LDX, ZERO, Y, LDX ) 894* 895* Test 1: Check the decomposition A := Q * B * PT 896* 2: Check the orthogonality of Q 897* 3: Check the orthogonality of PT 898* 899 CALL DBDT01( M, N, 1, A, LDA, Q, LDQ, BD, BE, PT, LDPT, 900 $ WORK, RESULT( 1 ) ) 901 CALL DORT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK, 902 $ RESULT( 2 ) ) 903 CALL DORT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK, 904 $ RESULT( 3 ) ) 905 END IF 906* 907* Use DBDSQR to form the SVD of the bidiagonal matrix B: 908* B := U * S1 * VT, and compute Z = U' * Y. 909* 910 CALL DCOPY( MNMIN, BD, 1, S1, 1 ) 911 IF( MNMIN.GT.0 ) 912 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 913 CALL DLACPY( ' ', M, NRHS, Y, LDX, Z, LDX ) 914 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, U, LDPT ) 915 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, VT, LDPT ) 916* 917 CALL DBDSQR( UPLO, MNMIN, MNMIN, MNMIN, NRHS, S1, WORK, VT, 918 $ LDPT, U, LDPT, Z, LDX, WORK( MNMIN+1 ), IINFO ) 919* 920* Check error code from DBDSQR. 921* 922 IF( IINFO.NE.0 ) THEN 923 WRITE( NOUT, FMT = 9998 )'DBDSQR(vects)', IINFO, M, N, 924 $ JTYPE, IOLDSD 925 INFO = ABS( IINFO ) 926 IF( IINFO.LT.0 ) THEN 927 RETURN 928 ELSE 929 RESULT( 4 ) = ULPINV 930 GO TO 270 931 END IF 932 END IF 933* 934* Use DBDSQR to compute only the singular values of the 935* bidiagonal matrix B; U, VT, and Z should not be modified. 936* 937 CALL DCOPY( MNMIN, BD, 1, S2, 1 ) 938 IF( MNMIN.GT.0 ) 939 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 940* 941 CALL DBDSQR( UPLO, MNMIN, 0, 0, 0, S2, WORK, VT, LDPT, U, 942 $ LDPT, Z, LDX, WORK( MNMIN+1 ), IINFO ) 943* 944* Check error code from DBDSQR. 945* 946 IF( IINFO.NE.0 ) THEN 947 WRITE( NOUT, FMT = 9998 )'DBDSQR(values)', IINFO, M, N, 948 $ JTYPE, IOLDSD 949 INFO = ABS( IINFO ) 950 IF( IINFO.LT.0 ) THEN 951 RETURN 952 ELSE 953 RESULT( 9 ) = ULPINV 954 GO TO 270 955 END IF 956 END IF 957* 958* Test 4: Check the decomposition B := U * S1 * VT 959* 5: Check the computation Z := U' * Y 960* 6: Check the orthogonality of U 961* 7: Check the orthogonality of VT 962* 963 CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT, 964 $ WORK, RESULT( 4 ) ) 965 CALL DBDT02( MNMIN, NRHS, Y, LDX, Z, LDX, U, LDPT, WORK, 966 $ RESULT( 5 ) ) 967 CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK, 968 $ RESULT( 6 ) ) 969 CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK, 970 $ RESULT( 7 ) ) 971* 972* Test 8: Check that the singular values are sorted in 973* non-increasing order and are non-negative 974* 975 RESULT( 8 ) = ZERO 976 DO 110 I = 1, MNMIN - 1 977 IF( S1( I ).LT.S1( I+1 ) ) 978 $ RESULT( 8 ) = ULPINV 979 IF( S1( I ).LT.ZERO ) 980 $ RESULT( 8 ) = ULPINV 981 110 CONTINUE 982 IF( MNMIN.GE.1 ) THEN 983 IF( S1( MNMIN ).LT.ZERO ) 984 $ RESULT( 8 ) = ULPINV 985 END IF 986* 987* Test 9: Compare DBDSQR with and without singular vectors 988* 989 TEMP2 = ZERO 990* 991 DO 120 J = 1, MNMIN 992 TEMP1 = ABS( S1( J )-S2( J ) ) / 993 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ), 994 $ ULP*MAX( ABS( S1( J ) ), ABS( S2( J ) ) ) ) 995 TEMP2 = MAX( TEMP1, TEMP2 ) 996 120 CONTINUE 997* 998 RESULT( 9 ) = TEMP2 999* 1000* Test 10: Sturm sequence test of singular values 1001* Go up by factors of two until it succeeds 1002* 1003 TEMP1 = THRESH*( HALF-ULP ) 1004* 1005 DO 130 J = 0, LOG2UI 1006* CALL DSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO ) 1007 IF( IINFO.EQ.0 ) 1008 $ GO TO 140 1009 TEMP1 = TEMP1*TWO 1010 130 CONTINUE 1011* 1012 140 CONTINUE 1013 RESULT( 10 ) = TEMP1 1014* 1015* Use DBDSQR to form the decomposition A := (QU) S (VT PT) 1016* from the bidiagonal form A := Q B PT. 1017* 1018 IF( .NOT.BIDIAG ) THEN 1019 CALL DCOPY( MNMIN, BD, 1, S2, 1 ) 1020 IF( MNMIN.GT.0 ) 1021 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 1022* 1023 CALL DBDSQR( UPLO, MNMIN, N, M, NRHS, S2, WORK, PT, LDPT, 1024 $ Q, LDQ, Y, LDX, WORK( MNMIN+1 ), IINFO ) 1025* 1026* Test 11: Check the decomposition A := Q*U * S2 * VT*PT 1027* 12: Check the computation Z := U' * Q' * X 1028* 13: Check the orthogonality of Q*U 1029* 14: Check the orthogonality of VT*PT 1030* 1031 CALL DBDT01( M, N, 0, A, LDA, Q, LDQ, S2, DUMMA, PT, 1032 $ LDPT, WORK, RESULT( 11 ) ) 1033 CALL DBDT02( M, NRHS, X, LDX, Y, LDX, Q, LDQ, WORK, 1034 $ RESULT( 12 ) ) 1035 CALL DORT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK, 1036 $ RESULT( 13 ) ) 1037 CALL DORT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK, 1038 $ RESULT( 14 ) ) 1039 END IF 1040* 1041* Use DBDSDC to form the SVD of the bidiagonal matrix B: 1042* B := U * S1 * VT 1043* 1044 CALL DCOPY( MNMIN, BD, 1, S1, 1 ) 1045 IF( MNMIN.GT.0 ) 1046 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 1047 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, U, LDPT ) 1048 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, VT, LDPT ) 1049* 1050 CALL DBDSDC( UPLO, 'I', MNMIN, S1, WORK, U, LDPT, VT, LDPT, 1051 $ DUM, IDUM, WORK( MNMIN+1 ), IWORK, IINFO ) 1052* 1053* Check error code from DBDSDC. 1054* 1055 IF( IINFO.NE.0 ) THEN 1056 WRITE( NOUT, FMT = 9998 )'DBDSDC(vects)', IINFO, M, N, 1057 $ JTYPE, IOLDSD 1058 INFO = ABS( IINFO ) 1059 IF( IINFO.LT.0 ) THEN 1060 RETURN 1061 ELSE 1062 RESULT( 15 ) = ULPINV 1063 GO TO 270 1064 END IF 1065 END IF 1066* 1067* Use DBDSDC to compute only the singular values of the 1068* bidiagonal matrix B; U and VT should not be modified. 1069* 1070 CALL DCOPY( MNMIN, BD, 1, S2, 1 ) 1071 IF( MNMIN.GT.0 ) 1072 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 1073* 1074 CALL DBDSDC( UPLO, 'N', MNMIN, S2, WORK, DUM, 1, DUM, 1, 1075 $ DUM, IDUM, WORK( MNMIN+1 ), IWORK, IINFO ) 1076* 1077* Check error code from DBDSDC. 1078* 1079 IF( IINFO.NE.0 ) THEN 1080 WRITE( NOUT, FMT = 9998 )'DBDSDC(values)', IINFO, M, N, 1081 $ JTYPE, IOLDSD 1082 INFO = ABS( IINFO ) 1083 IF( IINFO.LT.0 ) THEN 1084 RETURN 1085 ELSE 1086 RESULT( 18 ) = ULPINV 1087 GO TO 270 1088 END IF 1089 END IF 1090* 1091* Test 15: Check the decomposition B := U * S1 * VT 1092* 16: Check the orthogonality of U 1093* 17: Check the orthogonality of VT 1094* 1095 CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT, 1096 $ WORK, RESULT( 15 ) ) 1097 CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK, 1098 $ RESULT( 16 ) ) 1099 CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK, 1100 $ RESULT( 17 ) ) 1101* 1102* Test 18: Check that the singular values are sorted in 1103* non-increasing order and are non-negative 1104* 1105 RESULT( 18 ) = ZERO 1106 DO 150 I = 1, MNMIN - 1 1107 IF( S1( I ).LT.S1( I+1 ) ) 1108 $ RESULT( 18 ) = ULPINV 1109 IF( S1( I ).LT.ZERO ) 1110 $ RESULT( 18 ) = ULPINV 1111 150 CONTINUE 1112 IF( MNMIN.GE.1 ) THEN 1113 IF( S1( MNMIN ).LT.ZERO ) 1114 $ RESULT( 18 ) = ULPINV 1115 END IF 1116* 1117* Test 19: Compare DBDSQR with and without singular vectors 1118* 1119 TEMP2 = ZERO 1120* 1121 DO 160 J = 1, MNMIN 1122 TEMP1 = ABS( S1( J )-S2( J ) ) / 1123 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ), 1124 $ ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) ) 1125 TEMP2 = MAX( TEMP1, TEMP2 ) 1126 160 CONTINUE 1127* 1128 RESULT( 19 ) = TEMP2 1129* 1130* 1131* Use DBDSVDX to compute the SVD of the bidiagonal matrix B: 1132* B := U * S1 * VT 1133* 1134 IF( JTYPE.EQ.10 .OR. JTYPE.EQ.16 ) THEN 1135* ================================= 1136* Matrix types temporarily disabled 1137* ================================= 1138 RESULT( 20:34 ) = ZERO 1139 GO TO 270 1140 END IF 1141* 1142 IWBS = 1 1143 IWBD = IWBS + MNMIN 1144 IWBE = IWBD + MNMIN 1145 IWBZ = IWBE + MNMIN 1146 IWWORK = IWBZ + 2*MNMIN*(MNMIN+1) 1147 MNMIN2 = MAX( 1,MNMIN*2 ) 1148* 1149 CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 ) 1150 IF( MNMIN.GT.0 ) 1151 $ CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 ) 1152* 1153 CALL DBDSVDX( UPLO, 'V', 'A', MNMIN, WORK( IWBD ), 1154 $ WORK( IWBE ), ZERO, ZERO, 0, 0, NS1, S1, 1155 $ WORK( IWBZ ), MNMIN2, WORK( IWWORK ), 1156 $ IWORK, IINFO) 1157* 1158* Check error code from DBDSVDX. 1159* 1160 IF( IINFO.NE.0 ) THEN 1161 WRITE( NOUT, FMT = 9998 )'DBDSVDX(vects,A)', IINFO, M, N, 1162 $ JTYPE, IOLDSD 1163 INFO = ABS( IINFO ) 1164 IF( IINFO.LT.0 ) THEN 1165 RETURN 1166 ELSE 1167 RESULT( 20 ) = ULPINV 1168 GO TO 270 1169 END IF 1170 END IF 1171* 1172 J = IWBZ 1173 DO 170 I = 1, NS1 1174 CALL DCOPY( MNMIN, WORK( J ), 1, U( 1,I ), 1 ) 1175 J = J + MNMIN 1176 CALL DCOPY( MNMIN, WORK( J ), 1, VT( I,1 ), LDPT ) 1177 J = J + MNMIN 1178 170 CONTINUE 1179* 1180* Use DBDSVDX to compute only the singular values of the 1181* bidiagonal matrix B; U and VT should not be modified. 1182* 1183 IF( JTYPE.EQ.9 ) THEN 1184* ================================= 1185* Matrix types temporarily disabled 1186* ================================= 1187 RESULT( 24 ) = ZERO 1188 GO TO 270 1189 END IF 1190* 1191 CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 ) 1192 IF( MNMIN.GT.0 ) 1193 $ CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 ) 1194* 1195 CALL DBDSVDX( UPLO, 'N', 'A', MNMIN, WORK( IWBD ), 1196 $ WORK( IWBE ), ZERO, ZERO, 0, 0, NS2, S2, 1197 $ WORK( IWBZ ), MNMIN2, WORK( IWWORK ), 1198 $ IWORK, IINFO ) 1199* 1200* Check error code from DBDSVDX. 1201* 1202 IF( IINFO.NE.0 ) THEN 1203 WRITE( NOUT, FMT = 9998 )'DBDSVDX(values,A)', IINFO, 1204 $ M, N, JTYPE, IOLDSD 1205 INFO = ABS( IINFO ) 1206 IF( IINFO.LT.0 ) THEN 1207 RETURN 1208 ELSE 1209 RESULT( 24 ) = ULPINV 1210 GO TO 270 1211 END IF 1212 END IF 1213* 1214* Save S1 for tests 30-34. 1215* 1216 CALL DCOPY( MNMIN, S1, 1, WORK( IWBS ), 1 ) 1217* 1218* Test 20: Check the decomposition B := U * S1 * VT 1219* 21: Check the orthogonality of U 1220* 22: Check the orthogonality of VT 1221* 23: Check that the singular values are sorted in 1222* non-increasing order and are non-negative 1223* 24: Compare DBDSVDX with and without singular vectors 1224* 1225 CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, 1226 $ LDPT, WORK( IWBS+MNMIN ), RESULT( 20 ) ) 1227 CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, 1228 $ WORK( IWBS+MNMIN ), LWORK-MNMIN, 1229 $ RESULT( 21 ) ) 1230 CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, 1231 $ WORK( IWBS+MNMIN ), LWORK-MNMIN, 1232 $ RESULT( 22) ) 1233* 1234 RESULT( 23 ) = ZERO 1235 DO 180 I = 1, MNMIN - 1 1236 IF( S1( I ).LT.S1( I+1 ) ) 1237 $ RESULT( 23 ) = ULPINV 1238 IF( S1( I ).LT.ZERO ) 1239 $ RESULT( 23 ) = ULPINV 1240 180 CONTINUE 1241 IF( MNMIN.GE.1 ) THEN 1242 IF( S1( MNMIN ).LT.ZERO ) 1243 $ RESULT( 23 ) = ULPINV 1244 END IF 1245* 1246 TEMP2 = ZERO 1247 DO 190 J = 1, MNMIN 1248 TEMP1 = ABS( S1( J )-S2( J ) ) / 1249 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ), 1250 $ ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) ) 1251 TEMP2 = MAX( TEMP1, TEMP2 ) 1252 190 CONTINUE 1253 RESULT( 24 ) = TEMP2 1254 ANORM = S1( 1 ) 1255* 1256* Use DBDSVDX with RANGE='I': choose random values for IL and 1257* IU, and ask for the IL-th through IU-th singular values 1258* and corresponding vectors. 1259* 1260 DO 200 I = 1, 4 1261 ISEED2( I ) = ISEED( I ) 1262 200 CONTINUE 1263 IF( MNMIN.LE.1 ) THEN 1264 IL = 1 1265 IU = MNMIN 1266 ELSE 1267 IL = 1 + INT( ( MNMIN-1 )*DLARND( 1, ISEED2 ) ) 1268 IU = 1 + INT( ( MNMIN-1 )*DLARND( 1, ISEED2 ) ) 1269 IF( IU.LT.IL ) THEN 1270 ITEMP = IU 1271 IU = IL 1272 IL = ITEMP 1273 END IF 1274 END IF 1275* 1276 CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 ) 1277 IF( MNMIN.GT.0 ) 1278 $ CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 ) 1279* 1280 CALL DBDSVDX( UPLO, 'V', 'I', MNMIN, WORK( IWBD ), 1281 $ WORK( IWBE ), ZERO, ZERO, IL, IU, NS1, S1, 1282 $ WORK( IWBZ ), MNMIN2, WORK( IWWORK ), 1283 $ IWORK, IINFO) 1284* 1285* Check error code from DBDSVDX. 1286* 1287 IF( IINFO.NE.0 ) THEN 1288 WRITE( NOUT, FMT = 9998 )'DBDSVDX(vects,I)', IINFO, 1289 $ M, N, JTYPE, IOLDSD 1290 INFO = ABS( IINFO ) 1291 IF( IINFO.LT.0 ) THEN 1292 RETURN 1293 ELSE 1294 RESULT( 25 ) = ULPINV 1295 GO TO 270 1296 END IF 1297 END IF 1298* 1299 J = IWBZ 1300 DO 210 I = 1, NS1 1301 CALL DCOPY( MNMIN, WORK( J ), 1, U( 1,I ), 1 ) 1302 J = J + MNMIN 1303 CALL DCOPY( MNMIN, WORK( J ), 1, VT( I,1 ), LDPT ) 1304 J = J + MNMIN 1305 210 CONTINUE 1306* 1307* Use DBDSVDX to compute only the singular values of the 1308* bidiagonal matrix B; U and VT should not be modified. 1309* 1310 CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 ) 1311 IF( MNMIN.GT.0 ) 1312 $ CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 ) 1313* 1314 CALL DBDSVDX( UPLO, 'N', 'I', MNMIN, WORK( IWBD ), 1315 $ WORK( IWBE ), ZERO, ZERO, IL, IU, NS2, S2, 1316 $ WORK( IWBZ ), MNMIN2, WORK( IWWORK ), 1317 $ IWORK, IINFO ) 1318* 1319* Check error code from DBDSVDX. 1320* 1321 IF( IINFO.NE.0 ) THEN 1322 WRITE( NOUT, FMT = 9998 )'DBDSVDX(values,I)', IINFO, 1323 $ M, N, JTYPE, IOLDSD 1324 INFO = ABS( IINFO ) 1325 IF( IINFO.LT.0 ) THEN 1326 RETURN 1327 ELSE 1328 RESULT( 29 ) = ULPINV 1329 GO TO 270 1330 END IF 1331 END IF 1332* 1333* Test 25: Check S1 - U' * B * VT' 1334* 26: Check the orthogonality of U 1335* 27: Check the orthogonality of VT 1336* 28: Check that the singular values are sorted in 1337* non-increasing order and are non-negative 1338* 29: Compare DBDSVDX with and without singular vectors 1339* 1340 CALL DBDT04( UPLO, MNMIN, BD, BE, S1, NS1, U, 1341 $ LDPT, VT, LDPT, WORK( IWBS+MNMIN ), 1342 $ RESULT( 25 ) ) 1343 CALL DORT01( 'Columns', MNMIN, NS1, U, LDPT, 1344 $ WORK( IWBS+MNMIN ), LWORK-MNMIN, 1345 $ RESULT( 26 ) ) 1346 CALL DORT01( 'Rows', NS1, MNMIN, VT, LDPT, 1347 $ WORK( IWBS+MNMIN ), LWORK-MNMIN, 1348 $ RESULT( 27 ) ) 1349* 1350 RESULT( 28 ) = ZERO 1351 DO 220 I = 1, NS1 - 1 1352 IF( S1( I ).LT.S1( I+1 ) ) 1353 $ RESULT( 28 ) = ULPINV 1354 IF( S1( I ).LT.ZERO ) 1355 $ RESULT( 28 ) = ULPINV 1356 220 CONTINUE 1357 IF( NS1.GE.1 ) THEN 1358 IF( S1( NS1 ).LT.ZERO ) 1359 $ RESULT( 28 ) = ULPINV 1360 END IF 1361* 1362 TEMP2 = ZERO 1363 DO 230 J = 1, NS1 1364 TEMP1 = ABS( S1( J )-S2( J ) ) / 1365 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ), 1366 $ ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) ) 1367 TEMP2 = MAX( TEMP1, TEMP2 ) 1368 230 CONTINUE 1369 RESULT( 29 ) = TEMP2 1370* 1371* Use DBDSVDX with RANGE='V': determine the values VL and VU 1372* of the IL-th and IU-th singular values and ask for all 1373* singular values in this range. 1374* 1375 CALL DCOPY( MNMIN, WORK( IWBS ), 1, S1, 1 ) 1376* 1377 IF( MNMIN.GT.0 ) THEN 1378 IF( IL.NE.1 ) THEN 1379 VU = S1( IL ) + MAX( HALF*ABS( S1( IL )-S1( IL-1 ) ), 1380 $ ULP*ANORM, TWO*RTUNFL ) 1381 ELSE 1382 VU = S1( 1 ) + MAX( HALF*ABS( S1( MNMIN )-S1( 1 ) ), 1383 $ ULP*ANORM, TWO*RTUNFL ) 1384 END IF 1385 IF( IU.NE.NS1 ) THEN 1386 VL = S1( IU ) - MAX( ULP*ANORM, TWO*RTUNFL, 1387 $ HALF*ABS( S1( IU+1 )-S1( IU ) ) ) 1388 ELSE 1389 VL = S1( NS1 ) - MAX( ULP*ANORM, TWO*RTUNFL, 1390 $ HALF*ABS( S1( MNMIN )-S1( 1 ) ) ) 1391 END IF 1392 VL = MAX( VL,ZERO ) 1393 VU = MAX( VU,ZERO ) 1394 IF( VL.GE.VU ) VU = MAX( VU*2, VU+VL+HALF ) 1395 ELSE 1396 VL = ZERO 1397 VU = ONE 1398 END IF 1399* 1400 CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 ) 1401 IF( MNMIN.GT.0 ) 1402 $ CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 ) 1403* 1404 CALL DBDSVDX( UPLO, 'V', 'V', MNMIN, WORK( IWBD ), 1405 $ WORK( IWBE ), VL, VU, 0, 0, NS1, S1, 1406 $ WORK( IWBZ ), MNMIN2, WORK( IWWORK ), 1407 $ IWORK, IINFO ) 1408* 1409* Check error code from DBDSVDX. 1410* 1411 IF( IINFO.NE.0 ) THEN 1412 WRITE( NOUT, FMT = 9998 )'DBDSVDX(vects,V)', IINFO, 1413 $ M, N, JTYPE, IOLDSD 1414 INFO = ABS( IINFO ) 1415 IF( IINFO.LT.0 ) THEN 1416 RETURN 1417 ELSE 1418 RESULT( 30 ) = ULPINV 1419 GO TO 270 1420 END IF 1421 END IF 1422* 1423 J = IWBZ 1424 DO 240 I = 1, NS1 1425 CALL DCOPY( MNMIN, WORK( J ), 1, U( 1,I ), 1 ) 1426 J = J + MNMIN 1427 CALL DCOPY( MNMIN, WORK( J ), 1, VT( I,1 ), LDPT ) 1428 J = J + MNMIN 1429 240 CONTINUE 1430* 1431* Use DBDSVDX to compute only the singular values of the 1432* bidiagonal matrix B; U and VT should not be modified. 1433* 1434 CALL DCOPY( MNMIN, BD, 1, WORK( IWBD ), 1 ) 1435 IF( MNMIN.GT.0 ) 1436 $ CALL DCOPY( MNMIN-1, BE, 1, WORK( IWBE ), 1 ) 1437* 1438 CALL DBDSVDX( UPLO, 'N', 'V', MNMIN, WORK( IWBD ), 1439 $ WORK( IWBE ), VL, VU, 0, 0, NS2, S2, 1440 $ WORK( IWBZ ), MNMIN2, WORK( IWWORK ), 1441 $ IWORK, IINFO ) 1442* 1443* Check error code from DBDSVDX. 1444* 1445 IF( IINFO.NE.0 ) THEN 1446 WRITE( NOUT, FMT = 9998 )'DBDSVDX(values,V)', IINFO, 1447 $ M, N, JTYPE, IOLDSD 1448 INFO = ABS( IINFO ) 1449 IF( IINFO.LT.0 ) THEN 1450 RETURN 1451 ELSE 1452 RESULT( 34 ) = ULPINV 1453 GO TO 270 1454 END IF 1455 END IF 1456* 1457* Test 30: Check S1 - U' * B * VT' 1458* 31: Check the orthogonality of U 1459* 32: Check the orthogonality of VT 1460* 33: Check that the singular values are sorted in 1461* non-increasing order and are non-negative 1462* 34: Compare DBDSVDX with and without singular vectors 1463* 1464 CALL DBDT04( UPLO, MNMIN, BD, BE, S1, NS1, U, 1465 $ LDPT, VT, LDPT, WORK( IWBS+MNMIN ), 1466 $ RESULT( 30 ) ) 1467 CALL DORT01( 'Columns', MNMIN, NS1, U, LDPT, 1468 $ WORK( IWBS+MNMIN ), LWORK-MNMIN, 1469 $ RESULT( 31 ) ) 1470 CALL DORT01( 'Rows', NS1, MNMIN, VT, LDPT, 1471 $ WORK( IWBS+MNMIN ), LWORK-MNMIN, 1472 $ RESULT( 32 ) ) 1473* 1474 RESULT( 33 ) = ZERO 1475 DO 250 I = 1, NS1 - 1 1476 IF( S1( I ).LT.S1( I+1 ) ) 1477 $ RESULT( 28 ) = ULPINV 1478 IF( S1( I ).LT.ZERO ) 1479 $ RESULT( 28 ) = ULPINV 1480 250 CONTINUE 1481 IF( NS1.GE.1 ) THEN 1482 IF( S1( NS1 ).LT.ZERO ) 1483 $ RESULT( 28 ) = ULPINV 1484 END IF 1485* 1486 TEMP2 = ZERO 1487 DO 260 J = 1, NS1 1488 TEMP1 = ABS( S1( J )-S2( J ) ) / 1489 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ), 1490 $ ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) ) 1491 TEMP2 = MAX( TEMP1, TEMP2 ) 1492 260 CONTINUE 1493 RESULT( 34 ) = TEMP2 1494* 1495* End of Loop -- Check for RESULT(j) > THRESH 1496* 1497 270 CONTINUE 1498* 1499 DO 280 J = 1, 34 1500 IF( RESULT( J ).GE.THRESH ) THEN 1501 IF( NFAIL.EQ.0 ) 1502 $ CALL DLAHD2( NOUT, PATH ) 1503 WRITE( NOUT, FMT = 9999 )M, N, JTYPE, IOLDSD, J, 1504 $ RESULT( J ) 1505 NFAIL = NFAIL + 1 1506 END IF 1507 280 CONTINUE 1508 IF( .NOT.BIDIAG ) THEN 1509 NTEST = NTEST + 34 1510 ELSE 1511 NTEST = NTEST + 30 1512 END IF 1513* 1514 290 CONTINUE 1515 300 CONTINUE 1516* 1517* Summary 1518* 1519 CALL ALASUM( PATH, NOUT, NFAIL, NTEST, 0 ) 1520* 1521 RETURN 1522* 1523* End of DCHKBD 1524* 1525 9999 FORMAT( ' M=', I5, ', N=', I5, ', type ', I2, ', seed=', 1526 $ 4( I4, ',' ), ' test(', I2, ')=', G11.4 ) 1527 9998 FORMAT( ' DCHKBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=', 1528 $ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), 1529 $ I5, ')' ) 1530* 1531 END 1532