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