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