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