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