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