![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CDRVLS 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 CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB, 00012 * NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B, 00013 * COPYB, C, S, COPYS, WORK, RWORK, IWORK, 00014 * NOUT ) 00015 * 00016 * .. Scalar Arguments .. 00017 * LOGICAL TSTERR 00018 * INTEGER NM, NN, NNB, NNS, NOUT 00019 * REAL THRESH 00020 * .. 00021 * .. Array Arguments .. 00022 * LOGICAL DOTYPE( * ) 00023 * INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ), 00024 * $ NVAL( * ), NXVAL( * ) 00025 * REAL COPYS( * ), RWORK( * ), S( * ) 00026 * COMPLEX A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ), 00027 * $ WORK( * ) 00028 * .. 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> CDRVLS tests the least squares driver routines CGELS, CGELSX, CGELSS, 00037 *> CGELSY and CGELSD. 00038 *> \endverbatim 00039 * 00040 * Arguments: 00041 * ========== 00042 * 00043 *> \param[in] DOTYPE 00044 *> \verbatim 00045 *> DOTYPE is LOGICAL array, dimension (NTYPES) 00046 *> The matrix types to be used for testing. Matrices of type j 00047 *> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 00048 *> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 00049 *> The matrix of type j is generated as follows: 00050 *> j=1: A = U*D*V where U and V are random unitary matrices 00051 *> and D has random entries (> 0.1) taken from a uniform 00052 *> distribution (0,1). A is full rank. 00053 *> j=2: The same of 1, but A is scaled up. 00054 *> j=3: The same of 1, but A is scaled down. 00055 *> j=4: A = U*D*V where U and V are random unitary matrices 00056 *> and D has 3*min(M,N)/4 random entries (> 0.1) taken 00057 *> from a uniform distribution (0,1) and the remaining 00058 *> entries set to 0. A is rank-deficient. 00059 *> j=5: The same of 4, but A is scaled up. 00060 *> j=6: The same of 5, but A is scaled down. 00061 *> \endverbatim 00062 *> 00063 *> \param[in] NM 00064 *> \verbatim 00065 *> NM is INTEGER 00066 *> The number of values of M contained in the vector MVAL. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] MVAL 00070 *> \verbatim 00071 *> MVAL is INTEGER array, dimension (NM) 00072 *> The values of the matrix row dimension M. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] NN 00076 *> \verbatim 00077 *> NN is INTEGER 00078 *> The number of values of N contained in the vector NVAL. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] NVAL 00082 *> \verbatim 00083 *> NVAL is INTEGER array, dimension (NN) 00084 *> The values of the matrix column dimension N. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] NNB 00088 *> \verbatim 00089 *> NNB is INTEGER 00090 *> The number of values of NB and NX contained in the 00091 *> vectors NBVAL and NXVAL. The blocking parameters are used 00092 *> in pairs (NB,NX). 00093 *> \endverbatim 00094 *> 00095 *> \param[in] NBVAL 00096 *> \verbatim 00097 *> NBVAL is INTEGER array, dimension (NNB) 00098 *> The values of the blocksize NB. 00099 *> \endverbatim 00100 *> 00101 *> \param[in] NXVAL 00102 *> \verbatim 00103 *> NXVAL is INTEGER array, dimension (NNB) 00104 *> The values of the crossover point NX. 00105 *> \endverbatim 00106 *> 00107 *> \param[in] NNS 00108 *> \verbatim 00109 *> NNS is INTEGER 00110 *> The number of values of NRHS contained in the vector NSVAL. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] NSVAL 00114 *> \verbatim 00115 *> NSVAL is INTEGER array, dimension (NNS) 00116 *> The values of the number of right hand sides NRHS. 00117 *> \endverbatim 00118 *> 00119 *> \param[in] THRESH 00120 *> \verbatim 00121 *> THRESH is REAL 00122 *> The threshold value for the test ratios. A result is 00123 *> included in the output file if RESULT >= THRESH. To have 00124 *> every test ratio printed, use THRESH = 0. 00125 *> \endverbatim 00126 *> 00127 *> \param[in] TSTERR 00128 *> \verbatim 00129 *> TSTERR is LOGICAL 00130 *> Flag that indicates whether error exits are to be tested. 00131 *> \endverbatim 00132 *> 00133 *> \param[out] A 00134 *> \verbatim 00135 *> A is COMPLEX array, dimension (MMAX*NMAX) 00136 *> where MMAX is the maximum value of M in MVAL and NMAX is the 00137 *> maximum value of N in NVAL. 00138 *> \endverbatim 00139 *> 00140 *> \param[out] COPYA 00141 *> \verbatim 00142 *> COPYA is COMPLEX array, dimension (MMAX*NMAX) 00143 *> \endverbatim 00144 *> 00145 *> \param[out] B 00146 *> \verbatim 00147 *> B is COMPLEX array, dimension (MMAX*NSMAX) 00148 *> where MMAX is the maximum value of M in MVAL and NSMAX is the 00149 *> maximum value of NRHS in NSVAL. 00150 *> \endverbatim 00151 *> 00152 *> \param[out] COPYB 00153 *> \verbatim 00154 *> COPYB is COMPLEX array, dimension (MMAX*NSMAX) 00155 *> \endverbatim 00156 *> 00157 *> \param[out] C 00158 *> \verbatim 00159 *> C is COMPLEX array, dimension (MMAX*NSMAX) 00160 *> \endverbatim 00161 *> 00162 *> \param[out] S 00163 *> \verbatim 00164 *> S is REAL array, dimension 00165 *> (min(MMAX,NMAX)) 00166 *> \endverbatim 00167 *> 00168 *> \param[out] COPYS 00169 *> \verbatim 00170 *> COPYS is REAL array, dimension 00171 *> (min(MMAX,NMAX)) 00172 *> \endverbatim 00173 *> 00174 *> \param[out] WORK 00175 *> \verbatim 00176 *> WORK is COMPLEX array, dimension 00177 *> (MMAX*NMAX + 4*NMAX + MMAX). 00178 *> \endverbatim 00179 *> 00180 *> \param[out] RWORK 00181 *> \verbatim 00182 *> RWORK is REAL array, dimension (5*NMAX-1) 00183 *> \endverbatim 00184 *> 00185 *> \param[out] IWORK 00186 *> \verbatim 00187 *> IWORK is INTEGER array, dimension (15*NMAX) 00188 *> \endverbatim 00189 *> 00190 *> \param[in] NOUT 00191 *> \verbatim 00192 *> NOUT is INTEGER 00193 *> The unit number for output. 00194 *> \endverbatim 00195 * 00196 * Authors: 00197 * ======== 00198 * 00199 *> \author Univ. of Tennessee 00200 *> \author Univ. of California Berkeley 00201 *> \author Univ. of Colorado Denver 00202 *> \author NAG Ltd. 00203 * 00204 *> \date November 2011 00205 * 00206 *> \ingroup complex_lin 00207 * 00208 * ===================================================================== 00209 SUBROUTINE CDRVLS( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL, NNB, 00210 $ NBVAL, NXVAL, THRESH, TSTERR, A, COPYA, B, 00211 $ COPYB, C, S, COPYS, WORK, RWORK, IWORK, 00212 $ NOUT ) 00213 * 00214 * -- LAPACK test routine (version 3.4.0) -- 00215 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00216 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00217 * November 2011 00218 * 00219 * .. Scalar Arguments .. 00220 LOGICAL TSTERR 00221 INTEGER NM, NN, NNB, NNS, NOUT 00222 REAL THRESH 00223 * .. 00224 * .. Array Arguments .. 00225 LOGICAL DOTYPE( * ) 00226 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ), 00227 $ NVAL( * ), NXVAL( * ) 00228 REAL COPYS( * ), RWORK( * ), S( * ) 00229 COMPLEX A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ), 00230 $ WORK( * ) 00231 * .. 00232 * 00233 * ===================================================================== 00234 * 00235 * .. Parameters .. 00236 INTEGER NTESTS 00237 PARAMETER ( NTESTS = 18 ) 00238 INTEGER SMLSIZ 00239 PARAMETER ( SMLSIZ = 25 ) 00240 REAL ONE, ZERO 00241 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00242 COMPLEX CONE, CZERO 00243 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), 00244 $ CZERO = ( 0.0E+0, 0.0E+0 ) ) 00245 * .. 00246 * .. Local Scalars .. 00247 CHARACTER TRANS 00248 CHARACTER*3 PATH 00249 INTEGER CRANK, I, IM, IN, INB, INFO, INS, IRANK, 00250 $ ISCALE, ITRAN, ITYPE, J, K, LDA, LDB, LDWORK, 00251 $ LWLSY, LWORK, M, MNMIN, N, NB, NCOLS, NERRS, 00252 $ NFAIL, NRHS, NROWS, NRUN, RANK 00253 REAL EPS, NORMA, NORMB, RCOND 00254 * .. 00255 * .. Local Arrays .. 00256 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00257 REAL RESULT( NTESTS ) 00258 * .. 00259 * .. External Functions .. 00260 REAL CQRT12, CQRT14, CQRT17, SASUM, SLAMCH 00261 EXTERNAL CQRT12, CQRT14, CQRT17, SASUM, SLAMCH 00262 * .. 00263 * .. External Subroutines .. 00264 EXTERNAL ALAERH, ALAHD, ALASVM, CERRLS, CGELS, CGELSD, 00265 $ CGELSS, CGELSX, CGELSY, CGEMM, CLACPY, CLARNV, 00266 $ CQRT13, CQRT15, CQRT16, CSSCAL, SAXPY, 00267 $ XLAENV 00268 * .. 00269 * .. Intrinsic Functions .. 00270 INTRINSIC MAX, MIN, REAL, SQRT 00271 * .. 00272 * .. Scalars in Common .. 00273 LOGICAL LERR, OK 00274 CHARACTER*32 SRNAMT 00275 INTEGER INFOT, IOUNIT 00276 * .. 00277 * .. Common blocks .. 00278 COMMON / INFOC / INFOT, IOUNIT, OK, LERR 00279 COMMON / SRNAMC / SRNAMT 00280 * .. 00281 * .. Data statements .. 00282 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00283 * .. 00284 * .. Executable Statements .. 00285 * 00286 * Initialize constants and the random number seed. 00287 * 00288 PATH( 1: 1 ) = 'Complex precision' 00289 PATH( 2: 3 ) = 'LS' 00290 NRUN = 0 00291 NFAIL = 0 00292 NERRS = 0 00293 DO 10 I = 1, 4 00294 ISEED( I ) = ISEEDY( I ) 00295 10 CONTINUE 00296 EPS = SLAMCH( 'Epsilon' ) 00297 * 00298 * Threshold for rank estimation 00299 * 00300 RCOND = SQRT( EPS ) - ( SQRT( EPS )-EPS ) / 2 00301 * 00302 * Test the error exits 00303 * 00304 CALL XLAENV( 9, SMLSIZ ) 00305 IF( TSTERR ) 00306 $ CALL CERRLS( PATH, NOUT ) 00307 * 00308 * Print the header if NM = 0 or NN = 0 and THRESH = 0. 00309 * 00310 IF( ( NM.EQ.0 .OR. NN.EQ.0 ) .AND. THRESH.EQ.ZERO ) 00311 $ CALL ALAHD( NOUT, PATH ) 00312 INFOT = 0 00313 * 00314 DO 140 IM = 1, NM 00315 M = MVAL( IM ) 00316 LDA = MAX( 1, M ) 00317 * 00318 DO 130 IN = 1, NN 00319 N = NVAL( IN ) 00320 MNMIN = MIN( M, N ) 00321 LDB = MAX( 1, M, N ) 00322 * 00323 DO 120 INS = 1, NNS 00324 NRHS = NSVAL( INS ) 00325 LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ), 00326 $ M*N+4*MNMIN+MAX( M, N ), 2*N+M ) 00327 * 00328 DO 110 IRANK = 1, 2 00329 DO 100 ISCALE = 1, 3 00330 ITYPE = ( IRANK-1 )*3 + ISCALE 00331 IF( .NOT.DOTYPE( ITYPE ) ) 00332 $ GO TO 100 00333 * 00334 IF( IRANK.EQ.1 ) THEN 00335 * 00336 * Test CGELS 00337 * 00338 * Generate a matrix of scaling type ISCALE 00339 * 00340 CALL CQRT13( ISCALE, M, N, COPYA, LDA, NORMA, 00341 $ ISEED ) 00342 DO 40 INB = 1, NNB 00343 NB = NBVAL( INB ) 00344 CALL XLAENV( 1, NB ) 00345 CALL XLAENV( 3, NXVAL( INB ) ) 00346 * 00347 DO 30 ITRAN = 1, 2 00348 IF( ITRAN.EQ.1 ) THEN 00349 TRANS = 'N' 00350 NROWS = M 00351 NCOLS = N 00352 ELSE 00353 TRANS = 'C' 00354 NROWS = N 00355 NCOLS = M 00356 END IF 00357 LDWORK = MAX( 1, NCOLS ) 00358 * 00359 * Set up a consistent rhs 00360 * 00361 IF( NCOLS.GT.0 ) THEN 00362 CALL CLARNV( 2, ISEED, NCOLS*NRHS, 00363 $ WORK ) 00364 CALL CSSCAL( NCOLS*NRHS, 00365 $ ONE / REAL( NCOLS ), WORK, 00366 $ 1 ) 00367 END IF 00368 CALL CGEMM( TRANS, 'No transpose', NROWS, 00369 $ NRHS, NCOLS, CONE, COPYA, LDA, 00370 $ WORK, LDWORK, CZERO, B, LDB ) 00371 CALL CLACPY( 'Full', NROWS, NRHS, B, LDB, 00372 $ COPYB, LDB ) 00373 * 00374 * Solve LS or overdetermined system 00375 * 00376 IF( M.GT.0 .AND. N.GT.0 ) THEN 00377 CALL CLACPY( 'Full', M, N, COPYA, LDA, 00378 $ A, LDA ) 00379 CALL CLACPY( 'Full', NROWS, NRHS, 00380 $ COPYB, LDB, B, LDB ) 00381 END IF 00382 SRNAMT = 'CGELS ' 00383 CALL CGELS( TRANS, M, N, NRHS, A, LDA, B, 00384 $ LDB, WORK, LWORK, INFO ) 00385 * 00386 IF( INFO.NE.0 ) 00387 $ CALL ALAERH( PATH, 'CGELS ', INFO, 0, 00388 $ TRANS, M, N, NRHS, -1, NB, 00389 $ ITYPE, NFAIL, NERRS, 00390 $ NOUT ) 00391 * 00392 * Check correctness of results 00393 * 00394 LDWORK = MAX( 1, NROWS ) 00395 IF( NROWS.GT.0 .AND. NRHS.GT.0 ) 00396 $ CALL CLACPY( 'Full', NROWS, NRHS, 00397 $ COPYB, LDB, C, LDB ) 00398 CALL CQRT16( TRANS, M, N, NRHS, COPYA, 00399 $ LDA, B, LDB, C, LDB, RWORK, 00400 $ RESULT( 1 ) ) 00401 * 00402 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. 00403 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN 00404 * 00405 * Solving LS system 00406 * 00407 RESULT( 2 ) = CQRT17( TRANS, 1, M, N, 00408 $ NRHS, COPYA, LDA, B, LDB, 00409 $ COPYB, LDB, C, WORK, 00410 $ LWORK ) 00411 ELSE 00412 * 00413 * Solving overdetermined system 00414 * 00415 RESULT( 2 ) = CQRT14( TRANS, M, N, 00416 $ NRHS, COPYA, LDA, B, LDB, 00417 $ WORK, LWORK ) 00418 END IF 00419 * 00420 * Print information about the tests that 00421 * did not pass the threshold. 00422 * 00423 DO 20 K = 1, 2 00424 IF( RESULT( K ).GE.THRESH ) THEN 00425 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00426 $ CALL ALAHD( NOUT, PATH ) 00427 WRITE( NOUT, FMT = 9999 )TRANS, M, 00428 $ N, NRHS, NB, ITYPE, K, 00429 $ RESULT( K ) 00430 NFAIL = NFAIL + 1 00431 END IF 00432 20 CONTINUE 00433 NRUN = NRUN + 2 00434 30 CONTINUE 00435 40 CONTINUE 00436 END IF 00437 * 00438 * Generate a matrix of scaling type ISCALE and rank 00439 * type IRANK. 00440 * 00441 CALL CQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA, 00442 $ COPYB, LDB, COPYS, RANK, NORMA, NORMB, 00443 $ ISEED, WORK, LWORK ) 00444 * 00445 * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M) 00446 * 00447 DO 50 J = 1, N 00448 IWORK( J ) = 0 00449 50 CONTINUE 00450 LDWORK = MAX( 1, M ) 00451 * 00452 * Test CGELSX 00453 * 00454 * CGELSX: Compute the minimum-norm solution X 00455 * to min( norm( A * X - B ) ) 00456 * using a complete orthogonal factorization. 00457 * 00458 CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 00459 CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB ) 00460 * 00461 SRNAMT = 'CGELSX' 00462 CALL CGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK, 00463 $ RCOND, CRANK, WORK, RWORK, INFO ) 00464 * 00465 IF( INFO.NE.0 ) 00466 $ CALL ALAERH( PATH, 'CGELSX', INFO, 0, ' ', M, N, 00467 $ NRHS, -1, NB, ITYPE, NFAIL, NERRS, 00468 $ NOUT ) 00469 * 00470 * workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS ) 00471 * 00472 * Test 3: Compute relative error in svd 00473 * workspace: M*N + 4*MIN(M,N) + MAX(M,N) 00474 * 00475 RESULT( 3 ) = CQRT12( CRANK, CRANK, A, LDA, COPYS, 00476 $ WORK, LWORK, RWORK ) 00477 * 00478 * Test 4: Compute error in solution 00479 * workspace: M*NRHS + M 00480 * 00481 CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 00482 $ LDWORK ) 00483 CALL CQRT16( 'No transpose', M, N, NRHS, COPYA, 00484 $ LDA, B, LDB, WORK, LDWORK, RWORK, 00485 $ RESULT( 4 ) ) 00486 * 00487 * Test 5: Check norm of r'*A 00488 * workspace: NRHS*(M+N) 00489 * 00490 RESULT( 5 ) = ZERO 00491 IF( M.GT.CRANK ) 00492 $ RESULT( 5 ) = CQRT17( 'No transpose', 1, M, N, 00493 $ NRHS, COPYA, LDA, B, LDB, COPYB, 00494 $ LDB, C, WORK, LWORK ) 00495 * 00496 * Test 6: Check if x is in the rowspace of A 00497 * workspace: (M+NRHS)*(N+2) 00498 * 00499 RESULT( 6 ) = ZERO 00500 * 00501 IF( N.GT.CRANK ) 00502 $ RESULT( 6 ) = CQRT14( 'No transpose', M, N, 00503 $ NRHS, COPYA, LDA, B, LDB, WORK, 00504 $ LWORK ) 00505 * 00506 * Print information about the tests that did not 00507 * pass the threshold. 00508 * 00509 DO 60 K = 3, 6 00510 IF( RESULT( K ).GE.THRESH ) THEN 00511 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00512 $ CALL ALAHD( NOUT, PATH ) 00513 WRITE( NOUT, FMT = 9998 )M, N, NRHS, 0, 00514 $ ITYPE, K, RESULT( K ) 00515 NFAIL = NFAIL + 1 00516 END IF 00517 60 CONTINUE 00518 NRUN = NRUN + 4 00519 * 00520 * Loop for testing different block sizes. 00521 * 00522 DO 90 INB = 1, NNB 00523 NB = NBVAL( INB ) 00524 CALL XLAENV( 1, NB ) 00525 CALL XLAENV( 3, NXVAL( INB ) ) 00526 * 00527 * Test CGELSY 00528 * 00529 * CGELSY: Compute the minimum-norm solution 00530 * X to min( norm( A * X - B ) ) 00531 * using the rank-revealing orthogonal 00532 * factorization. 00533 * 00534 CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 00535 CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B, 00536 $ LDB ) 00537 * 00538 * Initialize vector IWORK. 00539 * 00540 DO 70 J = 1, N 00541 IWORK( J ) = 0 00542 70 CONTINUE 00543 * 00544 * Set LWLSY to the adequate value. 00545 * 00546 LWLSY = MNMIN + MAX( 2*MNMIN, NB*( N+1 ), 00547 $ MNMIN+NB*NRHS ) 00548 LWLSY = MAX( 1, LWLSY ) 00549 * 00550 SRNAMT = 'CGELSY' 00551 CALL CGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK, 00552 $ RCOND, CRANK, WORK, LWLSY, RWORK, 00553 $ INFO ) 00554 IF( INFO.NE.0 ) 00555 $ CALL ALAERH( PATH, 'CGELSY', INFO, 0, ' ', M, 00556 $ N, NRHS, -1, NB, ITYPE, NFAIL, 00557 $ NERRS, NOUT ) 00558 * 00559 * workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS) 00560 * 00561 * Test 7: Compute relative error in svd 00562 * workspace: M*N + 4*MIN(M,N) + MAX(M,N) 00563 * 00564 RESULT( 7 ) = CQRT12( CRANK, CRANK, A, LDA, 00565 $ COPYS, WORK, LWORK, RWORK ) 00566 * 00567 * Test 8: Compute error in solution 00568 * workspace: M*NRHS + M 00569 * 00570 CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 00571 $ LDWORK ) 00572 CALL CQRT16( 'No transpose', M, N, NRHS, COPYA, 00573 $ LDA, B, LDB, WORK, LDWORK, RWORK, 00574 $ RESULT( 8 ) ) 00575 * 00576 * Test 9: Check norm of r'*A 00577 * workspace: NRHS*(M+N) 00578 * 00579 RESULT( 9 ) = ZERO 00580 IF( M.GT.CRANK ) 00581 $ RESULT( 9 ) = CQRT17( 'No transpose', 1, M, 00582 $ N, NRHS, COPYA, LDA, B, LDB, 00583 $ COPYB, LDB, C, WORK, LWORK ) 00584 * 00585 * Test 10: Check if x is in the rowspace of A 00586 * workspace: (M+NRHS)*(N+2) 00587 * 00588 RESULT( 10 ) = ZERO 00589 * 00590 IF( N.GT.CRANK ) 00591 $ RESULT( 10 ) = CQRT14( 'No transpose', M, N, 00592 $ NRHS, COPYA, LDA, B, LDB, 00593 $ WORK, LWORK ) 00594 * 00595 * Test CGELSS 00596 * 00597 * CGELSS: Compute the minimum-norm solution 00598 * X to min( norm( A * X - B ) ) 00599 * using the SVD. 00600 * 00601 CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 00602 CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B, 00603 $ LDB ) 00604 SRNAMT = 'CGELSS' 00605 CALL CGELSS( M, N, NRHS, A, LDA, B, LDB, S, 00606 $ RCOND, CRANK, WORK, LWORK, RWORK, 00607 $ INFO ) 00608 * 00609 IF( INFO.NE.0 ) 00610 $ CALL ALAERH( PATH, 'CGELSS', INFO, 0, ' ', M, 00611 $ N, NRHS, -1, NB, ITYPE, NFAIL, 00612 $ NERRS, NOUT ) 00613 * 00614 * workspace used: 3*min(m,n) + 00615 * max(2*min(m,n),nrhs,max(m,n)) 00616 * 00617 * Test 11: Compute relative error in svd 00618 * 00619 IF( RANK.GT.0 ) THEN 00620 CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) 00621 RESULT( 11 ) = SASUM( MNMIN, S, 1 ) / 00622 $ SASUM( MNMIN, COPYS, 1 ) / 00623 $ ( EPS*REAL( MNMIN ) ) 00624 ELSE 00625 RESULT( 11 ) = ZERO 00626 END IF 00627 * 00628 * Test 12: Compute error in solution 00629 * 00630 CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 00631 $ LDWORK ) 00632 CALL CQRT16( 'No transpose', M, N, NRHS, COPYA, 00633 $ LDA, B, LDB, WORK, LDWORK, RWORK, 00634 $ RESULT( 12 ) ) 00635 * 00636 * Test 13: Check norm of r'*A 00637 * 00638 RESULT( 13 ) = ZERO 00639 IF( M.GT.CRANK ) 00640 $ RESULT( 13 ) = CQRT17( 'No transpose', 1, M, 00641 $ N, NRHS, COPYA, LDA, B, LDB, 00642 $ COPYB, LDB, C, WORK, LWORK ) 00643 * 00644 * Test 14: Check if x is in the rowspace of A 00645 * 00646 RESULT( 14 ) = ZERO 00647 IF( N.GT.CRANK ) 00648 $ RESULT( 14 ) = CQRT14( 'No transpose', M, N, 00649 $ NRHS, COPYA, LDA, B, LDB, 00650 $ WORK, LWORK ) 00651 * 00652 * Test CGELSD 00653 * 00654 * CGELSD: Compute the minimum-norm solution X 00655 * to min( norm( A * X - B ) ) using a 00656 * divide and conquer SVD. 00657 * 00658 CALL XLAENV( 9, 25 ) 00659 * 00660 CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 00661 CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B, 00662 $ LDB ) 00663 * 00664 SRNAMT = 'CGELSD' 00665 CALL CGELSD( M, N, NRHS, A, LDA, B, LDB, S, 00666 $ RCOND, CRANK, WORK, LWORK, RWORK, 00667 $ IWORK, INFO ) 00668 IF( INFO.NE.0 ) 00669 $ CALL ALAERH( PATH, 'CGELSD', INFO, 0, ' ', M, 00670 $ N, NRHS, -1, NB, ITYPE, NFAIL, 00671 $ NERRS, NOUT ) 00672 * 00673 * Test 15: Compute relative error in svd 00674 * 00675 IF( RANK.GT.0 ) THEN 00676 CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) 00677 RESULT( 15 ) = SASUM( MNMIN, S, 1 ) / 00678 $ SASUM( MNMIN, COPYS, 1 ) / 00679 $ ( EPS*REAL( MNMIN ) ) 00680 ELSE 00681 RESULT( 15 ) = ZERO 00682 END IF 00683 * 00684 * Test 16: Compute error in solution 00685 * 00686 CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 00687 $ LDWORK ) 00688 CALL CQRT16( 'No transpose', M, N, NRHS, COPYA, 00689 $ LDA, B, LDB, WORK, LDWORK, RWORK, 00690 $ RESULT( 16 ) ) 00691 * 00692 * Test 17: Check norm of r'*A 00693 * 00694 RESULT( 17 ) = ZERO 00695 IF( M.GT.CRANK ) 00696 $ RESULT( 17 ) = CQRT17( 'No transpose', 1, M, 00697 $ N, NRHS, COPYA, LDA, B, LDB, 00698 $ COPYB, LDB, C, WORK, LWORK ) 00699 * 00700 * Test 18: Check if x is in the rowspace of A 00701 * 00702 RESULT( 18 ) = ZERO 00703 IF( N.GT.CRANK ) 00704 $ RESULT( 18 ) = CQRT14( 'No transpose', M, N, 00705 $ NRHS, COPYA, LDA, B, LDB, 00706 $ WORK, LWORK ) 00707 * 00708 * Print information about the tests that did not 00709 * pass the threshold. 00710 * 00711 DO 80 K = 7, NTESTS 00712 IF( RESULT( K ).GE.THRESH ) THEN 00713 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00714 $ CALL ALAHD( NOUT, PATH ) 00715 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB, 00716 $ ITYPE, K, RESULT( K ) 00717 NFAIL = NFAIL + 1 00718 END IF 00719 80 CONTINUE 00720 NRUN = NRUN + 12 00721 * 00722 90 CONTINUE 00723 100 CONTINUE 00724 110 CONTINUE 00725 120 CONTINUE 00726 130 CONTINUE 00727 140 CONTINUE 00728 * 00729 * Print a summary of the results. 00730 * 00731 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00732 * 00733 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4, 00734 $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 ) 00735 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4, 00736 $ ', type', I2, ', test(', I2, ')=', G12.5 ) 00737 RETURN 00738 * 00739 * End of CDRVLS 00740 * 00741 END