|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH, 00002 $ A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S, 00003 $ SSAV, E, WORK, LWORK, RWORK, IWORK, NOUNIT, 00004 $ INFO ) 00005 * 00006 * -- LAPACK test routine (version 3.1) -- 00007 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES, 00012 $ NTYPES 00013 DOUBLE PRECISION THRESH 00014 * .. 00015 * .. Array Arguments .. 00016 LOGICAL DOTYPE( * ) 00017 INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * ) 00018 DOUBLE PRECISION E( * ), RWORK( * ), S( * ), SSAV( * ) 00019 COMPLEX*16 A( LDA, * ), ASAV( LDA, * ), U( LDU, * ), 00020 $ USAV( LDU, * ), VT( LDVT, * ), 00021 $ VTSAV( LDVT, * ), WORK( * ) 00022 * .. 00023 * 00024 * Purpose 00025 * ======= 00026 * 00027 * ZDRVBD checks the singular value decomposition (SVD) driver ZGESVD 00028 * and ZGESDD. 00029 * ZGESVD and CGESDD factors A = U diag(S) VT, where U and VT are 00030 * unitary and diag(S) is diagonal with the entries of the array S on 00031 * its diagonal. The entries of S are the singular values, nonnegative 00032 * and stored in decreasing order. U and VT can be optionally not 00033 * computed, overwritten on A, or computed partially. 00034 * 00035 * A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN. 00036 * U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N. 00037 * 00038 * When ZDRVBD is called, a number of matrix "sizes" (M's and N's) 00039 * and a number of matrix "types" are specified. For each size (M,N) 00040 * and each type of matrix, and for the minimal workspace as well as 00041 * workspace adequate to permit blocking, an M x N matrix "A" will be 00042 * generated and used to test the SVD routines. For each matrix, A will 00043 * be factored as A = U diag(S) VT and the following 12 tests computed: 00044 * 00045 * Test for ZGESVD: 00046 * 00047 * (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp ) 00048 * 00049 * (2) | I - U'U | / ( M ulp ) 00050 * 00051 * (3) | I - VT VT' | / ( N ulp ) 00052 * 00053 * (4) S contains MNMIN nonnegative values in decreasing order. 00054 * (Return 0 if true, 1/ULP if false.) 00055 * 00056 * (5) | U - Upartial | / ( M ulp ) where Upartial is a partially 00057 * computed U. 00058 * 00059 * (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially 00060 * computed VT. 00061 * 00062 * (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the 00063 * vector of singular values from the partial SVD 00064 * 00065 * Test for ZGESDD: 00066 * 00067 * (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp ) 00068 * 00069 * (2) | I - U'U | / ( M ulp ) 00070 * 00071 * (3) | I - VT VT' | / ( N ulp ) 00072 * 00073 * (4) S contains MNMIN nonnegative values in decreasing order. 00074 * (Return 0 if true, 1/ULP if false.) 00075 * 00076 * (5) | U - Upartial | / ( M ulp ) where Upartial is a partially 00077 * computed U. 00078 * 00079 * (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially 00080 * computed VT. 00081 * 00082 * (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the 00083 * vector of singular values from the partial SVD 00084 * 00085 * The "sizes" are specified by the arrays MM(1:NSIZES) and 00086 * NN(1:NSIZES); the value of each element pair (MM(j),NN(j)) 00087 * specifies one size. The "types" are specified by a logical array 00088 * DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j" 00089 * will be generated. 00090 * Currently, the list of possible types is: 00091 * 00092 * (1) The zero matrix. 00093 * (2) The identity matrix. 00094 * (3) A matrix of the form U D V, where U and V are unitary and 00095 * D has evenly spaced entries 1, ..., ULP with random signs 00096 * on the diagonal. 00097 * (4) Same as (3), but multiplied by the underflow-threshold / ULP. 00098 * (5) Same as (3), but multiplied by the overflow-threshold * ULP. 00099 * 00100 * Arguments 00101 * ========== 00102 * 00103 * NSIZES (input) INTEGER 00104 * The number of sizes of matrices to use. If it is zero, 00105 * ZDRVBD does nothing. It must be at least zero. 00106 * 00107 * MM (input) INTEGER array, dimension (NSIZES) 00108 * An array containing the matrix "heights" to be used. For 00109 * each j=1,...,NSIZES, if MM(j) is zero, then MM(j) and NN(j) 00110 * will be ignored. The MM(j) values must be at least zero. 00111 * 00112 * NN (input) INTEGER array, dimension (NSIZES) 00113 * An array containing the matrix "widths" to be used. For 00114 * each j=1,...,NSIZES, if NN(j) is zero, then MM(j) and NN(j) 00115 * will be ignored. The NN(j) values must be at least zero. 00116 * 00117 * NTYPES (input) INTEGER 00118 * The number of elements in DOTYPE. If it is zero, ZDRVBD 00119 * does nothing. It must be at least zero. If it is MAXTYP+1 00120 * and NSIZES is 1, then an additional type, MAXTYP+1 is 00121 * defined, which is to use whatever matrices are in A and B. 00122 * This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 00123 * DOTYPE(MAXTYP+1) is .TRUE. . 00124 * 00125 * DOTYPE (input) LOGICAL array, dimension (NTYPES) 00126 * If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix 00127 * of type j will be generated. If NTYPES is smaller than the 00128 * maximum number of types defined (PARAMETER MAXTYP), then 00129 * types NTYPES+1 through MAXTYP will not be generated. If 00130 * NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through 00131 * DOTYPE(NTYPES) will be ignored. 00132 * 00133 * ISEED (input/output) INTEGER array, dimension (4) 00134 * On entry ISEED specifies the seed of the random number 00135 * generator. The array elements should be between 0 and 4095; 00136 * if not they will be reduced mod 4096. Also, ISEED(4) must 00137 * be odd. The random number generator uses a linear 00138 * congruential sequence limited to small integers, and so 00139 * should produce machine independent random numbers. The 00140 * values of ISEED are changed on exit, and can be used in the 00141 * next call to ZDRVBD to continue the same random number 00142 * sequence. 00143 * 00144 * THRESH (input) DOUBLE PRECISION 00145 * A test will count as "failed" if the "error", computed as 00146 * described above, exceeds THRESH. Note that the error 00147 * is scaled to be O(1), so THRESH should be a reasonably 00148 * small multiple of 1, e.g., 10 or 100. In particular, 00149 * it should not depend on the precision (single vs. double) 00150 * or the size of the matrix. It must be at least zero. 00151 * 00152 * NOUNIT (input) INTEGER 00153 * The FORTRAN unit number for printing out error messages 00154 * (e.g., if a routine returns IINFO not equal to 0.) 00155 * 00156 * A (output) COMPLEX*16 array, dimension (LDA,max(NN)) 00157 * Used to hold the matrix whose singular values are to be 00158 * computed. On exit, A contains the last matrix actually 00159 * used. 00160 * 00161 * LDA (input) INTEGER 00162 * The leading dimension of A. It must be at 00163 * least 1 and at least max( MM ). 00164 * 00165 * U (output) COMPLEX*16 array, dimension (LDU,max(MM)) 00166 * Used to hold the computed matrix of right singular vectors. 00167 * On exit, U contains the last such vectors actually computed. 00168 * 00169 * LDU (input) INTEGER 00170 * The leading dimension of U. It must be at 00171 * least 1 and at least max( MM ). 00172 * 00173 * VT (output) COMPLEX*16 array, dimension (LDVT,max(NN)) 00174 * Used to hold the computed matrix of left singular vectors. 00175 * On exit, VT contains the last such vectors actually computed. 00176 * 00177 * LDVT (input) INTEGER 00178 * The leading dimension of VT. It must be at 00179 * least 1 and at least max( NN ). 00180 * 00181 * ASAV (output) COMPLEX*16 array, dimension (LDA,max(NN)) 00182 * Used to hold a different copy of the matrix whose singular 00183 * values are to be computed. On exit, A contains the last 00184 * matrix actually used. 00185 * 00186 * USAV (output) COMPLEX*16 array, dimension (LDU,max(MM)) 00187 * Used to hold a different copy of the computed matrix of 00188 * right singular vectors. On exit, USAV contains the last such 00189 * vectors actually computed. 00190 * 00191 * VTSAV (output) COMPLEX*16 array, dimension (LDVT,max(NN)) 00192 * Used to hold a different copy of the computed matrix of 00193 * left singular vectors. On exit, VTSAV contains the last such 00194 * vectors actually computed. 00195 * 00196 * S (output) DOUBLE PRECISION array, dimension (max(min(MM,NN))) 00197 * Contains the computed singular values. 00198 * 00199 * SSAV (output) DOUBLE PRECISION array, dimension (max(min(MM,NN))) 00200 * Contains another copy of the computed singular values. 00201 * 00202 * E (output) DOUBLE PRECISION array, dimension (max(min(MM,NN))) 00203 * Workspace for ZGESVD. 00204 * 00205 * WORK (workspace) COMPLEX*16 array, dimension (LWORK) 00206 * 00207 * LWORK (input) INTEGER 00208 * The number of entries in WORK. This must be at least 00209 * MAX(3*MIN(M,N)+MAX(M,N)**2,5*MIN(M,N),3*MAX(M,N)) for all 00210 * pairs (M,N)=(MM(j),NN(j)) 00211 * 00212 * RWORK (workspace) DOUBLE PRECISION array, 00213 * dimension ( 5*max(max(MM,NN)) ) 00214 * 00215 * IWORK (workspace) INTEGER array, dimension at least 8*min(M,N) 00216 * 00217 * RESULT (output) DOUBLE PRECISION array, dimension (7) 00218 * The values computed by the 7 tests described above. 00219 * The values are currently limited to 1/ULP, to avoid 00220 * overflow. 00221 * 00222 * INFO (output) INTEGER 00223 * If 0, then everything ran OK. 00224 * -1: NSIZES < 0 00225 * -2: Some MM(j) < 0 00226 * -3: Some NN(j) < 0 00227 * -4: NTYPES < 0 00228 * -7: THRESH < 0 00229 * -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ). 00230 * -12: LDU < 1 or LDU < MMAX. 00231 * -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ). 00232 * -21: LWORK too small. 00233 * If ZLATMS, or ZGESVD returns an error code, the 00234 * absolute value of it is returned. 00235 * 00236 * ===================================================================== 00237 * 00238 * .. Parameters .. 00239 DOUBLE PRECISION ZERO, ONE 00240 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00241 COMPLEX*16 CZERO, CONE 00242 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00243 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00244 INTEGER MAXTYP 00245 PARAMETER ( MAXTYP = 5 ) 00246 * .. 00247 * .. Local Scalars .. 00248 LOGICAL BADMM, BADNN 00249 CHARACTER JOBQ, JOBU, JOBVT 00250 INTEGER I, IINFO, IJQ, IJU, IJVT, IWSPC, IWTMP, J, 00251 $ JSIZE, JTYPE, LSWORK, M, MINWRK, MMAX, MNMAX, 00252 $ MNMIN, MTYPES, N, NERRS, NFAIL, NMAX, NTEST, 00253 $ NTESTF, NTESTT 00254 DOUBLE PRECISION ANORM, DIF, DIV, OVFL, ULP, ULPINV, UNFL 00255 * .. 00256 * .. Local Arrays .. 00257 CHARACTER CJOB( 4 ) 00258 INTEGER IOLDSD( 4 ) 00259 DOUBLE PRECISION RESULT( 14 ) 00260 * .. 00261 * .. External Functions .. 00262 DOUBLE PRECISION DLAMCH 00263 EXTERNAL DLAMCH 00264 * .. 00265 * .. External Subroutines .. 00266 EXTERNAL ALASVM, XERBLA, ZBDT01, ZGESDD, ZGESVD, ZLACPY, 00267 $ ZLASET, ZLATMS, ZUNT01, ZUNT03 00268 * .. 00269 * .. Intrinsic Functions .. 00270 INTRINSIC ABS, DBLE, MAX, MIN 00271 * .. 00272 * .. Data statements .. 00273 DATA CJOB / 'N', 'O', 'S', 'A' / 00274 * .. 00275 * .. Executable Statements .. 00276 * 00277 * Check for errors 00278 * 00279 INFO = 0 00280 * 00281 * Important constants 00282 * 00283 NERRS = 0 00284 NTESTT = 0 00285 NTESTF = 0 00286 BADMM = .FALSE. 00287 BADNN = .FALSE. 00288 MMAX = 1 00289 NMAX = 1 00290 MNMAX = 1 00291 MINWRK = 1 00292 DO 10 J = 1, NSIZES 00293 MMAX = MAX( MMAX, MM( J ) ) 00294 IF( MM( J ).LT.0 ) 00295 $ BADMM = .TRUE. 00296 NMAX = MAX( NMAX, NN( J ) ) 00297 IF( NN( J ).LT.0 ) 00298 $ BADNN = .TRUE. 00299 MNMAX = MAX( MNMAX, MIN( MM( J ), NN( J ) ) ) 00300 MINWRK = MAX( MINWRK, MAX( 3*MIN( MM( J ), 00301 $ NN( J ) )+MAX( MM( J ), NN( J ) )**2, 5*MIN( MM( J ), 00302 $ NN( J ) ), 3*MAX( MM( J ), NN( J ) ) ) ) 00303 10 CONTINUE 00304 * 00305 * Check for errors 00306 * 00307 IF( NSIZES.LT.0 ) THEN 00308 INFO = -1 00309 ELSE IF( BADMM ) THEN 00310 INFO = -2 00311 ELSE IF( BADNN ) THEN 00312 INFO = -3 00313 ELSE IF( NTYPES.LT.0 ) THEN 00314 INFO = -4 00315 ELSE IF( LDA.LT.MAX( 1, MMAX ) ) THEN 00316 INFO = -10 00317 ELSE IF( LDU.LT.MAX( 1, MMAX ) ) THEN 00318 INFO = -12 00319 ELSE IF( LDVT.LT.MAX( 1, NMAX ) ) THEN 00320 INFO = -14 00321 ELSE IF( MINWRK.GT.LWORK ) THEN 00322 INFO = -21 00323 END IF 00324 * 00325 IF( INFO.NE.0 ) THEN 00326 CALL XERBLA( 'ZDRVBD', -INFO ) 00327 RETURN 00328 END IF 00329 * 00330 * Quick return if nothing to do 00331 * 00332 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 ) 00333 $ RETURN 00334 * 00335 * More Important constants 00336 * 00337 UNFL = DLAMCH( 'S' ) 00338 OVFL = ONE / UNFL 00339 ULP = DLAMCH( 'E' ) 00340 ULPINV = ONE / ULP 00341 * 00342 * Loop over sizes, types 00343 * 00344 NERRS = 0 00345 * 00346 DO 180 JSIZE = 1, NSIZES 00347 M = MM( JSIZE ) 00348 N = NN( JSIZE ) 00349 MNMIN = MIN( M, N ) 00350 * 00351 IF( NSIZES.NE.1 ) THEN 00352 MTYPES = MIN( MAXTYP, NTYPES ) 00353 ELSE 00354 MTYPES = MIN( MAXTYP+1, NTYPES ) 00355 END IF 00356 * 00357 DO 170 JTYPE = 1, MTYPES 00358 IF( .NOT.DOTYPE( JTYPE ) ) 00359 $ GO TO 170 00360 NTEST = 0 00361 * 00362 DO 20 J = 1, 4 00363 IOLDSD( J ) = ISEED( J ) 00364 20 CONTINUE 00365 * 00366 * Compute "A" 00367 * 00368 IF( MTYPES.GT.MAXTYP ) 00369 $ GO TO 50 00370 * 00371 IF( JTYPE.EQ.1 ) THEN 00372 * 00373 * Zero matrix 00374 * 00375 CALL ZLASET( 'Full', M, N, CZERO, CZERO, A, LDA ) 00376 DO 30 I = 1, MIN( M, N ) 00377 S( I ) = ZERO 00378 30 CONTINUE 00379 * 00380 ELSE IF( JTYPE.EQ.2 ) THEN 00381 * 00382 * Identity matrix 00383 * 00384 CALL ZLASET( 'Full', M, N, CZERO, CONE, A, LDA ) 00385 DO 40 I = 1, MIN( M, N ) 00386 S( I ) = ONE 00387 40 CONTINUE 00388 * 00389 ELSE 00390 * 00391 * (Scaled) random matrix 00392 * 00393 IF( JTYPE.EQ.3 ) 00394 $ ANORM = ONE 00395 IF( JTYPE.EQ.4 ) 00396 $ ANORM = UNFL / ULP 00397 IF( JTYPE.EQ.5 ) 00398 $ ANORM = OVFL*ULP 00399 CALL ZLATMS( M, N, 'U', ISEED, 'N', S, 4, DBLE( MNMIN ), 00400 $ ANORM, M-1, N-1, 'N', A, LDA, WORK, IINFO ) 00401 IF( IINFO.NE.0 ) THEN 00402 WRITE( NOUNIT, FMT = 9996 )'Generator', IINFO, M, N, 00403 $ JTYPE, IOLDSD 00404 INFO = ABS( IINFO ) 00405 RETURN 00406 END IF 00407 END IF 00408 * 00409 50 CONTINUE 00410 CALL ZLACPY( 'F', M, N, A, LDA, ASAV, LDA ) 00411 * 00412 * Do for minimal and adequate (for blocking) workspace 00413 * 00414 DO 160 IWSPC = 1, 4 00415 * 00416 * Test for ZGESVD 00417 * 00418 IWTMP = 2*MIN( M, N )+MAX( M, N ) 00419 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3 00420 LSWORK = MIN( LSWORK, LWORK ) 00421 LSWORK = MAX( LSWORK, 1 ) 00422 IF( IWSPC.EQ.4 ) 00423 $ LSWORK = LWORK 00424 * 00425 DO 60 J = 1, 14 00426 RESULT( J ) = -ONE 00427 60 CONTINUE 00428 * 00429 * Factorize A 00430 * 00431 IF( IWSPC.GT.1 ) 00432 $ CALL ZLACPY( 'F', M, N, ASAV, LDA, A, LDA ) 00433 CALL ZGESVD( 'A', 'A', M, N, A, LDA, SSAV, USAV, LDU, 00434 $ VTSAV, LDVT, WORK, LSWORK, RWORK, IINFO ) 00435 IF( IINFO.NE.0 ) THEN 00436 WRITE( NOUNIT, FMT = 9995 )'GESVD', IINFO, M, N, 00437 $ JTYPE, LSWORK, IOLDSD 00438 INFO = ABS( IINFO ) 00439 RETURN 00440 END IF 00441 * 00442 * Do tests 1--4 00443 * 00444 CALL ZBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E, 00445 $ VTSAV, LDVT, WORK, RWORK, RESULT( 1 ) ) 00446 IF( M.NE.0 .AND. N.NE.0 ) THEN 00447 CALL ZUNT01( 'Columns', MNMIN, M, USAV, LDU, WORK, 00448 $ LWORK, RWORK, RESULT( 2 ) ) 00449 CALL ZUNT01( 'Rows', MNMIN, N, VTSAV, LDVT, WORK, 00450 $ LWORK, RWORK, RESULT( 3 ) ) 00451 END IF 00452 RESULT( 4 ) = 0 00453 DO 70 I = 1, MNMIN - 1 00454 IF( SSAV( I ).LT.SSAV( I+1 ) ) 00455 $ RESULT( 4 ) = ULPINV 00456 IF( SSAV( I ).LT.ZERO ) 00457 $ RESULT( 4 ) = ULPINV 00458 70 CONTINUE 00459 IF( MNMIN.GE.1 ) THEN 00460 IF( SSAV( MNMIN ).LT.ZERO ) 00461 $ RESULT( 4 ) = ULPINV 00462 END IF 00463 * 00464 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV 00465 * 00466 RESULT( 5 ) = ZERO 00467 RESULT( 6 ) = ZERO 00468 RESULT( 7 ) = ZERO 00469 DO 100 IJU = 0, 3 00470 DO 90 IJVT = 0, 3 00471 IF( ( IJU.EQ.3 .AND. IJVT.EQ.3 ) .OR. 00472 $ ( IJU.EQ.1 .AND. IJVT.EQ.1 ) )GO TO 90 00473 JOBU = CJOB( IJU+1 ) 00474 JOBVT = CJOB( IJVT+1 ) 00475 CALL ZLACPY( 'F', M, N, ASAV, LDA, A, LDA ) 00476 CALL ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, 00477 $ VT, LDVT, WORK, LSWORK, RWORK, IINFO ) 00478 * 00479 * Compare U 00480 * 00481 DIF = ZERO 00482 IF( M.GT.0 .AND. N.GT.0 ) THEN 00483 IF( IJU.EQ.1 ) THEN 00484 CALL ZUNT03( 'C', M, MNMIN, M, MNMIN, USAV, 00485 $ LDU, A, LDA, WORK, LWORK, RWORK, 00486 $ DIF, IINFO ) 00487 ELSE IF( IJU.EQ.2 ) THEN 00488 CALL ZUNT03( 'C', M, MNMIN, M, MNMIN, USAV, 00489 $ LDU, U, LDU, WORK, LWORK, RWORK, 00490 $ DIF, IINFO ) 00491 ELSE IF( IJU.EQ.3 ) THEN 00492 CALL ZUNT03( 'C', M, M, M, MNMIN, USAV, LDU, 00493 $ U, LDU, WORK, LWORK, RWORK, DIF, 00494 $ IINFO ) 00495 END IF 00496 END IF 00497 RESULT( 5 ) = MAX( RESULT( 5 ), DIF ) 00498 * 00499 * Compare VT 00500 * 00501 DIF = ZERO 00502 IF( M.GT.0 .AND. N.GT.0 ) THEN 00503 IF( IJVT.EQ.1 ) THEN 00504 CALL ZUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV, 00505 $ LDVT, A, LDA, WORK, LWORK, 00506 $ RWORK, DIF, IINFO ) 00507 ELSE IF( IJVT.EQ.2 ) THEN 00508 CALL ZUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV, 00509 $ LDVT, VT, LDVT, WORK, LWORK, 00510 $ RWORK, DIF, IINFO ) 00511 ELSE IF( IJVT.EQ.3 ) THEN 00512 CALL ZUNT03( 'R', N, N, N, MNMIN, VTSAV, 00513 $ LDVT, VT, LDVT, WORK, LWORK, 00514 $ RWORK, DIF, IINFO ) 00515 END IF 00516 END IF 00517 RESULT( 6 ) = MAX( RESULT( 6 ), DIF ) 00518 * 00519 * Compare S 00520 * 00521 DIF = ZERO 00522 DIV = MAX( DBLE( MNMIN )*ULP*S( 1 ), 00523 $ DLAMCH( 'Safe minimum' ) ) 00524 DO 80 I = 1, MNMIN - 1 00525 IF( SSAV( I ).LT.SSAV( I+1 ) ) 00526 $ DIF = ULPINV 00527 IF( SSAV( I ).LT.ZERO ) 00528 $ DIF = ULPINV 00529 DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV ) 00530 80 CONTINUE 00531 RESULT( 7 ) = MAX( RESULT( 7 ), DIF ) 00532 90 CONTINUE 00533 100 CONTINUE 00534 * 00535 * Test for ZGESDD 00536 * 00537 IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N ) 00538 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3 00539 LSWORK = MIN( LSWORK, LWORK ) 00540 LSWORK = MAX( LSWORK, 1 ) 00541 IF( IWSPC.EQ.4 ) 00542 $ LSWORK = LWORK 00543 * 00544 * Factorize A 00545 * 00546 CALL ZLACPY( 'F', M, N, ASAV, LDA, A, LDA ) 00547 CALL ZGESDD( 'A', M, N, A, LDA, SSAV, USAV, LDU, VTSAV, 00548 $ LDVT, WORK, LSWORK, RWORK, IWORK, IINFO ) 00549 IF( IINFO.NE.0 ) THEN 00550 WRITE( NOUNIT, FMT = 9995 )'GESDD', IINFO, M, N, 00551 $ JTYPE, LSWORK, IOLDSD 00552 INFO = ABS( IINFO ) 00553 RETURN 00554 END IF 00555 * 00556 * Do tests 1--4 00557 * 00558 CALL ZBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E, 00559 $ VTSAV, LDVT, WORK, RWORK, RESULT( 8 ) ) 00560 IF( M.NE.0 .AND. N.NE.0 ) THEN 00561 CALL ZUNT01( 'Columns', MNMIN, M, USAV, LDU, WORK, 00562 $ LWORK, RWORK, RESULT( 9 ) ) 00563 CALL ZUNT01( 'Rows', MNMIN, N, VTSAV, LDVT, WORK, 00564 $ LWORK, RWORK, RESULT( 10 ) ) 00565 END IF 00566 RESULT( 11 ) = 0 00567 DO 110 I = 1, MNMIN - 1 00568 IF( SSAV( I ).LT.SSAV( I+1 ) ) 00569 $ RESULT( 11 ) = ULPINV 00570 IF( SSAV( I ).LT.ZERO ) 00571 $ RESULT( 11 ) = ULPINV 00572 110 CONTINUE 00573 IF( MNMIN.GE.1 ) THEN 00574 IF( SSAV( MNMIN ).LT.ZERO ) 00575 $ RESULT( 11 ) = ULPINV 00576 END IF 00577 * 00578 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV 00579 * 00580 RESULT( 12 ) = ZERO 00581 RESULT( 13 ) = ZERO 00582 RESULT( 14 ) = ZERO 00583 DO 130 IJQ = 0, 2 00584 JOBQ = CJOB( IJQ+1 ) 00585 CALL ZLACPY( 'F', M, N, ASAV, LDA, A, LDA ) 00586 CALL ZGESDD( JOBQ, M, N, A, LDA, S, U, LDU, VT, LDVT, 00587 $ WORK, LSWORK, RWORK, IWORK, IINFO ) 00588 * 00589 * Compare U 00590 * 00591 DIF = ZERO 00592 IF( M.GT.0 .AND. N.GT.0 ) THEN 00593 IF( IJQ.EQ.1 ) THEN 00594 IF( M.GE.N ) THEN 00595 CALL ZUNT03( 'C', M, MNMIN, M, MNMIN, USAV, 00596 $ LDU, A, LDA, WORK, LWORK, RWORK, 00597 $ DIF, IINFO ) 00598 ELSE 00599 CALL ZUNT03( 'C', M, MNMIN, M, MNMIN, USAV, 00600 $ LDU, U, LDU, WORK, LWORK, RWORK, 00601 $ DIF, IINFO ) 00602 END IF 00603 ELSE IF( IJQ.EQ.2 ) THEN 00604 CALL ZUNT03( 'C', M, MNMIN, M, MNMIN, USAV, LDU, 00605 $ U, LDU, WORK, LWORK, RWORK, DIF, 00606 $ IINFO ) 00607 END IF 00608 END IF 00609 RESULT( 12 ) = MAX( RESULT( 12 ), DIF ) 00610 * 00611 * Compare VT 00612 * 00613 DIF = ZERO 00614 IF( M.GT.0 .AND. N.GT.0 ) THEN 00615 IF( IJQ.EQ.1 ) THEN 00616 IF( M.GE.N ) THEN 00617 CALL ZUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV, 00618 $ LDVT, VT, LDVT, WORK, LWORK, 00619 $ RWORK, DIF, IINFO ) 00620 ELSE 00621 CALL ZUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV, 00622 $ LDVT, A, LDA, WORK, LWORK, 00623 $ RWORK, DIF, IINFO ) 00624 END IF 00625 ELSE IF( IJQ.EQ.2 ) THEN 00626 CALL ZUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV, 00627 $ LDVT, VT, LDVT, WORK, LWORK, RWORK, 00628 $ DIF, IINFO ) 00629 END IF 00630 END IF 00631 RESULT( 13 ) = MAX( RESULT( 13 ), DIF ) 00632 * 00633 * Compare S 00634 * 00635 DIF = ZERO 00636 DIV = MAX( DBLE( MNMIN )*ULP*S( 1 ), 00637 $ DLAMCH( 'Safe minimum' ) ) 00638 DO 120 I = 1, MNMIN - 1 00639 IF( SSAV( I ).LT.SSAV( I+1 ) ) 00640 $ DIF = ULPINV 00641 IF( SSAV( I ).LT.ZERO ) 00642 $ DIF = ULPINV 00643 DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV ) 00644 120 CONTINUE 00645 RESULT( 14 ) = MAX( RESULT( 14 ), DIF ) 00646 130 CONTINUE 00647 * 00648 * End of Loop -- Check for RESULT(j) > THRESH 00649 * 00650 NTEST = 0 00651 NFAIL = 0 00652 DO 140 J = 1, 14 00653 IF( RESULT( J ).GE.ZERO ) 00654 $ NTEST = NTEST + 1 00655 IF( RESULT( J ).GE.THRESH ) 00656 $ NFAIL = NFAIL + 1 00657 140 CONTINUE 00658 * 00659 IF( NFAIL.GT.0 ) 00660 $ NTESTF = NTESTF + 1 00661 IF( NTESTF.EQ.1 ) THEN 00662 WRITE( NOUNIT, FMT = 9999 ) 00663 WRITE( NOUNIT, FMT = 9998 )THRESH 00664 NTESTF = 2 00665 END IF 00666 * 00667 DO 150 J = 1, 14 00668 IF( RESULT( J ).GE.THRESH ) THEN 00669 WRITE( NOUNIT, FMT = 9997 )M, N, JTYPE, IWSPC, 00670 $ IOLDSD, J, RESULT( J ) 00671 END IF 00672 150 CONTINUE 00673 * 00674 NERRS = NERRS + NFAIL 00675 NTESTT = NTESTT + NTEST 00676 * 00677 160 CONTINUE 00678 * 00679 170 CONTINUE 00680 180 CONTINUE 00681 * 00682 * Summary 00683 * 00684 CALL ALASVM( 'ZBD', NOUNIT, NERRS, NTESTT, 0 ) 00685 * 00686 9999 FORMAT( ' SVD -- Complex Singular Value Decomposition Driver ', 00687 $ / ' Matrix types (see ZDRVBD for details):', 00688 $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix', 00689 $ / ' 3 = Evenly spaced singular values near 1', 00690 $ / ' 4 = Evenly spaced singular values near underflow', 00691 $ / ' 5 = Evenly spaced singular values near overflow', 00692 $ / / ' Tests performed: ( A is dense, U and V are unitary,', 00693 $ / 19X, ' S is an array, and Upartial, VTpartial, and', 00694 $ / 19X, ' Spartial are partially computed U, VT and S),', / ) 00695 9998 FORMAT( ' Tests performed with Test Threshold = ', F8.2, 00696 $ / ' ZGESVD: ', / 00697 $ ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ', 00698 $ / ' 2 = | I - U**T U | / ( M ulp ) ', 00699 $ / ' 3 = | I - VT VT**T | / ( N ulp ) ', 00700 $ / ' 4 = 0 if S contains min(M,N) nonnegative values in', 00701 $ ' decreasing order, else 1/ulp', 00702 $ / ' 5 = | U - Upartial | / ( M ulp )', 00703 $ / ' 6 = | VT - VTpartial | / ( N ulp )', 00704 $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )', 00705 $ / ' ZGESDD: ', / 00706 $ ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ', 00707 $ / ' 9 = | I - U**T U | / ( M ulp ) ', 00708 $ / '10 = | I - VT VT**T | / ( N ulp ) ', 00709 $ / '11 = 0 if S contains min(M,N) nonnegative values in', 00710 $ ' decreasing order, else 1/ulp', 00711 $ / '12 = | U - Upartial | / ( M ulp )', 00712 $ / '13 = | VT - VTpartial | / ( N ulp )', 00713 $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )', / / ) 00714 9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1, 00715 $ ', seed=', 4( I4, ',' ), ' test(', I1, ')=', G11.4 ) 00716 9996 FORMAT( ' ZDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=', 00717 $ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), 00718 $ I5, ')' ) 00719 9995 FORMAT( ' ZDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=', 00720 $ I6, ', N=', I6, ', JTYPE=', I6, ', LSWORK=', I6, / 9X, 00721 $ 'ISEED=(', 3( I5, ',' ), I5, ')' ) 00722 * 00723 RETURN 00724 * 00725 * End of ZDRVBD 00726 * 00727 END
1.7.3 
