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