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