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