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