|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DDRVPP( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, 00002 $ A, AFAC, ASAV, B, BSAV, X, XACT, S, WORK, 00003 $ RWORK, IWORK, NOUT ) 00004 * 00005 * -- LAPACK test routine (version 3.1) -- 00006 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 LOGICAL TSTERR 00011 INTEGER NMAX, NN, NOUT, NRHS 00012 DOUBLE PRECISION THRESH 00013 * .. 00014 * .. Array Arguments .. 00015 LOGICAL DOTYPE( * ) 00016 INTEGER IWORK( * ), NVAL( * ) 00017 DOUBLE PRECISION A( * ), AFAC( * ), ASAV( * ), B( * ), 00018 $ BSAV( * ), RWORK( * ), S( * ), WORK( * ), 00019 $ X( * ), XACT( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * DDRVPP tests the driver routines DPPSV and -SVX. 00026 * 00027 * Arguments 00028 * ========= 00029 * 00030 * DOTYPE (input) LOGICAL array, dimension (NTYPES) 00031 * The matrix types to be used for testing. Matrices of type j 00032 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 00033 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 00034 * 00035 * NN (input) INTEGER 00036 * The number of values of N contained in the vector NVAL. 00037 * 00038 * NVAL (input) INTEGER array, dimension (NN) 00039 * The values of the matrix dimension N. 00040 * 00041 * NRHS (input) INTEGER 00042 * The number of right hand side vectors to be generated for 00043 * each linear system. 00044 * 00045 * THRESH (input) DOUBLE PRECISION 00046 * The threshold value for the test ratios. A result is 00047 * included in the output file if RESULT >= THRESH. To have 00048 * every test ratio printed, use THRESH = 0. 00049 * 00050 * TSTERR (input) LOGICAL 00051 * Flag that indicates whether error exits are to be tested. 00052 * 00053 * NMAX (input) INTEGER 00054 * The maximum value permitted for N, used in dimensioning the 00055 * work arrays. 00056 * 00057 * A (workspace) DOUBLE PRECISION array, dimension 00058 * (NMAX*(NMAX+1)/2) 00059 * 00060 * AFAC (workspace) DOUBLE PRECISION array, dimension 00061 * (NMAX*(NMAX+1)/2) 00062 * 00063 * ASAV (workspace) DOUBLE PRECISION array, dimension 00064 * (NMAX*(NMAX+1)/2) 00065 * 00066 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS) 00067 * 00068 * BSAV (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS) 00069 * 00070 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS) 00071 * 00072 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS) 00073 * 00074 * S (workspace) DOUBLE PRECISION array, dimension (NMAX) 00075 * 00076 * WORK (workspace) DOUBLE PRECISION array, dimension 00077 * (NMAX*max(3,NRHS)) 00078 * 00079 * RWORK (workspace) DOUBLE PRECISION array, dimension (NMAX+2*NRHS) 00080 * 00081 * IWORK (workspace) INTEGER array, dimension (NMAX) 00082 * 00083 * NOUT (input) INTEGER 00084 * The unit number for output. 00085 * 00086 * ===================================================================== 00087 * 00088 * .. Parameters .. 00089 DOUBLE PRECISION ONE, ZERO 00090 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00091 INTEGER NTYPES 00092 PARAMETER ( NTYPES = 9 ) 00093 INTEGER NTESTS 00094 PARAMETER ( NTESTS = 6 ) 00095 * .. 00096 * .. Local Scalars .. 00097 LOGICAL EQUIL, NOFACT, PREFAC, ZEROT 00098 CHARACTER DIST, EQUED, FACT, PACKIT, TYPE, UPLO, XTYPE 00099 CHARACTER*3 PATH 00100 INTEGER I, IEQUED, IFACT, IMAT, IN, INFO, IOFF, IUPLO, 00101 $ IZERO, K, K1, KL, KU, LDA, MODE, N, NERRS, 00102 $ NFACT, NFAIL, NIMAT, NPP, NRUN, NT 00103 DOUBLE PRECISION AINVNM, AMAX, ANORM, CNDNUM, RCOND, RCONDC, 00104 $ ROLDC, SCOND 00105 * .. 00106 * .. Local Arrays .. 00107 CHARACTER EQUEDS( 2 ), FACTS( 3 ), PACKS( 2 ), UPLOS( 2 ) 00108 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00109 DOUBLE PRECISION RESULT( NTESTS ) 00110 * .. 00111 * .. External Functions .. 00112 LOGICAL LSAME 00113 DOUBLE PRECISION DGET06, DLANSP 00114 EXTERNAL LSAME, DGET06, DLANSP 00115 * .. 00116 * .. External Subroutines .. 00117 EXTERNAL ALADHD, ALAERH, ALASVM, DCOPY, DERRVX, DGET04, 00118 $ DLACPY, DLAQSP, DLARHS, DLASET, DLATB4, DLATMS, 00119 $ DPPEQU, DPPSV, DPPSVX, DPPT01, DPPT02, DPPT05, 00120 $ DPPTRF, DPPTRI 00121 * .. 00122 * .. Scalars in Common .. 00123 LOGICAL LERR, OK 00124 CHARACTER*32 SRNAMT 00125 INTEGER INFOT, NUNIT 00126 * .. 00127 * .. Common blocks .. 00128 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00129 COMMON / SRNAMC / SRNAMT 00130 * .. 00131 * .. Intrinsic Functions .. 00132 INTRINSIC MAX 00133 * .. 00134 * .. Data statements .. 00135 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00136 DATA UPLOS / 'U', 'L' / , FACTS / 'F', 'N', 'E' / , 00137 $ PACKS / 'C', 'R' / , EQUEDS / 'N', 'Y' / 00138 * .. 00139 * .. Executable Statements .. 00140 * 00141 * Initialize constants and the random number seed. 00142 * 00143 PATH( 1: 1 ) = 'Double precision' 00144 PATH( 2: 3 ) = 'PP' 00145 NRUN = 0 00146 NFAIL = 0 00147 NERRS = 0 00148 DO 10 I = 1, 4 00149 ISEED( I ) = ISEEDY( I ) 00150 10 CONTINUE 00151 * 00152 * Test the error exits 00153 * 00154 IF( TSTERR ) 00155 $ CALL DERRVX( PATH, NOUT ) 00156 INFOT = 0 00157 * 00158 * Do for each value of N in NVAL 00159 * 00160 DO 140 IN = 1, NN 00161 N = NVAL( IN ) 00162 LDA = MAX( N, 1 ) 00163 NPP = N*( N+1 ) / 2 00164 XTYPE = 'N' 00165 NIMAT = NTYPES 00166 IF( N.LE.0 ) 00167 $ NIMAT = 1 00168 * 00169 DO 130 IMAT = 1, NIMAT 00170 * 00171 * Do the tests only if DOTYPE( IMAT ) is true. 00172 * 00173 IF( .NOT.DOTYPE( IMAT ) ) 00174 $ GO TO 130 00175 * 00176 * Skip types 3, 4, or 5 if the matrix size is too small. 00177 * 00178 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5 00179 IF( ZEROT .AND. N.LT.IMAT-2 ) 00180 $ GO TO 130 00181 * 00182 * Do first for UPLO = 'U', then for UPLO = 'L' 00183 * 00184 DO 120 IUPLO = 1, 2 00185 UPLO = UPLOS( IUPLO ) 00186 PACKIT = PACKS( IUPLO ) 00187 * 00188 * Set up parameters with DLATB4 and generate a test matrix 00189 * with DLATMS. 00190 * 00191 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00192 $ CNDNUM, DIST ) 00193 RCONDC = ONE / CNDNUM 00194 * 00195 SRNAMT = 'DLATMS' 00196 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 00197 $ CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK, 00198 $ INFO ) 00199 * 00200 * Check error code from DLATMS. 00201 * 00202 IF( INFO.NE.0 ) THEN 00203 CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1, 00204 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00205 GO TO 120 00206 END IF 00207 * 00208 * For types 3-5, zero one row and column of the matrix to 00209 * test that INFO is returned correctly. 00210 * 00211 IF( ZEROT ) THEN 00212 IF( IMAT.EQ.3 ) THEN 00213 IZERO = 1 00214 ELSE IF( IMAT.EQ.4 ) THEN 00215 IZERO = N 00216 ELSE 00217 IZERO = N / 2 + 1 00218 END IF 00219 * 00220 * Set row and column IZERO of A to 0. 00221 * 00222 IF( IUPLO.EQ.1 ) THEN 00223 IOFF = ( IZERO-1 )*IZERO / 2 00224 DO 20 I = 1, IZERO - 1 00225 A( IOFF+I ) = ZERO 00226 20 CONTINUE 00227 IOFF = IOFF + IZERO 00228 DO 30 I = IZERO, N 00229 A( IOFF ) = ZERO 00230 IOFF = IOFF + I 00231 30 CONTINUE 00232 ELSE 00233 IOFF = IZERO 00234 DO 40 I = 1, IZERO - 1 00235 A( IOFF ) = ZERO 00236 IOFF = IOFF + N - I 00237 40 CONTINUE 00238 IOFF = IOFF - IZERO 00239 DO 50 I = IZERO, N 00240 A( IOFF+I ) = ZERO 00241 50 CONTINUE 00242 END IF 00243 ELSE 00244 IZERO = 0 00245 END IF 00246 * 00247 * Save a copy of the matrix A in ASAV. 00248 * 00249 CALL DCOPY( NPP, A, 1, ASAV, 1 ) 00250 * 00251 DO 110 IEQUED = 1, 2 00252 EQUED = EQUEDS( IEQUED ) 00253 IF( IEQUED.EQ.1 ) THEN 00254 NFACT = 3 00255 ELSE 00256 NFACT = 1 00257 END IF 00258 * 00259 DO 100 IFACT = 1, NFACT 00260 FACT = FACTS( IFACT ) 00261 PREFAC = LSAME( FACT, 'F' ) 00262 NOFACT = LSAME( FACT, 'N' ) 00263 EQUIL = LSAME( FACT, 'E' ) 00264 * 00265 IF( ZEROT ) THEN 00266 IF( PREFAC ) 00267 $ GO TO 100 00268 RCONDC = ZERO 00269 * 00270 ELSE IF( .NOT.LSAME( FACT, 'N' ) ) THEN 00271 * 00272 * Compute the condition number for comparison with 00273 * the value returned by DPPSVX (FACT = 'N' reuses 00274 * the condition number from the previous iteration 00275 * with FACT = 'F'). 00276 * 00277 CALL DCOPY( NPP, ASAV, 1, AFAC, 1 ) 00278 IF( EQUIL .OR. IEQUED.GT.1 ) THEN 00279 * 00280 * Compute row and column scale factors to 00281 * equilibrate the matrix A. 00282 * 00283 CALL DPPEQU( UPLO, N, AFAC, S, SCOND, AMAX, 00284 $ INFO ) 00285 IF( INFO.EQ.0 .AND. N.GT.0 ) THEN 00286 IF( IEQUED.GT.1 ) 00287 $ SCOND = ZERO 00288 * 00289 * Equilibrate the matrix. 00290 * 00291 CALL DLAQSP( UPLO, N, AFAC, S, SCOND, 00292 $ AMAX, EQUED ) 00293 END IF 00294 END IF 00295 * 00296 * Save the condition number of the 00297 * non-equilibrated system for use in DGET04. 00298 * 00299 IF( EQUIL ) 00300 $ ROLDC = RCONDC 00301 * 00302 * Compute the 1-norm of A. 00303 * 00304 ANORM = DLANSP( '1', UPLO, N, AFAC, RWORK ) 00305 * 00306 * Factor the matrix A. 00307 * 00308 CALL DPPTRF( UPLO, N, AFAC, INFO ) 00309 * 00310 * Form the inverse of A. 00311 * 00312 CALL DCOPY( NPP, AFAC, 1, A, 1 ) 00313 CALL DPPTRI( UPLO, N, A, INFO ) 00314 * 00315 * Compute the 1-norm condition number of A. 00316 * 00317 AINVNM = DLANSP( '1', UPLO, N, A, RWORK ) 00318 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00319 RCONDC = ONE 00320 ELSE 00321 RCONDC = ( ONE / ANORM ) / AINVNM 00322 END IF 00323 END IF 00324 * 00325 * Restore the matrix A. 00326 * 00327 CALL DCOPY( NPP, ASAV, 1, A, 1 ) 00328 * 00329 * Form an exact solution and set the right hand side. 00330 * 00331 SRNAMT = 'DLARHS' 00332 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, 00333 $ NRHS, A, LDA, XACT, LDA, B, LDA, 00334 $ ISEED, INFO ) 00335 XTYPE = 'C' 00336 CALL DLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA ) 00337 * 00338 IF( NOFACT ) THEN 00339 * 00340 * --- Test DPPSV --- 00341 * 00342 * Compute the L*L' or U'*U factorization of the 00343 * matrix and solve the system. 00344 * 00345 CALL DCOPY( NPP, A, 1, AFAC, 1 ) 00346 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00347 * 00348 SRNAMT = 'DPPSV ' 00349 CALL DPPSV( UPLO, N, NRHS, AFAC, X, LDA, INFO ) 00350 * 00351 * Check error code from DPPSV . 00352 * 00353 IF( INFO.NE.IZERO ) THEN 00354 CALL ALAERH( PATH, 'DPPSV ', INFO, IZERO, 00355 $ UPLO, N, N, -1, -1, NRHS, IMAT, 00356 $ NFAIL, NERRS, NOUT ) 00357 GO TO 70 00358 ELSE IF( INFO.NE.0 ) THEN 00359 GO TO 70 00360 END IF 00361 * 00362 * Reconstruct matrix from factors and compute 00363 * residual. 00364 * 00365 CALL DPPT01( UPLO, N, A, AFAC, RWORK, 00366 $ RESULT( 1 ) ) 00367 * 00368 * Compute residual of the computed solution. 00369 * 00370 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, 00371 $ LDA ) 00372 CALL DPPT02( UPLO, N, NRHS, A, X, LDA, WORK, 00373 $ LDA, RWORK, RESULT( 2 ) ) 00374 * 00375 * Check solution from generated exact solution. 00376 * 00377 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00378 $ RESULT( 3 ) ) 00379 NT = 3 00380 * 00381 * Print information about the tests that did not 00382 * pass the threshold. 00383 * 00384 DO 60 K = 1, NT 00385 IF( RESULT( K ).GE.THRESH ) THEN 00386 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00387 $ CALL ALADHD( NOUT, PATH ) 00388 WRITE( NOUT, FMT = 9999 )'DPPSV ', UPLO, 00389 $ N, IMAT, K, RESULT( K ) 00390 NFAIL = NFAIL + 1 00391 END IF 00392 60 CONTINUE 00393 NRUN = NRUN + NT 00394 70 CONTINUE 00395 END IF 00396 * 00397 * --- Test DPPSVX --- 00398 * 00399 IF( .NOT.PREFAC .AND. NPP.GT.0 ) 00400 $ CALL DLASET( 'Full', NPP, 1, ZERO, ZERO, AFAC, 00401 $ NPP ) 00402 CALL DLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA ) 00403 IF( IEQUED.GT.1 .AND. N.GT.0 ) THEN 00404 * 00405 * Equilibrate the matrix if FACT='F' and 00406 * EQUED='Y'. 00407 * 00408 CALL DLAQSP( UPLO, N, A, S, SCOND, AMAX, EQUED ) 00409 END IF 00410 * 00411 * Solve the system and compute the condition number 00412 * and error bounds using DPPSVX. 00413 * 00414 SRNAMT = 'DPPSVX' 00415 CALL DPPSVX( FACT, UPLO, N, NRHS, A, AFAC, EQUED, 00416 $ S, B, LDA, X, LDA, RCOND, RWORK, 00417 $ RWORK( NRHS+1 ), WORK, IWORK, INFO ) 00418 * 00419 * Check the error code from DPPSVX. 00420 * 00421 IF( INFO.NE.IZERO ) THEN 00422 CALL ALAERH( PATH, 'DPPSVX', INFO, IZERO, 00423 $ FACT // UPLO, N, N, -1, -1, NRHS, 00424 $ IMAT, NFAIL, NERRS, NOUT ) 00425 GO TO 90 00426 END IF 00427 * 00428 IF( INFO.EQ.0 ) THEN 00429 IF( .NOT.PREFAC ) THEN 00430 * 00431 * Reconstruct matrix from factors and compute 00432 * residual. 00433 * 00434 CALL DPPT01( UPLO, N, A, AFAC, 00435 $ RWORK( 2*NRHS+1 ), RESULT( 1 ) ) 00436 K1 = 1 00437 ELSE 00438 K1 = 2 00439 END IF 00440 * 00441 * Compute residual of the computed solution. 00442 * 00443 CALL DLACPY( 'Full', N, NRHS, BSAV, LDA, WORK, 00444 $ LDA ) 00445 CALL DPPT02( UPLO, N, NRHS, ASAV, X, LDA, WORK, 00446 $ LDA, RWORK( 2*NRHS+1 ), 00447 $ RESULT( 2 ) ) 00448 * 00449 * Check solution from generated exact solution. 00450 * 00451 IF( NOFACT .OR. ( PREFAC .AND. LSAME( EQUED, 00452 $ 'N' ) ) ) THEN 00453 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, 00454 $ RCONDC, RESULT( 3 ) ) 00455 ELSE 00456 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, 00457 $ ROLDC, RESULT( 3 ) ) 00458 END IF 00459 * 00460 * Check the error bounds from iterative 00461 * refinement. 00462 * 00463 CALL DPPT05( UPLO, N, NRHS, ASAV, B, LDA, X, 00464 $ LDA, XACT, LDA, RWORK, 00465 $ RWORK( NRHS+1 ), RESULT( 4 ) ) 00466 ELSE 00467 K1 = 6 00468 END IF 00469 * 00470 * Compare RCOND from DPPSVX with the computed value 00471 * in RCONDC. 00472 * 00473 RESULT( 6 ) = DGET06( RCOND, RCONDC ) 00474 * 00475 * Print information about the tests that did not pass 00476 * the threshold. 00477 * 00478 DO 80 K = K1, 6 00479 IF( RESULT( K ).GE.THRESH ) THEN 00480 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00481 $ CALL ALADHD( NOUT, PATH ) 00482 IF( PREFAC ) THEN 00483 WRITE( NOUT, FMT = 9997 )'DPPSVX', FACT, 00484 $ UPLO, N, EQUED, IMAT, K, RESULT( K ) 00485 ELSE 00486 WRITE( NOUT, FMT = 9998 )'DPPSVX', FACT, 00487 $ UPLO, N, IMAT, K, RESULT( K ) 00488 END IF 00489 NFAIL = NFAIL + 1 00490 END IF 00491 80 CONTINUE 00492 NRUN = NRUN + 7 - K1 00493 90 CONTINUE 00494 100 CONTINUE 00495 110 CONTINUE 00496 120 CONTINUE 00497 130 CONTINUE 00498 140 CONTINUE 00499 * 00500 * Print a summary of the results. 00501 * 00502 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00503 * 00504 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I1, 00505 $ ', test(', I1, ')=', G12.5 ) 00506 9998 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5, 00507 $ ', type ', I1, ', test(', I1, ')=', G12.5 ) 00508 9997 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5, 00509 $ ', EQUED=''', A1, ''', type ', I1, ', test(', I1, ')=', 00510 $ G12.5 ) 00511 RETURN 00512 * 00513 * End of DDRVPP 00514 * 00515 END
1.7.3 
