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