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