 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DCHKBD( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS, 00002 $ ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX, 00003 $ Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK, 00004 $ IWORK, NOUT, 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, LDPT, LDQ, LDX, LWORK, NOUT, NRHS, 00012 $ NSIZES, NTYPES 00013 DOUBLE PRECISION THRESH 00014 * .. 00015 * .. Array Arguments .. 00016 LOGICAL DOTYPE( * ) 00017 INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * ) 00018 DOUBLE PRECISION A( LDA, * ), BD( * ), BE( * ), PT( LDPT, * ), 00019 $ Q( LDQ, * ), S1( * ), S2( * ), U( LDPT, * ), 00020 $ VT( LDPT, * ), WORK( * ), X( LDX, * ), 00021 $ Y( LDX, * ), Z( LDX, * ) 00022 * .. 00023 * 00024 * Purpose 00025 * ======= 00026 * 00027 * DCHKBD checks the singular value decomposition (SVD) routines. 00028 * 00029 * DGEBRD reduces a real general m by n matrix A to upper or lower 00030 * bidiagonal form B by an orthogonal transformation: Q' * A * P = B 00031 * (or A = Q * B * P'). The matrix B is upper bidiagonal if m >= n 00032 * and lower bidiagonal if m < n. 00033 * 00034 * DORGBR generates the orthogonal matrices Q and P' from DGEBRD. 00035 * Note that Q and P are not necessarily square. 00036 * 00037 * DBDSQR computes the singular value decomposition of the bidiagonal 00038 * matrix B as B = U S V'. It is called three times to compute 00039 * 1) B = U S1 V', where S1 is the diagonal matrix of singular 00040 * values and the columns of the matrices U and V are the left 00041 * and right singular vectors, respectively, of B. 00042 * 2) Same as 1), but the singular values are stored in S2 and the 00043 * singular vectors are not computed. 00044 * 3) A = (UQ) S (P'V'), the SVD of the original matrix A. 00045 * In addition, DBDSQR has an option to apply the left orthogonal matrix 00046 * U to a matrix X, useful in least squares applications. 00047 * 00048 * DBDSDC computes the singular value decomposition of the bidiagonal 00049 * matrix B as B = U S V' using divide-and-conquer. It is called twice 00050 * to compute 00051 * 1) B = U S1 V', where S1 is the diagonal matrix of singular 00052 * values and the columns of the matrices U and V are the left 00053 * and right singular vectors, respectively, of B. 00054 * 2) Same as 1), but the singular values are stored in S2 and the 00055 * singular vectors are not computed. 00056 * 00057 * For each pair of matrix dimensions (M,N) and each selected matrix 00058 * type, an M by N matrix A and an M by NRHS matrix X are generated. 00059 * The problem dimensions are as follows 00060 * A: M x N 00061 * Q: M x min(M,N) (but M x M if NRHS > 0) 00062 * P: min(M,N) x N 00063 * B: min(M,N) x min(M,N) 00064 * U, V: min(M,N) x min(M,N) 00065 * S1, S2 diagonal, order min(M,N) 00066 * X: M x NRHS 00067 * 00068 * For each generated matrix, 14 tests are performed: 00069 * 00070 * Test DGEBRD and DORGBR 00071 * 00072 * (1) | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P' 00073 * 00074 * (2) | I - Q' Q | / ( M ulp ) 00075 * 00076 * (3) | I - PT PT' | / ( N ulp ) 00077 * 00078 * Test DBDSQR on bidiagonal matrix B 00079 * 00080 * (4) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V' 00081 * 00082 * (5) | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X 00083 * and Z = U' Y. 00084 * (6) | I - U' U | / ( min(M,N) ulp ) 00085 * 00086 * (7) | I - VT VT' | / ( min(M,N) ulp ) 00087 * 00088 * (8) S1 contains min(M,N) nonnegative values in decreasing order. 00089 * (Return 0 if true, 1/ULP if false.) 00090 * 00091 * (9) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without 00092 * computing U and V. 00093 * 00094 * (10) 0 if the true singular values of B are within THRESH of 00095 * those in S1. 2*THRESH if they are not. (Tested using 00096 * DSVDCH) 00097 * 00098 * Test DBDSQR on matrix A 00099 * 00100 * (11) | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp ) 00101 * 00102 * (12) | X - (QU) Z | / ( |X| max(M,k) ulp ) 00103 * 00104 * (13) | I - (QU)'(QU) | / ( M ulp ) 00105 * 00106 * (14) | I - (VT PT) (PT'VT') | / ( N ulp ) 00107 * 00108 * Test DBDSDC on bidiagonal matrix B 00109 * 00110 * (15) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V' 00111 * 00112 * (16) | I - U' U | / ( min(M,N) ulp ) 00113 * 00114 * (17) | I - VT VT' | / ( min(M,N) ulp ) 00115 * 00116 * (18) S1 contains min(M,N) nonnegative values in decreasing order. 00117 * (Return 0 if true, 1/ULP if false.) 00118 * 00119 * (19) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without 00120 * computing U and V. 00121 * The possible matrix types are 00122 * 00123 * (1) The zero matrix. 00124 * (2) The identity matrix. 00125 * 00126 * (3) A diagonal matrix with evenly spaced entries 00127 * 1, ..., ULP and random signs. 00128 * (ULP = (first number larger than 1) - 1 ) 00129 * (4) A diagonal matrix with geometrically spaced entries 00130 * 1, ..., ULP and random signs. 00131 * (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP 00132 * and random signs. 00133 * 00134 * (6) Same as (3), but multiplied by SQRT( overflow threshold ) 00135 * (7) Same as (3), but multiplied by SQRT( underflow threshold ) 00136 * 00137 * (8) A matrix of the form U D V, where U and V are orthogonal and 00138 * D has evenly spaced entries 1, ..., ULP with random signs 00139 * on the diagonal. 00140 * 00141 * (9) A matrix of the form U D V, where U and V are orthogonal and 00142 * D has geometrically spaced entries 1, ..., ULP with random 00143 * signs on the diagonal. 00144 * 00145 * (10) A matrix of the form U D V, where U and V are orthogonal and 00146 * D has "clustered" entries 1, ULP,..., ULP with random 00147 * signs on the diagonal. 00148 * 00149 * (11) Same as (8), but multiplied by SQRT( overflow threshold ) 00150 * (12) Same as (8), but multiplied by SQRT( underflow threshold ) 00151 * 00152 * (13) Rectangular matrix with random entries chosen from (-1,1). 00153 * (14) Same as (13), but multiplied by SQRT( overflow threshold ) 00154 * (15) Same as (13), but multiplied by SQRT( underflow threshold ) 00155 * 00156 * Special case: 00157 * (16) A bidiagonal matrix with random entries chosen from a 00158 * logarithmic distribution on [ulp^2,ulp^(-2)] (I.e., each 00159 * entry is e^x, where x is chosen uniformly on 00160 * [ 2 log(ulp), -2 log(ulp) ] .) For *this* type: 00161 * (a) DGEBRD is not called to reduce it to bidiagonal form. 00162 * (b) the bidiagonal is min(M,N) x min(M,N); if M<N, the 00163 * matrix will be lower bidiagonal, otherwise upper. 00164 * (c) only tests 5--8 and 14 are performed. 00165 * 00166 * A subset of the full set of matrix types may be selected through 00167 * the logical array DOTYPE. 00168 * 00169 * Arguments 00170 * ========== 00171 * 00172 * NSIZES (input) INTEGER 00173 * The number of values of M and N contained in the vectors 00174 * MVAL and NVAL. The matrix sizes are used in pairs (M,N). 00175 * 00176 * MVAL (input) INTEGER array, dimension (NM) 00177 * The values of the matrix row dimension M. 00178 * 00179 * NVAL (input) INTEGER array, dimension (NM) 00180 * The values of the matrix column dimension N. 00181 * 00182 * NTYPES (input) INTEGER 00183 * The number of elements in DOTYPE. If it is zero, DCHKBD 00184 * does nothing. It must be at least zero. If it is MAXTYP+1 00185 * and NSIZES is 1, then an additional type, MAXTYP+1 is 00186 * defined, which is to use whatever matrices are in A and B. 00187 * This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and 00188 * DOTYPE(MAXTYP+1) is .TRUE. . 00189 * 00190 * DOTYPE (input) LOGICAL array, dimension (NTYPES) 00191 * If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix 00192 * of type j will be generated. If NTYPES is smaller than the 00193 * maximum number of types defined (PARAMETER MAXTYP), then 00194 * types NTYPES+1 through MAXTYP will not be generated. If 00195 * NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through 00196 * DOTYPE(NTYPES) will be ignored. 00197 * 00198 * NRHS (input) INTEGER 00199 * The number of columns in the "right-hand side" matrices X, Y, 00200 * and Z, used in testing DBDSQR. If NRHS = 0, then the 00201 * operations on the right-hand side will not be tested. 00202 * NRHS must be at least 0. 00203 * 00204 * ISEED (input/output) INTEGER array, dimension (4) 00205 * On entry ISEED specifies the seed of the random number 00206 * generator. The array elements should be between 0 and 4095; 00207 * if not they will be reduced mod 4096. Also, ISEED(4) must 00208 * be odd. The values of ISEED are changed on exit, and can be 00209 * used in the next call to DCHKBD to continue the same random 00210 * number sequence. 00211 * 00212 * THRESH (input) DOUBLE PRECISION 00213 * The threshold value for the test ratios. A result is 00214 * included in the output file if RESULT >= THRESH. To have 00215 * every test ratio printed, use THRESH = 0. Note that the 00216 * expected value of the test ratios is O(1), so THRESH should 00217 * be a reasonably small multiple of 1, e.g., 10 or 100. 00218 * 00219 * A (workspace) DOUBLE PRECISION array, dimension (LDA,NMAX) 00220 * where NMAX is the maximum value of N in NVAL. 00221 * 00222 * LDA (input) INTEGER 00223 * The leading dimension of the array A. LDA >= max(1,MMAX), 00224 * where MMAX is the maximum value of M in MVAL. 00225 * 00226 * BD (workspace) DOUBLE PRECISION array, dimension 00227 * (max(min(MVAL(j),NVAL(j)))) 00228 * 00229 * BE (workspace) DOUBLE PRECISION array, dimension 00230 * (max(min(MVAL(j),NVAL(j)))) 00231 * 00232 * S1 (workspace) DOUBLE PRECISION array, dimension 00233 * (max(min(MVAL(j),NVAL(j)))) 00234 * 00235 * S2 (workspace) DOUBLE PRECISION array, dimension 00236 * (max(min(MVAL(j),NVAL(j)))) 00237 * 00238 * X (workspace) DOUBLE PRECISION array, dimension (LDX,NRHS) 00239 * 00240 * LDX (input) INTEGER 00241 * The leading dimension of the arrays X, Y, and Z. 00242 * LDX >= max(1,MMAX) 00243 * 00244 * Y (workspace) DOUBLE PRECISION array, dimension (LDX,NRHS) 00245 * 00246 * Z (workspace) DOUBLE PRECISION array, dimension (LDX,NRHS) 00247 * 00248 * Q (workspace) DOUBLE PRECISION array, dimension (LDQ,MMAX) 00249 * 00250 * LDQ (input) INTEGER 00251 * The leading dimension of the array Q. LDQ >= max(1,MMAX). 00252 * 00253 * PT (workspace) DOUBLE PRECISION array, dimension (LDPT,NMAX) 00254 * 00255 * LDPT (input) INTEGER 00256 * The leading dimension of the arrays PT, U, and V. 00257 * LDPT >= max(1, max(min(MVAL(j),NVAL(j)))). 00258 * 00259 * U (workspace) DOUBLE PRECISION array, dimension 00260 * (LDPT,max(min(MVAL(j),NVAL(j)))) 00261 * 00262 * V (workspace) DOUBLE PRECISION array, dimension 00263 * (LDPT,max(min(MVAL(j),NVAL(j)))) 00264 * 00265 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 00266 * 00267 * LWORK (input) INTEGER 00268 * The number of entries in WORK. This must be at least 00269 * 3(M+N) and M(M + max(M,N,k) + 1) + N*min(M,N) for all 00270 * pairs (M,N)=(MM(j),NN(j)) 00271 * 00272 * IWORK (workspace) INTEGER array, dimension at least 8*min(M,N) 00273 * 00274 * NOUT (input) INTEGER 00275 * The FORTRAN unit number for printing out error messages 00276 * (e.g., if a routine returns IINFO not equal to 0.) 00277 * 00278 * INFO (output) INTEGER 00279 * If 0, then everything ran OK. 00280 * -1: NSIZES < 0 00281 * -2: Some MM(j) < 0 00282 * -3: Some NN(j) < 0 00283 * -4: NTYPES < 0 00284 * -6: NRHS < 0 00285 * -8: THRESH < 0 00286 * -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ). 00287 * -17: LDB < 1 or LDB < MMAX. 00288 * -21: LDQ < 1 or LDQ < MMAX. 00289 * -23: LDPT< 1 or LDPT< MNMAX. 00290 * -27: LWORK too small. 00291 * If DLATMR, SLATMS, DGEBRD, DORGBR, or DBDSQR, 00292 * returns an error code, the 00293 * absolute value of it is returned. 00294 * 00295 *----------------------------------------------------------------------- 00296 * 00297 * Some Local Variables and Parameters: 00298 * ---- ----- --------- --- ---------- 00299 * 00300 * ZERO, ONE Real 0 and 1. 00301 * MAXTYP The number of types defined. 00302 * NTEST The number of tests performed, or which can 00303 * be performed so far, for the current matrix. 00304 * MMAX Largest value in NN. 00305 * NMAX Largest value in NN. 00306 * MNMIN min(MM(j), NN(j)) (the dimension of the bidiagonal 00307 * matrix.) 00308 * MNMAX The maximum value of MNMIN for j=1,...,NSIZES. 00309 * NFAIL The number of tests which have exceeded THRESH 00310 * COND, IMODE Values to be passed to the matrix generators. 00311 * ANORM Norm of A; passed to matrix generators. 00312 * 00313 * OVFL, UNFL Overflow and underflow thresholds. 00314 * RTOVFL, RTUNFL Square roots of the previous 2 values. 00315 * ULP, ULPINV Finest relative precision and its inverse. 00316 * 00317 * The following four arrays decode JTYPE: 00318 * KTYPE(j) The general type (1-10) for type "j". 00319 * KMODE(j) The MODE value to be passed to the matrix 00320 * generator for type "j". 00321 * KMAGN(j) The order of magnitude ( O(1), 00322 * O(overflow^(1/2) ), O(underflow^(1/2) ) 00323 * 00324 * ====================================================================== 00325 * 00326 * .. Parameters .. 00327 DOUBLE PRECISION ZERO, ONE, TWO, HALF 00328 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00329 $ HALF = 0.5D0 ) 00330 INTEGER MAXTYP 00331 PARAMETER ( MAXTYP = 16 ) 00332 * .. 00333 * .. Local Scalars .. 00334 LOGICAL BADMM, BADNN, BIDIAG 00335 CHARACTER UPLO 00336 CHARACTER*3 PATH 00337 INTEGER I, IINFO, IMODE, ITYPE, J, JCOL, JSIZE, JTYPE, 00338 $ LOG2UI, M, MINWRK, MMAX, MNMAX, MNMIN, MQ, 00339 $ MTYPES, N, NFAIL, NMAX, NTEST 00340 DOUBLE PRECISION AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL, 00341 $ TEMP1, TEMP2, ULP, ULPINV, UNFL 00342 * .. 00343 * .. Local Arrays .. 00344 INTEGER IDUM( 1 ), IOLDSD( 4 ), KMAGN( MAXTYP ), 00345 $ KMODE( MAXTYP ), KTYPE( MAXTYP ) 00346 DOUBLE PRECISION DUM( 1 ), DUMMA( 1 ), RESULT( 19 ) 00347 * .. 00348 * .. External Functions .. 00349 DOUBLE PRECISION DLAMCH, DLARND 00350 EXTERNAL DLAMCH, DLARND 00351 * .. 00352 * .. External Subroutines .. 00353 EXTERNAL ALASUM, DBDSDC, DBDSQR, DBDT01, DBDT02, DBDT03, 00354 $ DCOPY, DGEBRD, DGEMM, DLABAD, DLACPY, DLAHD2, 00355 $ DLASET, DLATMR, DLATMS, DORGBR, DORT01, XERBLA 00356 * .. 00357 * .. Intrinsic Functions .. 00358 INTRINSIC ABS, EXP, INT, LOG, MAX, MIN, SQRT 00359 * .. 00360 * .. Scalars in Common .. 00361 LOGICAL LERR, OK 00362 CHARACTER*32 SRNAMT 00363 INTEGER INFOT, NUNIT 00364 * .. 00365 * .. Common blocks .. 00366 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00367 COMMON / SRNAMC / SRNAMT 00368 * .. 00369 * .. Data statements .. 00370 DATA KTYPE / 1, 2, 5*4, 5*6, 3*9, 10 / 00371 DATA KMAGN / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 / 00372 DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0, 00373 $ 0, 0, 0 / 00374 * .. 00375 * .. Executable Statements .. 00376 * 00377 * Check for errors 00378 * 00379 INFO = 0 00380 * 00381 BADMM = .FALSE. 00382 BADNN = .FALSE. 00383 MMAX = 1 00384 NMAX = 1 00385 MNMAX = 1 00386 MINWRK = 1 00387 DO 10 J = 1, NSIZES 00388 MMAX = MAX( MMAX, MVAL( J ) ) 00389 IF( MVAL( J ).LT.0 ) 00390 $ BADMM = .TRUE. 00391 NMAX = MAX( NMAX, NVAL( J ) ) 00392 IF( NVAL( J ).LT.0 ) 00393 $ BADNN = .TRUE. 00394 MNMAX = MAX( MNMAX, MIN( MVAL( J ), NVAL( J ) ) ) 00395 MINWRK = MAX( MINWRK, 3*( MVAL( J )+NVAL( J ) ), 00396 $ MVAL( J )*( MVAL( J )+MAX( MVAL( J ), NVAL( J ), 00397 $ NRHS )+1 )+NVAL( J )*MIN( NVAL( J ), MVAL( J ) ) ) 00398 10 CONTINUE 00399 * 00400 * Check for errors 00401 * 00402 IF( NSIZES.LT.0 ) THEN 00403 INFO = -1 00404 ELSE IF( BADMM ) THEN 00405 INFO = -2 00406 ELSE IF( BADNN ) THEN 00407 INFO = -3 00408 ELSE IF( NTYPES.LT.0 ) THEN 00409 INFO = -4 00410 ELSE IF( NRHS.LT.0 ) THEN 00411 INFO = -6 00412 ELSE IF( LDA.LT.MMAX ) THEN 00413 INFO = -11 00414 ELSE IF( LDX.LT.MMAX ) THEN 00415 INFO = -17 00416 ELSE IF( LDQ.LT.MMAX ) THEN 00417 INFO = -21 00418 ELSE IF( LDPT.LT.MNMAX ) THEN 00419 INFO = -23 00420 ELSE IF( MINWRK.GT.LWORK ) THEN 00421 INFO = -27 00422 END IF 00423 * 00424 IF( INFO.NE.0 ) THEN 00425 CALL XERBLA( 'DCHKBD', -INFO ) 00426 RETURN 00427 END IF 00428 * 00429 * Initialize constants 00430 * 00431 PATH( 1: 1 ) = 'Double precision' 00432 PATH( 2: 3 ) = 'BD' 00433 NFAIL = 0 00434 NTEST = 0 00435 UNFL = DLAMCH( 'Safe minimum' ) 00436 OVFL = DLAMCH( 'Overflow' ) 00437 CALL DLABAD( UNFL, OVFL ) 00438 ULP = DLAMCH( 'Precision' ) 00439 ULPINV = ONE / ULP 00440 LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) ) 00441 RTUNFL = SQRT( UNFL ) 00442 RTOVFL = SQRT( OVFL ) 00443 INFOT = 0 00444 * 00445 * Loop over sizes, types 00446 * 00447 DO 200 JSIZE = 1, NSIZES 00448 M = MVAL( JSIZE ) 00449 N = NVAL( JSIZE ) 00450 MNMIN = MIN( M, N ) 00451 AMNINV = ONE / MAX( M, N, 1 ) 00452 * 00453 IF( NSIZES.NE.1 ) THEN 00454 MTYPES = MIN( MAXTYP, NTYPES ) 00455 ELSE 00456 MTYPES = MIN( MAXTYP+1, NTYPES ) 00457 END IF 00458 * 00459 DO 190 JTYPE = 1, MTYPES 00460 IF( .NOT.DOTYPE( JTYPE ) ) 00461 $ GO TO 190 00462 * 00463 DO 20 J = 1, 4 00464 IOLDSD( J ) = ISEED( J ) 00465 20 CONTINUE 00466 * 00467 DO 30 J = 1, 14 00468 RESULT( J ) = -ONE 00469 30 CONTINUE 00470 * 00471 UPLO = ' ' 00472 * 00473 * Compute "A" 00474 * 00475 * Control parameters: 00476 * 00477 * KMAGN KMODE KTYPE 00478 * =1 O(1) clustered 1 zero 00479 * =2 large clustered 2 identity 00480 * =3 small exponential (none) 00481 * =4 arithmetic diagonal, (w/ eigenvalues) 00482 * =5 random symmetric, w/ eigenvalues 00483 * =6 nonsymmetric, w/ singular values 00484 * =7 random diagonal 00485 * =8 random symmetric 00486 * =9 random nonsymmetric 00487 * =10 random bidiagonal (log. distrib.) 00488 * 00489 IF( MTYPES.GT.MAXTYP ) 00490 $ GO TO 100 00491 * 00492 ITYPE = KTYPE( JTYPE ) 00493 IMODE = KMODE( JTYPE ) 00494 * 00495 * Compute norm 00496 * 00497 GO TO ( 40, 50, 60 )KMAGN( JTYPE ) 00498 * 00499 40 CONTINUE 00500 ANORM = ONE 00501 GO TO 70 00502 * 00503 50 CONTINUE 00504 ANORM = ( RTOVFL*ULP )*AMNINV 00505 GO TO 70 00506 * 00507 60 CONTINUE 00508 ANORM = RTUNFL*MAX( M, N )*ULPINV 00509 GO TO 70 00510 * 00511 70 CONTINUE 00512 * 00513 CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 00514 IINFO = 0 00515 COND = ULPINV 00516 * 00517 BIDIAG = .FALSE. 00518 IF( ITYPE.EQ.1 ) THEN 00519 * 00520 * Zero matrix 00521 * 00522 IINFO = 0 00523 * 00524 ELSE IF( ITYPE.EQ.2 ) THEN 00525 * 00526 * Identity 00527 * 00528 DO 80 JCOL = 1, MNMIN 00529 A( JCOL, JCOL ) = ANORM 00530 80 CONTINUE 00531 * 00532 ELSE IF( ITYPE.EQ.4 ) THEN 00533 * 00534 * Diagonal Matrix, [Eigen]values Specified 00535 * 00536 CALL DLATMS( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, IMODE, 00537 $ COND, ANORM, 0, 0, 'N', A, LDA, 00538 $ WORK( MNMIN+1 ), IINFO ) 00539 * 00540 ELSE IF( ITYPE.EQ.5 ) THEN 00541 * 00542 * Symmetric, eigenvalues specified 00543 * 00544 CALL DLATMS( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, IMODE, 00545 $ COND, ANORM, M, N, 'N', A, LDA, 00546 $ WORK( MNMIN+1 ), IINFO ) 00547 * 00548 ELSE IF( ITYPE.EQ.6 ) THEN 00549 * 00550 * Nonsymmetric, singular values specified 00551 * 00552 CALL DLATMS( M, N, 'S', ISEED, 'N', WORK, IMODE, COND, 00553 $ ANORM, M, N, 'N', A, LDA, WORK( MNMIN+1 ), 00554 $ IINFO ) 00555 * 00556 ELSE IF( ITYPE.EQ.7 ) THEN 00557 * 00558 * Diagonal, random entries 00559 * 00560 CALL DLATMR( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, 6, ONE, 00561 $ ONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE, 00562 $ WORK( 2*MNMIN+1 ), 1, ONE, 'N', IWORK, 0, 0, 00563 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00564 * 00565 ELSE IF( ITYPE.EQ.8 ) THEN 00566 * 00567 * Symmetric, random entries 00568 * 00569 CALL DLATMR( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, 6, ONE, 00570 $ ONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE, 00571 $ WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N, 00572 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00573 * 00574 ELSE IF( ITYPE.EQ.9 ) THEN 00575 * 00576 * Nonsymmetric, random entries 00577 * 00578 CALL DLATMR( M, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE, 00579 $ 'T', 'N', WORK( MNMIN+1 ), 1, ONE, 00580 $ WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N, 00581 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) 00582 * 00583 ELSE IF( ITYPE.EQ.10 ) THEN 00584 * 00585 * Bidiagonal, random entries 00586 * 00587 TEMP1 = -TWO*LOG( ULP ) 00588 DO 90 J = 1, MNMIN 00589 BD( J ) = EXP( TEMP1*DLARND( 2, ISEED ) ) 00590 IF( J.LT.MNMIN ) 00591 $ BE( J ) = EXP( TEMP1*DLARND( 2, ISEED ) ) 00592 90 CONTINUE 00593 * 00594 IINFO = 0 00595 BIDIAG = .TRUE. 00596 IF( M.GE.N ) THEN 00597 UPLO = 'U' 00598 ELSE 00599 UPLO = 'L' 00600 END IF 00601 ELSE 00602 IINFO = 1 00603 END IF 00604 * 00605 IF( IINFO.EQ.0 ) THEN 00606 * 00607 * Generate Right-Hand Side 00608 * 00609 IF( BIDIAG ) THEN 00610 CALL DLATMR( MNMIN, NRHS, 'S', ISEED, 'N', WORK, 6, 00611 $ ONE, ONE, 'T', 'N', WORK( MNMIN+1 ), 1, 00612 $ ONE, WORK( 2*MNMIN+1 ), 1, ONE, 'N', 00613 $ IWORK, MNMIN, NRHS, ZERO, ONE, 'NO', Y, 00614 $ LDX, IWORK, IINFO ) 00615 ELSE 00616 CALL DLATMR( M, NRHS, 'S', ISEED, 'N', WORK, 6, ONE, 00617 $ ONE, 'T', 'N', WORK( M+1 ), 1, ONE, 00618 $ WORK( 2*M+1 ), 1, ONE, 'N', IWORK, M, 00619 $ NRHS, ZERO, ONE, 'NO', X, LDX, IWORK, 00620 $ IINFO ) 00621 END IF 00622 END IF 00623 * 00624 * Error Exit 00625 * 00626 IF( IINFO.NE.0 ) THEN 00627 WRITE( NOUT, FMT = 9998 )'Generator', IINFO, M, N, 00628 $ JTYPE, IOLDSD 00629 INFO = ABS( IINFO ) 00630 RETURN 00631 END IF 00632 * 00633 100 CONTINUE 00634 * 00635 * Call DGEBRD and DORGBR to compute B, Q, and P, do tests. 00636 * 00637 IF( .NOT.BIDIAG ) THEN 00638 * 00639 * Compute transformations to reduce A to bidiagonal form: 00640 * B := Q' * A * P. 00641 * 00642 CALL DLACPY( ' ', M, N, A, LDA, Q, LDQ ) 00643 CALL DGEBRD( M, N, Q, LDQ, BD, BE, WORK, WORK( MNMIN+1 ), 00644 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO ) 00645 * 00646 * Check error code from DGEBRD. 00647 * 00648 IF( IINFO.NE.0 ) THEN 00649 WRITE( NOUT, FMT = 9998 )'DGEBRD', IINFO, M, N, 00650 $ JTYPE, IOLDSD 00651 INFO = ABS( IINFO ) 00652 RETURN 00653 END IF 00654 * 00655 CALL DLACPY( ' ', M, N, Q, LDQ, PT, LDPT ) 00656 IF( M.GE.N ) THEN 00657 UPLO = 'U' 00658 ELSE 00659 UPLO = 'L' 00660 END IF 00661 * 00662 * Generate Q 00663 * 00664 MQ = M 00665 IF( NRHS.LE.0 ) 00666 $ MQ = MNMIN 00667 CALL DORGBR( 'Q', M, MQ, N, Q, LDQ, WORK, 00668 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO ) 00669 * 00670 * Check error code from DORGBR. 00671 * 00672 IF( IINFO.NE.0 ) THEN 00673 WRITE( NOUT, FMT = 9998 )'DORGBR(Q)', IINFO, M, N, 00674 $ JTYPE, IOLDSD 00675 INFO = ABS( IINFO ) 00676 RETURN 00677 END IF 00678 * 00679 * Generate P' 00680 * 00681 CALL DORGBR( 'P', MNMIN, N, M, PT, LDPT, WORK( MNMIN+1 ), 00682 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO ) 00683 * 00684 * Check error code from DORGBR. 00685 * 00686 IF( IINFO.NE.0 ) THEN 00687 WRITE( NOUT, FMT = 9998 )'DORGBR(P)', IINFO, M, N, 00688 $ JTYPE, IOLDSD 00689 INFO = ABS( IINFO ) 00690 RETURN 00691 END IF 00692 * 00693 * Apply Q' to an M by NRHS matrix X: Y := Q' * X. 00694 * 00695 CALL DGEMM( 'Transpose', 'No transpose', M, NRHS, M, ONE, 00696 $ Q, LDQ, X, LDX, ZERO, Y, LDX ) 00697 * 00698 * Test 1: Check the decomposition A := Q * B * PT 00699 * 2: Check the orthogonality of Q 00700 * 3: Check the orthogonality of PT 00701 * 00702 CALL DBDT01( M, N, 1, A, LDA, Q, LDQ, BD, BE, PT, LDPT, 00703 $ WORK, RESULT( 1 ) ) 00704 CALL DORT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK, 00705 $ RESULT( 2 ) ) 00706 CALL DORT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK, 00707 $ RESULT( 3 ) ) 00708 END IF 00709 * 00710 * Use DBDSQR to form the SVD of the bidiagonal matrix B: 00711 * B := U * S1 * VT, and compute Z = U' * Y. 00712 * 00713 CALL DCOPY( MNMIN, BD, 1, S1, 1 ) 00714 IF( MNMIN.GT.0 ) 00715 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 00716 CALL DLACPY( ' ', M, NRHS, Y, LDX, Z, LDX ) 00717 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, U, LDPT ) 00718 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, VT, LDPT ) 00719 * 00720 CALL DBDSQR( UPLO, MNMIN, MNMIN, MNMIN, NRHS, S1, WORK, VT, 00721 $ LDPT, U, LDPT, Z, LDX, WORK( MNMIN+1 ), IINFO ) 00722 * 00723 * Check error code from DBDSQR. 00724 * 00725 IF( IINFO.NE.0 ) THEN 00726 WRITE( NOUT, FMT = 9998 )'DBDSQR(vects)', IINFO, M, N, 00727 $ JTYPE, IOLDSD 00728 INFO = ABS( IINFO ) 00729 IF( IINFO.LT.0 ) THEN 00730 RETURN 00731 ELSE 00732 RESULT( 4 ) = ULPINV 00733 GO TO 170 00734 END IF 00735 END IF 00736 * 00737 * Use DBDSQR to compute only the singular values of the 00738 * bidiagonal matrix B; U, VT, and Z should not be modified. 00739 * 00740 CALL DCOPY( MNMIN, BD, 1, S2, 1 ) 00741 IF( MNMIN.GT.0 ) 00742 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 00743 * 00744 CALL DBDSQR( UPLO, MNMIN, 0, 0, 0, S2, WORK, VT, LDPT, U, 00745 $ LDPT, Z, LDX, WORK( MNMIN+1 ), IINFO ) 00746 * 00747 * Check error code from DBDSQR. 00748 * 00749 IF( IINFO.NE.0 ) THEN 00750 WRITE( NOUT, FMT = 9998 )'DBDSQR(values)', IINFO, M, N, 00751 $ JTYPE, IOLDSD 00752 INFO = ABS( IINFO ) 00753 IF( IINFO.LT.0 ) THEN 00754 RETURN 00755 ELSE 00756 RESULT( 9 ) = ULPINV 00757 GO TO 170 00758 END IF 00759 END IF 00760 * 00761 * Test 4: Check the decomposition B := U * S1 * VT 00762 * 5: Check the computation Z := U' * Y 00763 * 6: Check the orthogonality of U 00764 * 7: Check the orthogonality of VT 00765 * 00766 CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT, 00767 $ WORK, RESULT( 4 ) ) 00768 CALL DBDT02( MNMIN, NRHS, Y, LDX, Z, LDX, U, LDPT, WORK, 00769 $ RESULT( 5 ) ) 00770 CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK, 00771 $ RESULT( 6 ) ) 00772 CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK, 00773 $ RESULT( 7 ) ) 00774 * 00775 * Test 8: Check that the singular values are sorted in 00776 * non-increasing order and are non-negative 00777 * 00778 RESULT( 8 ) = ZERO 00779 DO 110 I = 1, MNMIN - 1 00780 IF( S1( I ).LT.S1( I+1 ) ) 00781 $ RESULT( 8 ) = ULPINV 00782 IF( S1( I ).LT.ZERO ) 00783 $ RESULT( 8 ) = ULPINV 00784 110 CONTINUE 00785 IF( MNMIN.GE.1 ) THEN 00786 IF( S1( MNMIN ).LT.ZERO ) 00787 $ RESULT( 8 ) = ULPINV 00788 END IF 00789 * 00790 * Test 9: Compare DBDSQR with and without singular vectors 00791 * 00792 TEMP2 = ZERO 00793 * 00794 DO 120 J = 1, MNMIN 00795 TEMP1 = ABS( S1( J )-S2( J ) ) / 00796 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ), 00797 $ ULP*MAX( ABS( S1( J ) ), ABS( S2( J ) ) ) ) 00798 TEMP2 = MAX( TEMP1, TEMP2 ) 00799 120 CONTINUE 00800 * 00801 RESULT( 9 ) = TEMP2 00802 * 00803 * Test 10: Sturm sequence test of singular values 00804 * Go up by factors of two until it succeeds 00805 * 00806 TEMP1 = THRESH*( HALF-ULP ) 00807 * 00808 DO 130 J = 0, LOG2UI 00809 * CALL DSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO ) 00810 IF( IINFO.EQ.0 ) 00811 $ GO TO 140 00812 TEMP1 = TEMP1*TWO 00813 130 CONTINUE 00814 * 00815 140 CONTINUE 00816 RESULT( 10 ) = TEMP1 00817 * 00818 * Use DBDSQR to form the decomposition A := (QU) S (VT PT) 00819 * from the bidiagonal form A := Q B PT. 00820 * 00821 IF( .NOT.BIDIAG ) THEN 00822 CALL DCOPY( MNMIN, BD, 1, S2, 1 ) 00823 IF( MNMIN.GT.0 ) 00824 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 00825 * 00826 CALL DBDSQR( UPLO, MNMIN, N, M, NRHS, S2, WORK, PT, LDPT, 00827 $ Q, LDQ, Y, LDX, WORK( MNMIN+1 ), IINFO ) 00828 * 00829 * Test 11: Check the decomposition A := Q*U * S2 * VT*PT 00830 * 12: Check the computation Z := U' * Q' * X 00831 * 13: Check the orthogonality of Q*U 00832 * 14: Check the orthogonality of VT*PT 00833 * 00834 CALL DBDT01( M, N, 0, A, LDA, Q, LDQ, S2, DUMMA, PT, 00835 $ LDPT, WORK, RESULT( 11 ) ) 00836 CALL DBDT02( M, NRHS, X, LDX, Y, LDX, Q, LDQ, WORK, 00837 $ RESULT( 12 ) ) 00838 CALL DORT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK, 00839 $ RESULT( 13 ) ) 00840 CALL DORT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK, 00841 $ RESULT( 14 ) ) 00842 END IF 00843 * 00844 * Use DBDSDC to form the SVD of the bidiagonal matrix B: 00845 * B := U * S1 * VT 00846 * 00847 CALL DCOPY( MNMIN, BD, 1, S1, 1 ) 00848 IF( MNMIN.GT.0 ) 00849 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 00850 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, U, LDPT ) 00851 CALL DLASET( 'Full', MNMIN, MNMIN, ZERO, ONE, VT, LDPT ) 00852 * 00853 CALL DBDSDC( UPLO, 'I', MNMIN, S1, WORK, U, LDPT, VT, LDPT, 00854 $ DUM, IDUM, WORK( MNMIN+1 ), IWORK, IINFO ) 00855 * 00856 * Check error code from DBDSDC. 00857 * 00858 IF( IINFO.NE.0 ) THEN 00859 WRITE( NOUT, FMT = 9998 )'DBDSDC(vects)', IINFO, M, N, 00860 $ JTYPE, IOLDSD 00861 INFO = ABS( IINFO ) 00862 IF( IINFO.LT.0 ) THEN 00863 RETURN 00864 ELSE 00865 RESULT( 15 ) = ULPINV 00866 GO TO 170 00867 END IF 00868 END IF 00869 * 00870 * Use DBDSDC to compute only the singular values of the 00871 * bidiagonal matrix B; U and VT should not be modified. 00872 * 00873 CALL DCOPY( MNMIN, BD, 1, S2, 1 ) 00874 IF( MNMIN.GT.0 ) 00875 $ CALL DCOPY( MNMIN-1, BE, 1, WORK, 1 ) 00876 * 00877 CALL DBDSDC( UPLO, 'N', MNMIN, S2, WORK, DUM, 1, DUM, 1, 00878 $ DUM, IDUM, WORK( MNMIN+1 ), IWORK, IINFO ) 00879 * 00880 * Check error code from DBDSDC. 00881 * 00882 IF( IINFO.NE.0 ) THEN 00883 WRITE( NOUT, FMT = 9998 )'DBDSDC(values)', IINFO, M, N, 00884 $ JTYPE, IOLDSD 00885 INFO = ABS( IINFO ) 00886 IF( IINFO.LT.0 ) THEN 00887 RETURN 00888 ELSE 00889 RESULT( 18 ) = ULPINV 00890 GO TO 170 00891 END IF 00892 END IF 00893 * 00894 * Test 15: Check the decomposition B := U * S1 * VT 00895 * 16: Check the orthogonality of U 00896 * 17: Check the orthogonality of VT 00897 * 00898 CALL DBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT, 00899 $ WORK, RESULT( 15 ) ) 00900 CALL DORT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK, 00901 $ RESULT( 16 ) ) 00902 CALL DORT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK, 00903 $ RESULT( 17 ) ) 00904 * 00905 * Test 18: Check that the singular values are sorted in 00906 * non-increasing order and are non-negative 00907 * 00908 RESULT( 18 ) = ZERO 00909 DO 150 I = 1, MNMIN - 1 00910 IF( S1( I ).LT.S1( I+1 ) ) 00911 $ RESULT( 18 ) = ULPINV 00912 IF( S1( I ).LT.ZERO ) 00913 $ RESULT( 18 ) = ULPINV 00914 150 CONTINUE 00915 IF( MNMIN.GE.1 ) THEN 00916 IF( S1( MNMIN ).LT.ZERO ) 00917 $ RESULT( 18 ) = ULPINV 00918 END IF 00919 * 00920 * Test 19: Compare DBDSQR with and without singular vectors 00921 * 00922 TEMP2 = ZERO 00923 * 00924 DO 160 J = 1, MNMIN 00925 TEMP1 = ABS( S1( J )-S2( J ) ) / 00926 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ), 00927 $ ULP*MAX( ABS( S1( 1 ) ), ABS( S2( 1 ) ) ) ) 00928 TEMP2 = MAX( TEMP1, TEMP2 ) 00929 160 CONTINUE 00930 * 00931 RESULT( 19 ) = TEMP2 00932 * 00933 * End of Loop -- Check for RESULT(j) > THRESH 00934 * 00935 170 CONTINUE 00936 DO 180 J = 1, 19 00937 IF( RESULT( J ).GE.THRESH ) THEN 00938 IF( NFAIL.EQ.0 ) 00939 $ CALL DLAHD2( NOUT, PATH ) 00940 WRITE( NOUT, FMT = 9999 )M, N, JTYPE, IOLDSD, J, 00941 $ RESULT( J ) 00942 NFAIL = NFAIL + 1 00943 END IF 00944 180 CONTINUE 00945 IF( .NOT.BIDIAG ) THEN 00946 NTEST = NTEST + 19 00947 ELSE 00948 NTEST = NTEST + 5 00949 END IF 00950 * 00951 190 CONTINUE 00952 200 CONTINUE 00953 * 00954 * Summary 00955 * 00956 CALL ALASUM( PATH, NOUT, NFAIL, NTEST, 0 ) 00957 * 00958 RETURN 00959 * 00960 * End of DCHKBD 00961 * 00962 9999 FORMAT( ' M=', I5, ', N=', I5, ', type ', I2, ', seed=', 00963 $ 4( I4, ',' ), ' test(', I2, ')=', G11.4 ) 00964 9998 FORMAT( ' DCHKBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=', 00965 $ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), 00966 $ I5, ')' ) 00967 * 00968 END
1.7.3 
