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