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