 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SDRVLS( 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 REAL THRESH 00013 * .. 00014 * .. Array Arguments .. 00015 LOGICAL DOTYPE( * ) 00016 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ), 00017 $ NVAL( * ), NXVAL( * ) 00018 REAL A( * ), B( * ), C( * ), COPYA( * ), COPYB( * ), 00019 $ COPYS( * ), S( * ), WORK( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * SDRVLS tests the least squares driver routines SGELS, SGELSS, SGELSX, 00026 * SGELSY and SGELSD. 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) REAL 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) REAL 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) REAL array, dimension (MMAX*NMAX) 00090 * 00091 * B (workspace) REAL 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) REAL array, dimension (MMAX*NSMAX) 00096 * 00097 * C (workspace) REAL array, dimension (MMAX*NSMAX) 00098 * 00099 * S (workspace) REAL array, dimension 00100 * (min(MMAX,NMAX)) 00101 * 00102 * COPYS (workspace) REAL array, dimension 00103 * (min(MMAX,NMAX)) 00104 * 00105 * WORK (workspace) REAL 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 REAL ONE, TWO, ZERO 00121 PARAMETER ( ONE = 1.0E0, TWO = 2.0E0, ZERO = 0.0E0 ) 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 REAL EPS, NORMA, NORMB, RCOND 00131 * .. 00132 * .. Local Arrays .. 00133 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00134 REAL RESULT( NTESTS ) 00135 * .. 00136 * .. External Functions .. 00137 REAL SASUM, SLAMCH, SQRT12, SQRT14, SQRT17 00138 EXTERNAL SASUM, SLAMCH, SQRT12, SQRT14, SQRT17 00139 * .. 00140 * .. External Subroutines .. 00141 EXTERNAL ALAERH, ALAHD, ALASVM, SAXPY, SERRLS, SGELS, 00142 $ SGELSD, SGELSS, SGELSX, SGELSY, SGEMM, SLACPY, 00143 $ SLARNV, SQRT13, SQRT15, SQRT16, SSCAL, 00144 $ XLAENV 00145 * .. 00146 * .. Intrinsic Functions .. 00147 INTRINSIC INT, LOG, MAX, MIN, REAL, 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 ) = 'Single 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 = SLAMCH( '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 SERRLS( 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 * 00192 DO 150 IM = 1, NM 00193 M = MVAL( IM ) 00194 LDA = MAX( 1, M ) 00195 * 00196 DO 140 IN = 1, NN 00197 N = NVAL( IN ) 00198 MNMIN = MIN( M, N ) 00199 LDB = MAX( 1, M, N ) 00200 * 00201 DO 130 INS = 1, NNS 00202 NRHS = NSVAL( INS ) 00203 NLVL = MAX( INT( LOG( MAX( ONE, REAL( MNMIN ) ) / 00204 $ REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1, 0 ) 00205 LWORK = MAX( 1, ( M+NRHS )*( N+2 ), ( N+NRHS )*( M+2 ), 00206 $ M*N+4*MNMIN+MAX( M, N ), 12*MNMIN+2*MNMIN*SMLSIZ+ 00207 $ 8*MNMIN*NLVL+MNMIN*NRHS+(SMLSIZ+1)**2 ) 00208 * 00209 DO 120 IRANK = 1, 2 00210 DO 110 ISCALE = 1, 3 00211 ITYPE = ( IRANK-1 )*3 + ISCALE 00212 IF( .NOT.DOTYPE( ITYPE ) ) 00213 $ GO TO 110 00214 * 00215 IF( IRANK.EQ.1 ) THEN 00216 * 00217 * Test SGELS 00218 * 00219 * Generate a matrix of scaling type ISCALE 00220 * 00221 CALL SQRT13( ISCALE, M, N, COPYA, LDA, NORMA, 00222 $ ISEED ) 00223 DO 40 INB = 1, NNB 00224 NB = NBVAL( INB ) 00225 CALL XLAENV( 1, NB ) 00226 CALL XLAENV( 3, NXVAL( INB ) ) 00227 * 00228 DO 30 ITRAN = 1, 2 00229 IF( ITRAN.EQ.1 ) THEN 00230 TRANS = 'N' 00231 NROWS = M 00232 NCOLS = N 00233 ELSE 00234 TRANS = 'T' 00235 NROWS = N 00236 NCOLS = M 00237 END IF 00238 LDWORK = MAX( 1, NCOLS ) 00239 * 00240 * Set up a consistent rhs 00241 * 00242 IF( NCOLS.GT.0 ) THEN 00243 CALL SLARNV( 2, ISEED, NCOLS*NRHS, 00244 $ WORK ) 00245 CALL SSCAL( NCOLS*NRHS, 00246 $ ONE / REAL( NCOLS ), WORK, 00247 $ 1 ) 00248 END IF 00249 CALL SGEMM( TRANS, 'No transpose', NROWS, 00250 $ NRHS, NCOLS, ONE, COPYA, LDA, 00251 $ WORK, LDWORK, ZERO, B, LDB ) 00252 CALL SLACPY( 'Full', NROWS, NRHS, B, LDB, 00253 $ COPYB, LDB ) 00254 * 00255 * Solve LS or overdetermined system 00256 * 00257 IF( M.GT.0 .AND. N.GT.0 ) THEN 00258 CALL SLACPY( 'Full', M, N, COPYA, LDA, 00259 $ A, LDA ) 00260 CALL SLACPY( 'Full', NROWS, NRHS, 00261 $ COPYB, LDB, B, LDB ) 00262 END IF 00263 SRNAMT = 'SGELS ' 00264 CALL SGELS( TRANS, M, N, NRHS, A, LDA, B, 00265 $ LDB, WORK, LWORK, INFO ) 00266 IF( INFO.NE.0 ) 00267 $ CALL ALAERH( PATH, 'SGELS ', INFO, 0, 00268 $ TRANS, M, N, NRHS, -1, NB, 00269 $ ITYPE, NFAIL, NERRS, 00270 $ NOUT ) 00271 * 00272 * Check correctness of results 00273 * 00274 LDWORK = MAX( 1, NROWS ) 00275 IF( NROWS.GT.0 .AND. NRHS.GT.0 ) 00276 $ CALL SLACPY( 'Full', NROWS, NRHS, 00277 $ COPYB, LDB, C, LDB ) 00278 CALL SQRT16( TRANS, M, N, NRHS, COPYA, 00279 $ LDA, B, LDB, C, LDB, WORK, 00280 $ RESULT( 1 ) ) 00281 * 00282 IF( ( ITRAN.EQ.1 .AND. M.GE.N ) .OR. 00283 $ ( ITRAN.EQ.2 .AND. M.LT.N ) ) THEN 00284 * 00285 * Solving LS system 00286 * 00287 RESULT( 2 ) = SQRT17( TRANS, 1, M, N, 00288 $ NRHS, COPYA, LDA, B, LDB, 00289 $ COPYB, LDB, C, WORK, 00290 $ LWORK ) 00291 ELSE 00292 * 00293 * Solving overdetermined system 00294 * 00295 RESULT( 2 ) = SQRT14( TRANS, M, N, 00296 $ NRHS, COPYA, LDA, B, LDB, 00297 $ WORK, LWORK ) 00298 END IF 00299 * 00300 * Print information about the tests that 00301 * did not pass the threshold. 00302 * 00303 DO 20 K = 1, 2 00304 IF( RESULT( K ).GE.THRESH ) THEN 00305 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00306 $ CALL ALAHD( NOUT, PATH ) 00307 WRITE( NOUT, FMT = 9999 )TRANS, M, 00308 $ N, NRHS, NB, ITYPE, K, 00309 $ RESULT( K ) 00310 NFAIL = NFAIL + 1 00311 END IF 00312 20 CONTINUE 00313 NRUN = NRUN + 2 00314 30 CONTINUE 00315 40 CONTINUE 00316 END IF 00317 * 00318 * Generate a matrix of scaling type ISCALE and rank 00319 * type IRANK. 00320 * 00321 CALL SQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA, 00322 $ COPYB, LDB, COPYS, RANK, NORMA, NORMB, 00323 $ ISEED, WORK, LWORK ) 00324 * 00325 * workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M) 00326 * 00327 * Initialize vector IWORK. 00328 * 00329 DO 50 J = 1, N 00330 IWORK( J ) = 0 00331 50 CONTINUE 00332 LDWORK = MAX( 1, M ) 00333 * 00334 * Test SGELSX 00335 * 00336 * SGELSX: Compute the minimum-norm solution X 00337 * to min( norm( A * X - B ) ) using a complete 00338 * orthogonal factorization. 00339 * 00340 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 00341 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB ) 00342 * 00343 SRNAMT = 'SGELSX' 00344 CALL SGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK, 00345 $ RCOND, CRANK, WORK, INFO ) 00346 IF( INFO.NE.0 ) 00347 $ CALL ALAERH( PATH, 'SGELSX', INFO, 0, ' ', M, N, 00348 $ NRHS, -1, NB, ITYPE, NFAIL, NERRS, 00349 $ NOUT ) 00350 * 00351 * workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS ) 00352 * 00353 * Test 3: Compute relative error in svd 00354 * workspace: M*N + 4*MIN(M,N) + MAX(M,N) 00355 * 00356 RESULT( 3 ) = SQRT12( CRANK, CRANK, A, LDA, COPYS, 00357 $ WORK, LWORK ) 00358 * 00359 * Test 4: Compute error in solution 00360 * workspace: M*NRHS + M 00361 * 00362 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 00363 $ LDWORK ) 00364 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA, 00365 $ LDA, B, LDB, WORK, LDWORK, 00366 $ WORK( M*NRHS+1 ), RESULT( 4 ) ) 00367 * 00368 * Test 5: Check norm of r'*A 00369 * workspace: NRHS*(M+N) 00370 * 00371 RESULT( 5 ) = ZERO 00372 IF( M.GT.CRANK ) 00373 $ RESULT( 5 ) = SQRT17( 'No transpose', 1, M, N, 00374 $ NRHS, COPYA, LDA, B, LDB, COPYB, 00375 $ LDB, C, WORK, LWORK ) 00376 * 00377 * Test 6: Check if x is in the rowspace of A 00378 * workspace: (M+NRHS)*(N+2) 00379 * 00380 RESULT( 6 ) = ZERO 00381 * 00382 IF( N.GT.CRANK ) 00383 $ RESULT( 6 ) = SQRT14( 'No transpose', M, N, 00384 $ NRHS, COPYA, LDA, B, LDB, WORK, 00385 $ LWORK ) 00386 * 00387 * Print information about the tests that did not 00388 * pass the threshold. 00389 * 00390 DO 60 K = 3, 6 00391 IF( RESULT( K ).GE.THRESH ) THEN 00392 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00393 $ CALL ALAHD( NOUT, PATH ) 00394 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB, 00395 $ ITYPE, K, RESULT( K ) 00396 NFAIL = NFAIL + 1 00397 END IF 00398 60 CONTINUE 00399 NRUN = NRUN + 4 00400 * 00401 * Loop for testing different block sizes. 00402 * 00403 DO 100 INB = 1, NNB 00404 NB = NBVAL( INB ) 00405 CALL XLAENV( 1, NB ) 00406 CALL XLAENV( 3, NXVAL( INB ) ) 00407 * 00408 * Test SGELSY 00409 * 00410 * SGELSY: Compute the minimum-norm solution X 00411 * to min( norm( A * X - B ) ) 00412 * using the rank-revealing orthogonal 00413 * factorization. 00414 * 00415 * Initialize vector IWORK. 00416 * 00417 DO 70 J = 1, N 00418 IWORK( J ) = 0 00419 70 CONTINUE 00420 * 00421 * Set LWLSY to the adequate value. 00422 * 00423 LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ), 00424 $ 2*MNMIN+NB*NRHS ) 00425 * 00426 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 00427 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B, 00428 $ LDB ) 00429 * 00430 SRNAMT = 'SGELSY' 00431 CALL SGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK, 00432 $ RCOND, CRANK, WORK, LWLSY, INFO ) 00433 IF( INFO.NE.0 ) 00434 $ CALL ALAERH( PATH, 'SGELSY', INFO, 0, ' ', M, 00435 $ N, NRHS, -1, NB, ITYPE, NFAIL, 00436 $ NERRS, NOUT ) 00437 * 00438 * Test 7: Compute relative error in svd 00439 * workspace: M*N + 4*MIN(M,N) + MAX(M,N) 00440 * 00441 RESULT( 7 ) = SQRT12( CRANK, CRANK, A, LDA, 00442 $ COPYS, WORK, LWORK ) 00443 * 00444 * Test 8: Compute error in solution 00445 * workspace: M*NRHS + M 00446 * 00447 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 00448 $ LDWORK ) 00449 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA, 00450 $ LDA, B, LDB, WORK, LDWORK, 00451 $ WORK( M*NRHS+1 ), RESULT( 8 ) ) 00452 * 00453 * Test 9: Check norm of r'*A 00454 * workspace: NRHS*(M+N) 00455 * 00456 RESULT( 9 ) = ZERO 00457 IF( M.GT.CRANK ) 00458 $ RESULT( 9 ) = SQRT17( 'No transpose', 1, M, 00459 $ N, NRHS, COPYA, LDA, B, LDB, 00460 $ COPYB, LDB, C, WORK, LWORK ) 00461 * 00462 * Test 10: Check if x is in the rowspace of A 00463 * workspace: (M+NRHS)*(N+2) 00464 * 00465 RESULT( 10 ) = ZERO 00466 * 00467 IF( N.GT.CRANK ) 00468 $ RESULT( 10 ) = SQRT14( 'No transpose', M, N, 00469 $ NRHS, COPYA, LDA, B, LDB, 00470 $ WORK, LWORK ) 00471 * 00472 * Test SGELSS 00473 * 00474 * SGELSS: Compute the minimum-norm solution X 00475 * to min( norm( A * X - B ) ) 00476 * using the SVD. 00477 * 00478 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 00479 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B, 00480 $ LDB ) 00481 SRNAMT = 'SGELSS' 00482 CALL SGELSS( M, N, NRHS, A, LDA, B, LDB, S, 00483 $ RCOND, CRANK, WORK, LWORK, INFO ) 00484 IF( INFO.NE.0 ) 00485 $ CALL ALAERH( PATH, 'SGELSS', INFO, 0, ' ', M, 00486 $ N, NRHS, -1, NB, ITYPE, NFAIL, 00487 $ NERRS, NOUT ) 00488 * 00489 * workspace used: 3*min(m,n) + 00490 * max(2*min(m,n),nrhs,max(m,n)) 00491 * 00492 * Test 11: Compute relative error in svd 00493 * 00494 IF( RANK.GT.0 ) THEN 00495 CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) 00496 RESULT( 11 ) = SASUM( MNMIN, S, 1 ) / 00497 $ SASUM( MNMIN, COPYS, 1 ) / 00498 $ ( EPS*REAL( MNMIN ) ) 00499 ELSE 00500 RESULT( 11 ) = ZERO 00501 END IF 00502 * 00503 * Test 12: Compute error in solution 00504 * 00505 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 00506 $ LDWORK ) 00507 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA, 00508 $ LDA, B, LDB, WORK, LDWORK, 00509 $ WORK( M*NRHS+1 ), RESULT( 12 ) ) 00510 * 00511 * Test 13: Check norm of r'*A 00512 * 00513 RESULT( 13 ) = ZERO 00514 IF( M.GT.CRANK ) 00515 $ RESULT( 13 ) = SQRT17( 'No transpose', 1, M, 00516 $ N, NRHS, COPYA, LDA, B, LDB, 00517 $ COPYB, LDB, C, WORK, LWORK ) 00518 * 00519 * Test 14: Check if x is in the rowspace of A 00520 * 00521 RESULT( 14 ) = ZERO 00522 IF( N.GT.CRANK ) 00523 $ RESULT( 14 ) = SQRT14( 'No transpose', M, N, 00524 $ NRHS, COPYA, LDA, B, LDB, 00525 $ WORK, LWORK ) 00526 * 00527 * Test SGELSD 00528 * 00529 * SGELSD: Compute the minimum-norm solution X 00530 * to min( norm( A * X - B ) ) using a 00531 * divide and conquer SVD. 00532 * 00533 * Initialize vector IWORK. 00534 * 00535 DO 80 J = 1, N 00536 IWORK( J ) = 0 00537 80 CONTINUE 00538 * 00539 CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA ) 00540 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B, 00541 $ LDB ) 00542 * 00543 SRNAMT = 'SGELSD' 00544 CALL SGELSD( M, N, NRHS, A, LDA, B, LDB, S, 00545 $ RCOND, CRANK, WORK, LWORK, IWORK, 00546 $ INFO ) 00547 IF( INFO.NE.0 ) 00548 $ CALL ALAERH( PATH, 'SGELSD', INFO, 0, ' ', M, 00549 $ N, NRHS, -1, NB, ITYPE, NFAIL, 00550 $ NERRS, NOUT ) 00551 * 00552 * Test 15: Compute relative error in svd 00553 * 00554 IF( RANK.GT.0 ) THEN 00555 CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 ) 00556 RESULT( 15 ) = SASUM( MNMIN, S, 1 ) / 00557 $ SASUM( MNMIN, COPYS, 1 ) / 00558 $ ( EPS*REAL( MNMIN ) ) 00559 ELSE 00560 RESULT( 15 ) = ZERO 00561 END IF 00562 * 00563 * Test 16: Compute error in solution 00564 * 00565 CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK, 00566 $ LDWORK ) 00567 CALL SQRT16( 'No transpose', M, N, NRHS, COPYA, 00568 $ LDA, B, LDB, WORK, LDWORK, 00569 $ WORK( M*NRHS+1 ), RESULT( 16 ) ) 00570 * 00571 * Test 17: Check norm of r'*A 00572 * 00573 RESULT( 17 ) = ZERO 00574 IF( M.GT.CRANK ) 00575 $ RESULT( 17 ) = SQRT17( 'No transpose', 1, M, 00576 $ N, NRHS, COPYA, LDA, B, LDB, 00577 $ COPYB, LDB, C, WORK, LWORK ) 00578 * 00579 * Test 18: Check if x is in the rowspace of A 00580 * 00581 RESULT( 18 ) = ZERO 00582 IF( N.GT.CRANK ) 00583 $ RESULT( 18 ) = SQRT14( 'No transpose', M, N, 00584 $ NRHS, COPYA, LDA, B, LDB, 00585 $ WORK, LWORK ) 00586 * 00587 * Print information about the tests that did not 00588 * pass the threshold. 00589 * 00590 DO 90 K = 7, NTESTS 00591 IF( RESULT( K ).GE.THRESH ) THEN 00592 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00593 $ CALL ALAHD( NOUT, PATH ) 00594 WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB, 00595 $ ITYPE, K, RESULT( K ) 00596 NFAIL = NFAIL + 1 00597 END IF 00598 90 CONTINUE 00599 NRUN = NRUN + 12 00600 * 00601 100 CONTINUE 00602 110 CONTINUE 00603 120 CONTINUE 00604 130 CONTINUE 00605 140 CONTINUE 00606 150 CONTINUE 00607 * 00608 * Print a summary of the results. 00609 * 00610 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00611 * 00612 9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4, 00613 $ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 ) 00614 9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4, 00615 $ ', type', I2, ', test(', I2, ')=', G12.5 ) 00616 RETURN 00617 * 00618 * End of SDRVLS 00619 * 00620 END
1.7.3 
