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