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