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