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