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