1*> \brief \b CDRVBD 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 CDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH, 12* A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S, 13* SSAV, E, WORK, LWORK, RWORK, IWORK, NOUNIT, 14* INFO ) 15* 16* .. Scalar Arguments .. 17* INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES, 18* $ NTYPES 19* REAL THRESH 20* .. 21* .. Array Arguments .. 22* LOGICAL DOTYPE( * ) 23* INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * ) 24* REAL E( * ), RWORK( * ), S( * ), SSAV( * ) 25* COMPLEX A( LDA, * ), ASAV( LDA, * ), U( LDU, * ), 26* $ USAV( LDU, * ), VT( LDVT, * ), 27* $ VTSAV( LDVT, * ), WORK( * ) 28* .. 29* 30* 31*> \par Purpose: 32* ============= 33*> 34*> \verbatim 35*> 36*> CDRVBD checks the singular value decomposition (SVD) driver CGESVD 37*> and CGESDD. 38*> CGESVD and CGESDD factors A = U diag(S) VT, where U and VT are 39*> unitary and diag(S) is diagonal with the entries of the array S on 40*> its diagonal. The entries of S are the singular values, nonnegative 41*> and stored in decreasing order. U and VT can be optionally not 42*> computed, overwritten on A, or computed partially. 43*> 44*> A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN. 45*> U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N. 46*> 47*> When CDRVBD is called, a number of matrix "sizes" (M's and N's) 48*> and a number of matrix "types" are specified. For each size (M,N) 49*> and each type of matrix, and for the minimal workspace as well as 50*> workspace adequate to permit blocking, an M x N matrix "A" will be 51*> generated and used to test the SVD routines. For each matrix, A will 52*> be factored as A = U diag(S) VT and the following 12 tests computed: 53*> 54*> Test for CGESVD: 55*> 56*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp ) 57*> 58*> (2) | I - U'U | / ( M ulp ) 59*> 60*> (3) | I - VT VT' | / ( N ulp ) 61*> 62*> (4) S contains MNMIN nonnegative values in decreasing order. 63*> (Return 0 if true, 1/ULP if false.) 64*> 65*> (5) | U - Upartial | / ( M ulp ) where Upartial is a partially 66*> computed U. 67*> 68*> (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially 69*> computed VT. 70*> 71*> (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the 72*> vector of singular values from the partial SVD 73*> 74*> Test for CGESDD: 75*> 76*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp ) 77*> 78*> (2) | I - U'U | / ( M ulp ) 79*> 80*> (3) | I - VT VT' | / ( N ulp ) 81*> 82*> (4) S contains MNMIN nonnegative values in decreasing order. 83*> (Return 0 if true, 1/ULP if false.) 84*> 85*> (5) | U - Upartial | / ( M ulp ) where Upartial is a partially 86*> computed U. 87*> 88*> (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially 89*> computed VT. 90*> 91*> (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the 92*> vector of singular values from the partial SVD 93*> 94*> The "sizes" are specified by the arrays MM(1:NSIZES) and 95*> NN(1:NSIZES); the value of each element pair (MM(j),NN(j)) 96*> specifies one size. The "types" are specified by a logical array 97*> DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j" 98*> will be generated. 99*> Currently, the list of possible types is: 100*> 101*> (1) The zero matrix. 102*> (2) The identity matrix. 103*> (3) A matrix of the form U D V, where U and V are unitary and 104*> D has evenly spaced entries 1, ..., ULP with random signs 105*> on the diagonal. 106*> (4) Same as (3), but multiplied by the underflow-threshold / ULP. 107*> (5) Same as (3), but multiplied by the overflow-threshold * ULP. 108*> \endverbatim 109* 110* Arguments: 111* ========== 112* 113*> \param[in] NSIZES 114*> \verbatim 115*> NSIZES is INTEGER 116*> The number of sizes of matrices to use. If it is zero, 117*> CDRVBD does nothing. It must be at least zero. 118*> \endverbatim 119*> 120*> \param[in] MM 121*> \verbatim 122*> MM is INTEGER array, dimension (NSIZES) 123*> An array containing the matrix "heights" to be used. For 124*> each j=1,...,NSIZES, if MM(j) is zero, then MM(j) and NN(j) 125*> will be ignored. The MM(j) values must be at least zero. 126*> \endverbatim 127*> 128*> \param[in] NN 129*> \verbatim 130*> NN is INTEGER array, dimension (NSIZES) 131*> An array containing the matrix "widths" to be used. For 132*> each j=1,...,NSIZES, if NN(j) is zero, then MM(j) and NN(j) 133*> will be ignored. The NN(j) values must be at least zero. 134*> \endverbatim 135*> 136*> \param[in] NTYPES 137*> \verbatim 138*> NTYPES is INTEGER 139*> The number of elements in DOTYPE. If it is zero, CDRVBD 140*> does nothing. It must be at least zero. If it is MAXTYP+1 141*> and NSIZES is 1, then an additional type, MAXTYP+1 is 142*> defined, which is to use whatever matrices are in A and B. 143*> This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 144*> DOTYPE(MAXTYP+1) is .TRUE. . 145*> \endverbatim 146*> 147*> \param[in] DOTYPE 148*> \verbatim 149*> DOTYPE is LOGICAL array, dimension (NTYPES) 150*> If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix 151*> of type j will be generated. If NTYPES is smaller than the 152*> maximum number of types defined (PARAMETER MAXTYP), then 153*> types NTYPES+1 through MAXTYP will not be generated. If 154*> NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through 155*> DOTYPE(NTYPES) will be ignored. 156*> \endverbatim 157*> 158*> \param[in,out] ISEED 159*> \verbatim 160*> ISEED is INTEGER array, dimension (4) 161*> On entry ISEED specifies the seed of the random number 162*> generator. The array elements should be between 0 and 4095; 163*> if not they will be reduced mod 4096. Also, ISEED(4) must 164*> be odd. The random number generator uses a linear 165*> congruential sequence limited to small integers, and so 166*> should produce machine independent random numbers. The 167*> values of ISEED are changed on exit, and can be used in the 168*> next call to CDRVBD to continue the same random number 169*> sequence. 170*> \endverbatim 171*> 172*> \param[in] THRESH 173*> \verbatim 174*> THRESH is REAL 175*> A test will count as "failed" if the "error", computed as 176*> described above, exceeds THRESH. Note that the error 177*> is scaled to be O(1), so THRESH should be a reasonably 178*> small multiple of 1, e.g., 10 or 100. In particular, 179*> it should not depend on the precision (single vs. double) 180*> or the size of the matrix. It must be at least zero. 181*> \endverbatim 182*> 183*> \param[out] A 184*> \verbatim 185*> A is COMPLEX array, dimension (LDA,max(NN)) 186*> Used to hold the matrix whose singular values are to be 187*> computed. On exit, A contains the last matrix actually 188*> used. 189*> \endverbatim 190*> 191*> \param[in] LDA 192*> \verbatim 193*> LDA is INTEGER 194*> The leading dimension of A. It must be at 195*> least 1 and at least max( MM ). 196*> \endverbatim 197*> 198*> \param[out] U 199*> \verbatim 200*> U is COMPLEX array, dimension (LDU,max(MM)) 201*> Used to hold the computed matrix of right singular vectors. 202*> On exit, U contains the last such vectors actually computed. 203*> \endverbatim 204*> 205*> \param[in] LDU 206*> \verbatim 207*> LDU is INTEGER 208*> The leading dimension of U. It must be at 209*> least 1 and at least max( MM ). 210*> \endverbatim 211*> 212*> \param[out] VT 213*> \verbatim 214*> VT is COMPLEX array, dimension (LDVT,max(NN)) 215*> Used to hold the computed matrix of left singular vectors. 216*> On exit, VT contains the last such vectors actually computed. 217*> \endverbatim 218*> 219*> \param[in] LDVT 220*> \verbatim 221*> LDVT is INTEGER 222*> The leading dimension of VT. It must be at 223*> least 1 and at least max( NN ). 224*> \endverbatim 225*> 226*> \param[out] ASAV 227*> \verbatim 228*> ASAV is COMPLEX array, dimension (LDA,max(NN)) 229*> Used to hold a different copy of the matrix whose singular 230*> values are to be computed. On exit, A contains the last 231*> matrix actually used. 232*> \endverbatim 233*> 234*> \param[out] USAV 235*> \verbatim 236*> USAV is COMPLEX array, dimension (LDU,max(MM)) 237*> Used to hold a different copy of the computed matrix of 238*> right singular vectors. On exit, USAV contains the last such 239*> vectors actually computed. 240*> \endverbatim 241*> 242*> \param[out] VTSAV 243*> \verbatim 244*> VTSAV is COMPLEX array, dimension (LDVT,max(NN)) 245*> Used to hold a different copy of the computed matrix of 246*> left singular vectors. On exit, VTSAV contains the last such 247*> vectors actually computed. 248*> \endverbatim 249*> 250*> \param[out] S 251*> \verbatim 252*> S is REAL array, dimension (max(min(MM,NN))) 253*> Contains the computed singular values. 254*> \endverbatim 255*> 256*> \param[out] SSAV 257*> \verbatim 258*> SSAV is REAL array, dimension (max(min(MM,NN))) 259*> Contains another copy of the computed singular values. 260*> \endverbatim 261*> 262*> \param[out] E 263*> \verbatim 264*> E is REAL array, dimension (max(min(MM,NN))) 265*> Workspace for CGESVD. 266*> \endverbatim 267*> 268*> \param[out] WORK 269*> \verbatim 270*> WORK is COMPLEX array, dimension (LWORK) 271*> \endverbatim 272*> 273*> \param[in] LWORK 274*> \verbatim 275*> LWORK is INTEGER 276*> The number of entries in WORK. This must be at least 277*> MAX(3*MIN(M,N)+MAX(M,N)**2,5*MIN(M,N),3*MAX(M,N)) for all 278*> pairs (M,N)=(MM(j),NN(j)) 279*> \endverbatim 280*> 281*> \param[out] RWORK 282*> \verbatim 283*> RWORK is REAL array, 284*> dimension ( 5*max(max(MM,NN)) ) 285*> \endverbatim 286*> 287*> \param[out] IWORK 288*> \verbatim 289*> IWORK is INTEGER array, dimension at least 8*min(M,N) 290*> \endverbatim 291*> 292*> \param[in] NOUNIT 293*> \verbatim 294*> NOUNIT is INTEGER 295*> The FORTRAN unit number for printing out error messages 296*> (e.g., if a routine returns IINFO not equal to 0.) 297*> \endverbatim 298*> 299*> \param[out] INFO 300*> \verbatim 301*> INFO is INTEGER 302*> If 0, then everything ran OK. 303*> -1: NSIZES < 0 304*> -2: Some MM(j) < 0 305*> -3: Some NN(j) < 0 306*> -4: NTYPES < 0 307*> -7: THRESH < 0 308*> -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ). 309*> -12: LDU < 1 or LDU < MMAX. 310*> -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ). 311*> -21: LWORK too small. 312*> If CLATMS, or CGESVD returns an error code, the 313*> absolute value of it is returned. 314*> \endverbatim 315* 316* Authors: 317* ======== 318* 319*> \author Univ. of Tennessee 320*> \author Univ. of California Berkeley 321*> \author Univ. of Colorado Denver 322*> \author NAG Ltd. 323* 324*> \date November 2011 325* 326*> \ingroup complex_eig 327* 328* ===================================================================== 329 SUBROUTINE CDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH, 330 $ A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S, 331 $ SSAV, E, WORK, LWORK, RWORK, IWORK, NOUNIT, 332 $ INFO ) 333* 334* -- LAPACK test routine (version 3.4.0) -- 335* -- LAPACK is a software package provided by Univ. of Tennessee, -- 336* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 337* November 2011 338* 339* .. Scalar Arguments .. 340 INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES, 341 $ NTYPES 342 REAL THRESH 343* .. 344* .. Array Arguments .. 345 LOGICAL DOTYPE( * ) 346 INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * ) 347 REAL E( * ), RWORK( * ), S( * ), SSAV( * ) 348 COMPLEX A( LDA, * ), ASAV( LDA, * ), U( LDU, * ), 349 $ USAV( LDU, * ), VT( LDVT, * ), 350 $ VTSAV( LDVT, * ), WORK( * ) 351* .. 352* 353* ===================================================================== 354* 355* .. Parameters .. 356 REAL ZERO, ONE 357 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 358 COMPLEX CZERO, CONE 359 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 360 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 361 INTEGER MAXTYP 362 PARAMETER ( MAXTYP = 5 ) 363* .. 364* .. Local Scalars .. 365 LOGICAL BADMM, BADNN 366 CHARACTER JOBQ, JOBU, JOBVT 367 INTEGER I, IINFO, IJQ, IJU, IJVT, IWSPC, IWTMP, J, 368 $ JSIZE, JTYPE, LSWORK, M, MINWRK, MMAX, MNMAX, 369 $ MNMIN, MTYPES, N, NERRS, NFAIL, NMAX, NTEST, 370 $ NTESTF, NTESTT 371 REAL ANORM, DIF, DIV, OVFL, ULP, ULPINV, UNFL 372* .. 373* .. Local Arrays .. 374 CHARACTER CJOB( 4 ) 375 INTEGER IOLDSD( 4 ) 376 REAL RESULT( 14 ) 377* .. 378* .. External Functions .. 379 REAL SLAMCH 380 EXTERNAL SLAMCH 381* .. 382* .. External Subroutines .. 383 EXTERNAL ALASVM, CBDT01, CGESDD, CGESVD, CLACPY, CLASET, 384 $ CLATMS, CUNT01, CUNT03, XERBLA 385* .. 386* .. Intrinsic Functions .. 387 INTRINSIC ABS, MAX, MIN, REAL 388* .. 389* .. Data statements .. 390 DATA CJOB / 'N', 'O', 'S', 'A' / 391* .. 392* .. Executable Statements .. 393* 394* Check for errors 395* 396 INFO = 0 397* 398* Important constants 399* 400 NERRS = 0 401 NTESTT = 0 402 NTESTF = 0 403 BADMM = .FALSE. 404 BADNN = .FALSE. 405 MMAX = 1 406 NMAX = 1 407 MNMAX = 1 408 MINWRK = 1 409 DO 10 J = 1, NSIZES 410 MMAX = MAX( MMAX, MM( J ) ) 411 IF( MM( J ).LT.0 ) 412 $ BADMM = .TRUE. 413 NMAX = MAX( NMAX, NN( J ) ) 414 IF( NN( J ).LT.0 ) 415 $ BADNN = .TRUE. 416 MNMAX = MAX( MNMAX, MIN( MM( J ), NN( J ) ) ) 417 MINWRK = MAX( MINWRK, MAX( 3*MIN( MM( J ), 418 $ NN( J ) )+MAX( MM( J ), NN( J ) )**2, 5*MIN( MM( J ), 419 $ NN( J ) ), 3*MAX( MM( J ), NN( J ) ) ) ) 420 10 CONTINUE 421* 422* Check for errors 423* 424 IF( NSIZES.LT.0 ) THEN 425 INFO = -1 426 ELSE IF( BADMM ) THEN 427 INFO = -2 428 ELSE IF( BADNN ) THEN 429 INFO = -3 430 ELSE IF( NTYPES.LT.0 ) THEN 431 INFO = -4 432 ELSE IF( LDA.LT.MAX( 1, MMAX ) ) THEN 433 INFO = -10 434 ELSE IF( LDU.LT.MAX( 1, MMAX ) ) THEN 435 INFO = -12 436 ELSE IF( LDVT.LT.MAX( 1, NMAX ) ) THEN 437 INFO = -14 438 ELSE IF( MINWRK.GT.LWORK ) THEN 439 INFO = -21 440 END IF 441* 442 IF( INFO.NE.0 ) THEN 443 CALL XERBLA( 'CDRVBD', -INFO ) 444 RETURN 445 END IF 446* 447* Quick return if nothing to do 448* 449 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 450 $ RETURN 451* 452* More Important constants 453* 454 UNFL = SLAMCH( 'S' ) 455 OVFL = ONE / UNFL 456 ULP = SLAMCH( 'E' ) 457 ULPINV = ONE / ULP 458* 459* Loop over sizes, types 460* 461 NERRS = 0 462* 463 DO 180 JSIZE = 1, NSIZES 464 M = MM( JSIZE ) 465 N = NN( JSIZE ) 466 MNMIN = MIN( M, N ) 467* 468 IF( NSIZES.NE.1 ) THEN 469 MTYPES = MIN( MAXTYP, NTYPES ) 470 ELSE 471 MTYPES = MIN( MAXTYP+1, NTYPES ) 472 END IF 473* 474 DO 170 JTYPE = 1, MTYPES 475 IF( .NOT.DOTYPE( JTYPE ) ) 476 $ GO TO 170 477 NTEST = 0 478* 479 DO 20 J = 1, 4 480 IOLDSD( J ) = ISEED( J ) 481 20 CONTINUE 482* 483* Compute "A" 484* 485 IF( MTYPES.GT.MAXTYP ) 486 $ GO TO 50 487* 488 IF( JTYPE.EQ.1 ) THEN 489* 490* Zero matrix 491* 492 CALL CLASET( 'Full', M, N, CZERO, CZERO, A, LDA ) 493 DO 30 I = 1, MIN( M, N ) 494 S( I ) = ZERO 495 30 CONTINUE 496* 497 ELSE IF( JTYPE.EQ.2 ) THEN 498* 499* Identity matrix 500* 501 CALL CLASET( 'Full', M, N, CZERO, CONE, A, LDA ) 502 DO 40 I = 1, MIN( M, N ) 503 S( I ) = ONE 504 40 CONTINUE 505* 506 ELSE 507* 508* (Scaled) random matrix 509* 510 IF( JTYPE.EQ.3 ) 511 $ ANORM = ONE 512 IF( JTYPE.EQ.4 ) 513 $ ANORM = UNFL / ULP 514 IF( JTYPE.EQ.5 ) 515 $ ANORM = OVFL*ULP 516 CALL CLATMS( M, N, 'U', ISEED, 'N', S, 4, REAL( MNMIN ), 517 $ ANORM, M-1, N-1, 'N', A, LDA, WORK, IINFO ) 518 IF( IINFO.NE.0 ) THEN 519 WRITE( NOUNIT, FMT = 9996 )'Generator', IINFO, M, N, 520 $ JTYPE, IOLDSD 521 INFO = ABS( IINFO ) 522 RETURN 523 END IF 524 END IF 525* 526 50 CONTINUE 527 CALL CLACPY( 'F', M, N, A, LDA, ASAV, LDA ) 528* 529* Do for minimal and adequate (for blocking) workspace 530* 531 DO 160 IWSPC = 1, 4 532* 533* Test for CGESVD 534* 535 IWTMP = 2*MIN( M, N )+MAX( M, N ) 536 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3 537 LSWORK = MIN( LSWORK, LWORK ) 538 LSWORK = MAX( LSWORK, 1 ) 539 IF( IWSPC.EQ.4 ) 540 $ LSWORK = LWORK 541* 542 DO 60 J = 1, 14 543 RESULT( J ) = -ONE 544 60 CONTINUE 545* 546* Factorize A 547* 548 IF( IWSPC.GT.1 ) 549 $ CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA ) 550 CALL CGESVD( 'A', 'A', M, N, A, LDA, SSAV, USAV, LDU, 551 $ VTSAV, LDVT, WORK, LSWORK, RWORK, IINFO ) 552 IF( IINFO.NE.0 ) THEN 553 WRITE( NOUNIT, FMT = 9995 )'GESVD', IINFO, M, N, 554 $ JTYPE, LSWORK, IOLDSD 555 INFO = ABS( IINFO ) 556 RETURN 557 END IF 558* 559* Do tests 1--4 560* 561 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E, 562 $ VTSAV, LDVT, WORK, RWORK, RESULT( 1 ) ) 563 IF( M.NE.0 .AND. N.NE.0 ) THEN 564 CALL CUNT01( 'Columns', MNMIN, M, USAV, LDU, WORK, 565 $ LWORK, RWORK, RESULT( 2 ) ) 566 CALL CUNT01( 'Rows', MNMIN, N, VTSAV, LDVT, WORK, 567 $ LWORK, RWORK, RESULT( 3 ) ) 568 END IF 569 RESULT( 4 ) = 0 570 DO 70 I = 1, MNMIN - 1 571 IF( SSAV( I ).LT.SSAV( I+1 ) ) 572 $ RESULT( 4 ) = ULPINV 573 IF( SSAV( I ).LT.ZERO ) 574 $ RESULT( 4 ) = ULPINV 575 70 CONTINUE 576 IF( MNMIN.GE.1 ) THEN 577 IF( SSAV( MNMIN ).LT.ZERO ) 578 $ RESULT( 4 ) = ULPINV 579 END IF 580* 581* Do partial SVDs, comparing to SSAV, USAV, and VTSAV 582* 583 RESULT( 5 ) = ZERO 584 RESULT( 6 ) = ZERO 585 RESULT( 7 ) = ZERO 586 DO 100 IJU = 0, 3 587 DO 90 IJVT = 0, 3 588 IF( ( IJU.EQ.3 .AND. IJVT.EQ.3 ) .OR. 589 $ ( IJU.EQ.1 .AND. IJVT.EQ.1 ) )GO TO 90 590 JOBU = CJOB( IJU+1 ) 591 JOBVT = CJOB( IJVT+1 ) 592 CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA ) 593 CALL CGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, 594 $ VT, LDVT, WORK, LSWORK, RWORK, IINFO ) 595* 596* Compare U 597* 598 DIF = ZERO 599 IF( M.GT.0 .AND. N.GT.0 ) THEN 600 IF( IJU.EQ.1 ) THEN 601 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV, 602 $ LDU, A, LDA, WORK, LWORK, RWORK, 603 $ DIF, IINFO ) 604 ELSE IF( IJU.EQ.2 ) THEN 605 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV, 606 $ LDU, U, LDU, WORK, LWORK, RWORK, 607 $ DIF, IINFO ) 608 ELSE IF( IJU.EQ.3 ) THEN 609 CALL CUNT03( 'C', M, M, M, MNMIN, USAV, LDU, 610 $ U, LDU, WORK, LWORK, RWORK, DIF, 611 $ IINFO ) 612 END IF 613 END IF 614 RESULT( 5 ) = MAX( RESULT( 5 ), DIF ) 615* 616* Compare VT 617* 618 DIF = ZERO 619 IF( M.GT.0 .AND. N.GT.0 ) THEN 620 IF( IJVT.EQ.1 ) THEN 621 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV, 622 $ LDVT, A, LDA, WORK, LWORK, 623 $ RWORK, DIF, IINFO ) 624 ELSE IF( IJVT.EQ.2 ) THEN 625 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV, 626 $ LDVT, VT, LDVT, WORK, LWORK, 627 $ RWORK, DIF, IINFO ) 628 ELSE IF( IJVT.EQ.3 ) THEN 629 CALL CUNT03( 'R', N, N, N, MNMIN, VTSAV, 630 $ LDVT, VT, LDVT, WORK, LWORK, 631 $ RWORK, DIF, IINFO ) 632 END IF 633 END IF 634 RESULT( 6 ) = MAX( RESULT( 6 ), DIF ) 635* 636* Compare S 637* 638 DIF = ZERO 639 DIV = MAX( REAL( MNMIN )*ULP*S( 1 ), 640 $ SLAMCH( 'Safe minimum' ) ) 641 DO 80 I = 1, MNMIN - 1 642 IF( SSAV( I ).LT.SSAV( I+1 ) ) 643 $ DIF = ULPINV 644 IF( SSAV( I ).LT.ZERO ) 645 $ DIF = ULPINV 646 DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV ) 647 80 CONTINUE 648 RESULT( 7 ) = MAX( RESULT( 7 ), DIF ) 649 90 CONTINUE 650 100 CONTINUE 651* 652* Test for CGESDD 653* 654 IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N ) 655 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3 656 LSWORK = MIN( LSWORK, LWORK ) 657 LSWORK = MAX( LSWORK, 1 ) 658 IF( IWSPC.EQ.4 ) 659 $ LSWORK = LWORK 660* 661* Factorize A 662* 663 CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA ) 664 CALL CGESDD( 'A', M, N, A, LDA, SSAV, USAV, LDU, VTSAV, 665 $ LDVT, WORK, LSWORK, RWORK, IWORK, IINFO ) 666 IF( IINFO.NE.0 ) THEN 667 WRITE( NOUNIT, FMT = 9995 )'GESDD', IINFO, M, N, 668 $ JTYPE, LSWORK, IOLDSD 669 INFO = ABS( IINFO ) 670 RETURN 671 END IF 672* 673* Do tests 1--4 674* 675 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E, 676 $ VTSAV, LDVT, WORK, RWORK, RESULT( 8 ) ) 677 IF( M.NE.0 .AND. N.NE.0 ) THEN 678 CALL CUNT01( 'Columns', MNMIN, M, USAV, LDU, WORK, 679 $ LWORK, RWORK, RESULT( 9 ) ) 680 CALL CUNT01( 'Rows', MNMIN, N, VTSAV, LDVT, WORK, 681 $ LWORK, RWORK, RESULT( 10 ) ) 682 END IF 683 RESULT( 11 ) = 0 684 DO 110 I = 1, MNMIN - 1 685 IF( SSAV( I ).LT.SSAV( I+1 ) ) 686 $ RESULT( 11 ) = ULPINV 687 IF( SSAV( I ).LT.ZERO ) 688 $ RESULT( 11 ) = ULPINV 689 110 CONTINUE 690 IF( MNMIN.GE.1 ) THEN 691 IF( SSAV( MNMIN ).LT.ZERO ) 692 $ RESULT( 11 ) = ULPINV 693 END IF 694* 695* Do partial SVDs, comparing to SSAV, USAV, and VTSAV 696* 697 RESULT( 12 ) = ZERO 698 RESULT( 13 ) = ZERO 699 RESULT( 14 ) = ZERO 700 DO 130 IJQ = 0, 2 701 JOBQ = CJOB( IJQ+1 ) 702 CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA ) 703 CALL CGESDD( JOBQ, M, N, A, LDA, S, U, LDU, VT, LDVT, 704 $ WORK, LSWORK, RWORK, IWORK, IINFO ) 705* 706* Compare U 707* 708 DIF = ZERO 709 IF( M.GT.0 .AND. N.GT.0 ) THEN 710 IF( IJQ.EQ.1 ) THEN 711 IF( M.GE.N ) THEN 712 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV, 713 $ LDU, A, LDA, WORK, LWORK, RWORK, 714 $ DIF, IINFO ) 715 ELSE 716 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV, 717 $ LDU, U, LDU, WORK, LWORK, RWORK, 718 $ DIF, IINFO ) 719 END IF 720 ELSE IF( IJQ.EQ.2 ) THEN 721 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV, LDU, 722 $ U, LDU, WORK, LWORK, RWORK, DIF, 723 $ IINFO ) 724 END IF 725 END IF 726 RESULT( 12 ) = MAX( RESULT( 12 ), DIF ) 727* 728* Compare VT 729* 730 DIF = ZERO 731 IF( M.GT.0 .AND. N.GT.0 ) THEN 732 IF( IJQ.EQ.1 ) THEN 733 IF( M.GE.N ) THEN 734 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV, 735 $ LDVT, VT, LDVT, WORK, LWORK, 736 $ RWORK, DIF, IINFO ) 737 ELSE 738 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV, 739 $ LDVT, A, LDA, WORK, LWORK, 740 $ RWORK, DIF, IINFO ) 741 END IF 742 ELSE IF( IJQ.EQ.2 ) THEN 743 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV, 744 $ LDVT, VT, LDVT, WORK, LWORK, RWORK, 745 $ DIF, IINFO ) 746 END IF 747 END IF 748 RESULT( 13 ) = MAX( RESULT( 13 ), DIF ) 749* 750* Compare S 751* 752 DIF = ZERO 753 DIV = MAX( REAL( MNMIN )*ULP*S( 1 ), 754 $ SLAMCH( 'Safe minimum' ) ) 755 DO 120 I = 1, MNMIN - 1 756 IF( SSAV( I ).LT.SSAV( I+1 ) ) 757 $ DIF = ULPINV 758 IF( SSAV( I ).LT.ZERO ) 759 $ DIF = ULPINV 760 DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV ) 761 120 CONTINUE 762 RESULT( 14 ) = MAX( RESULT( 14 ), DIF ) 763 130 CONTINUE 764* 765* End of Loop -- Check for RESULT(j) > THRESH 766* 767 NTEST = 0 768 NFAIL = 0 769 DO 140 J = 1, 14 770 IF( RESULT( J ).GE.ZERO ) 771 $ NTEST = NTEST + 1 772 IF( RESULT( J ).GE.THRESH ) 773 $ NFAIL = NFAIL + 1 774 140 CONTINUE 775* 776 IF( NFAIL.GT.0 ) 777 $ NTESTF = NTESTF + 1 778 IF( NTESTF.EQ.1 ) THEN 779 WRITE( NOUNIT, FMT = 9999 ) 780 WRITE( NOUNIT, FMT = 9998 )THRESH 781 NTESTF = 2 782 END IF 783* 784 DO 150 J = 1, 14 785 IF( RESULT( J ).GE.THRESH ) THEN 786 WRITE( NOUNIT, FMT = 9997 )M, N, JTYPE, IWSPC, 787 $ IOLDSD, J, RESULT( J ) 788 END IF 789 150 CONTINUE 790* 791 NERRS = NERRS + NFAIL 792 NTESTT = NTESTT + NTEST 793* 794 160 CONTINUE 795* 796 170 CONTINUE 797 180 CONTINUE 798* 799* Summary 800* 801 CALL ALASVM( 'CBD', NOUNIT, NERRS, NTESTT, 0 ) 802* 803 9999 FORMAT( ' SVD -- Complex Singular Value Decomposition Driver ', 804 $ / ' Matrix types (see CDRVBD for details):', 805 $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix', 806 $ / ' 3 = Evenly spaced singular values near 1', 807 $ / ' 4 = Evenly spaced singular values near underflow', 808 $ / ' 5 = Evenly spaced singular values near overflow', 809 $ / / ' Tests performed: ( A is dense, U and V are unitary,', 810 $ / 19X, ' S is an array, and Upartial, VTpartial, and', 811 $ / 19X, ' Spartial are partially computed U, VT and S),', / ) 812 9998 FORMAT( ' Tests performed with Test Threshold = ', F8.2, 813 $ / ' CGESVD: ', / 814 $ ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ', 815 $ / ' 2 = | I - U**T U | / ( M ulp ) ', 816 $ / ' 3 = | I - VT VT**T | / ( N ulp ) ', 817 $ / ' 4 = 0 if S contains min(M,N) nonnegative values in', 818 $ ' decreasing order, else 1/ulp', 819 $ / ' 5 = | U - Upartial | / ( M ulp )', 820 $ / ' 6 = | VT - VTpartial | / ( N ulp )', 821 $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )', 822 $ / ' CGESDD: ', / 823 $ ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ', 824 $ / ' 9 = | I - U**T U | / ( M ulp ) ', 825 $ / '10 = | I - VT VT**T | / ( N ulp ) ', 826 $ / '11 = 0 if S contains min(M,N) nonnegative values in', 827 $ ' decreasing order, else 1/ulp', 828 $ / '12 = | U - Upartial | / ( M ulp )', 829 $ / '13 = | VT - VTpartial | / ( N ulp )', 830 $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )', / / ) 831 9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1, 832 $ ', seed=', 4( I4, ',' ), ' test(', I1, ')=', G11.4 ) 833 9996 FORMAT( ' CDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=', 834 $ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), 835 $ I5, ')' ) 836 9995 FORMAT( ' CDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=', 837 $ I6, ', N=', I6, ', JTYPE=', I6, ', LSWORK=', I6, / 9X, 838 $ 'ISEED=(', 3( I5, ',' ), I5, ')' ) 839* 840 RETURN 841* 842* End of CDRVBD 843* 844 END 845