![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SDRVRFP 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 SDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL, 00012 * + THRESH, A, ASAV, AFAC, AINV, B, 00013 * + BSAV, XACT, X, ARF, ARFINV, 00014 * + S_WORK_SLATMS, S_WORK_SPOT01, S_TEMP_SPOT02, 00015 * + S_TEMP_SPOT03, S_WORK_SLANSY, 00016 * + S_WORK_SPOT02, S_WORK_SPOT03 ) 00017 * 00018 * .. Scalar Arguments .. 00019 * INTEGER NN, NNS, NNT, NOUT 00020 * REAL THRESH 00021 * .. 00022 * .. Array Arguments .. 00023 * INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT ) 00024 * REAL A( * ) 00025 * REAL AINV( * ) 00026 * REAL ASAV( * ) 00027 * REAL B( * ) 00028 * REAL BSAV( * ) 00029 * REAL AFAC( * ) 00030 * REAL ARF( * ) 00031 * REAL ARFINV( * ) 00032 * REAL XACT( * ) 00033 * REAL X( * ) 00034 * REAL S_WORK_SLATMS( * ) 00035 * REAL S_WORK_SPOT01( * ) 00036 * REAL S_TEMP_SPOT02( * ) 00037 * REAL S_TEMP_SPOT03( * ) 00038 * REAL S_WORK_SLANSY( * ) 00039 * REAL S_WORK_SPOT02( * ) 00040 * REAL S_WORK_SPOT03( * ) 00041 * .. 00042 * 00043 * 00044 *> \par Purpose: 00045 * ============= 00046 *> 00047 *> \verbatim 00048 *> 00049 *> SDRVRFP tests the LAPACK RFP routines: 00050 *> SPFTRF, SPFTRS, and SPFTRI. 00051 *> 00052 *> This testing routine follow the same tests as DDRVPO (test for the full 00053 *> format Symmetric Positive Definite solver). 00054 *> 00055 *> The tests are performed in Full Format, convertion back and forth from 00056 *> full format to RFP format are performed using the routines STRTTF and 00057 *> STFTTR. 00058 *> 00059 *> First, a specific matrix A of size N is created. There is nine types of 00060 *> different matrixes possible. 00061 *> 1. Diagonal 6. Random, CNDNUM = sqrt(0.1/EPS) 00062 *> 2. Random, CNDNUM = 2 7. Random, CNDNUM = 0.1/EPS 00063 *> *3. First row and column zero 8. Scaled near underflow 00064 *> *4. Last row and column zero 9. Scaled near overflow 00065 *> *5. Middle row and column zero 00066 *> (* - tests error exits from SPFTRF, no test ratios are computed) 00067 *> A solution XACT of size N-by-NRHS is created and the associated right 00068 *> hand side B as well. Then SPFTRF is called to compute L (or U), the 00069 *> Cholesky factor of A. Then L (or U) is used to solve the linear system 00070 *> of equations AX = B. This gives X. Then L (or U) is used to compute the 00071 *> inverse of A, AINV. The following four tests are then performed: 00072 *> (1) norm( L*L' - A ) / ( N * norm(A) * EPS ) or 00073 *> norm( U'*U - A ) / ( N * norm(A) * EPS ), 00074 *> (2) norm(B - A*X) / ( norm(A) * norm(X) * EPS ), 00075 *> (3) norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ), 00076 *> (4) ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ), 00077 *> where EPS is the machine precision, RCOND the condition number of A, and 00078 *> norm( . ) the 1-norm for (1,2,3) and the inf-norm for (4). 00079 *> Errors occur when INFO parameter is not as expected. Failures occur when 00080 *> a test ratios is greater than THRES. 00081 *> \endverbatim 00082 * 00083 * Arguments: 00084 * ========== 00085 * 00086 *> \param[in] NOUT 00087 *> \verbatim 00088 *> NOUT is INTEGER 00089 *> The unit number for output. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] NN 00093 *> \verbatim 00094 *> NN is INTEGER 00095 *> The number of values of N contained in the vector NVAL. 00096 *> \endverbatim 00097 *> 00098 *> \param[in] NVAL 00099 *> \verbatim 00100 *> NVAL is INTEGER array, dimension (NN) 00101 *> The values of the matrix dimension N. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] NNS 00105 *> \verbatim 00106 *> NNS is INTEGER 00107 *> The number of values of NRHS contained in the vector NSVAL. 00108 *> \endverbatim 00109 *> 00110 *> \param[in] NSVAL 00111 *> \verbatim 00112 *> NSVAL is INTEGER array, dimension (NNS) 00113 *> The values of the number of right-hand sides NRHS. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] NNT 00117 *> \verbatim 00118 *> NNT is INTEGER 00119 *> The number of values of MATRIX TYPE contained in the vector NTVAL. 00120 *> \endverbatim 00121 *> 00122 *> \param[in] NTVAL 00123 *> \verbatim 00124 *> NTVAL is INTEGER array, dimension (NNT) 00125 *> The values of matrix type (between 0 and 9 for PO/PP/PF matrices). 00126 *> \endverbatim 00127 *> 00128 *> \param[in] THRESH 00129 *> \verbatim 00130 *> THRESH is REAL 00131 *> The threshold value for the test ratios. A result is 00132 *> included in the output file if RESULT >= THRESH. To have 00133 *> every test ratio printed, use THRESH = 0. 00134 *> \endverbatim 00135 *> 00136 *> \param[out] A 00137 *> \verbatim 00138 *> A is REAL array, dimension (NMAX*NMAX) 00139 *> \endverbatim 00140 *> 00141 *> \param[out] ASAV 00142 *> \verbatim 00143 *> ASAV is REAL array, dimension (NMAX*NMAX) 00144 *> \endverbatim 00145 *> 00146 *> \param[out] AFAC 00147 *> \verbatim 00148 *> AFAC is REAL array, dimension (NMAX*NMAX) 00149 *> \endverbatim 00150 *> 00151 *> \param[out] AINV 00152 *> \verbatim 00153 *> AINV is REAL array, dimension (NMAX*NMAX) 00154 *> \endverbatim 00155 *> 00156 *> \param[out] B 00157 *> \verbatim 00158 *> B is REAL array, dimension (NMAX*MAXRHS) 00159 *> \endverbatim 00160 *> 00161 *> \param[out] BSAV 00162 *> \verbatim 00163 *> BSAV is REAL array, dimension (NMAX*MAXRHS) 00164 *> \endverbatim 00165 *> 00166 *> \param[out] XACT 00167 *> \verbatim 00168 *> XACT is REAL array, dimension (NMAX*MAXRHS) 00169 *> \endverbatim 00170 *> 00171 *> \param[out] X 00172 *> \verbatim 00173 *> X is REAL array, dimension (NMAX*MAXRHS) 00174 *> \endverbatim 00175 *> 00176 *> \param[out] ARF 00177 *> \verbatim 00178 *> ARF is REAL array, dimension ((NMAX*(NMAX+1))/2) 00179 *> \endverbatim 00180 *> 00181 *> \param[out] ARFINV 00182 *> \verbatim 00183 *> ARFINV is REAL array, dimension ((NMAX*(NMAX+1))/2) 00184 *> \endverbatim 00185 *> 00186 *> \param[out] S_WORK_SLATMS 00187 *> \verbatim 00188 *> S_WORK_SLATMS is REAL array, dimension ( 3*NMAX ) 00189 *> \endverbatim 00190 *> 00191 *> \param[out] S_WORK_SPOT01 00192 *> \verbatim 00193 *> S_WORK_SPOT01 is REAL array, dimension ( NMAX ) 00194 *> \endverbatim 00195 *> 00196 *> \param[out] S_TEMP_SPOT02 00197 *> \verbatim 00198 *> S_TEMP_SPOT02 is REAL array, dimension ( NMAX*MAXRHS ) 00199 *> \endverbatim 00200 *> 00201 *> \param[out] S_TEMP_SPOT03 00202 *> \verbatim 00203 *> S_TEMP_SPOT03 is REAL array, dimension ( NMAX*NMAX ) 00204 *> \endverbatim 00205 *> 00206 *> \param[out] S_WORK_SLATMS 00207 *> \verbatim 00208 *> S_WORK_SLATMS is REAL array, dimension ( NMAX ) 00209 *> \endverbatim 00210 *> 00211 *> \param[out] S_WORK_SLANSY 00212 *> \verbatim 00213 *> S_WORK_SLANSY is REAL array, dimension ( NMAX ) 00214 *> \endverbatim 00215 *> 00216 *> \param[out] S_WORK_SPOT02 00217 *> \verbatim 00218 *> S_WORK_SPOT02 is REAL array, dimension ( NMAX ) 00219 *> \endverbatim 00220 *> 00221 *> \param[out] S_WORK_SPOT03 00222 *> \verbatim 00223 *> S_WORK_SPOT03 is REAL array, dimension ( NMAX ) 00224 *> \endverbatim 00225 * 00226 * Authors: 00227 * ======== 00228 * 00229 *> \author Univ. of Tennessee 00230 *> \author Univ. of California Berkeley 00231 *> \author Univ. of Colorado Denver 00232 *> \author NAG Ltd. 00233 * 00234 *> \date November 2011 00235 * 00236 *> \ingroup single_lin 00237 * 00238 * ===================================================================== 00239 SUBROUTINE SDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL, 00240 + THRESH, A, ASAV, AFAC, AINV, B, 00241 + BSAV, XACT, X, ARF, ARFINV, 00242 + S_WORK_SLATMS, S_WORK_SPOT01, S_TEMP_SPOT02, 00243 + S_TEMP_SPOT03, S_WORK_SLANSY, 00244 + S_WORK_SPOT02, S_WORK_SPOT03 ) 00245 * 00246 * -- LAPACK test routine (version 3.4.0) -- 00247 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00248 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00249 * November 2011 00250 * 00251 * .. Scalar Arguments .. 00252 INTEGER NN, NNS, NNT, NOUT 00253 REAL THRESH 00254 * .. 00255 * .. Array Arguments .. 00256 INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT ) 00257 REAL A( * ) 00258 REAL AINV( * ) 00259 REAL ASAV( * ) 00260 REAL B( * ) 00261 REAL BSAV( * ) 00262 REAL AFAC( * ) 00263 REAL ARF( * ) 00264 REAL ARFINV( * ) 00265 REAL XACT( * ) 00266 REAL X( * ) 00267 REAL S_WORK_SLATMS( * ) 00268 REAL S_WORK_SPOT01( * ) 00269 REAL S_TEMP_SPOT02( * ) 00270 REAL S_TEMP_SPOT03( * ) 00271 REAL S_WORK_SLANSY( * ) 00272 REAL S_WORK_SPOT02( * ) 00273 REAL S_WORK_SPOT03( * ) 00274 * .. 00275 * 00276 * ===================================================================== 00277 * 00278 * .. Parameters .. 00279 REAL ONE, ZERO 00280 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00281 INTEGER NTESTS 00282 PARAMETER ( NTESTS = 4 ) 00283 * .. 00284 * .. Local Scalars .. 00285 LOGICAL ZEROT 00286 INTEGER I, INFO, IUPLO, LDA, LDB, IMAT, NERRS, NFAIL, 00287 + NRHS, NRUN, IZERO, IOFF, K, NT, N, IFORM, IIN, 00288 + IIT, IIS 00289 CHARACTER DIST, CTYPE, UPLO, CFORM 00290 INTEGER KL, KU, MODE 00291 REAL ANORM, AINVNM, CNDNUM, RCONDC 00292 * .. 00293 * .. Local Arrays .. 00294 CHARACTER UPLOS( 2 ), FORMS( 2 ) 00295 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00296 REAL RESULT( NTESTS ) 00297 * .. 00298 * .. External Functions .. 00299 REAL SLANSY 00300 EXTERNAL SLANSY 00301 * .. 00302 * .. External Subroutines .. 00303 EXTERNAL ALADHD, ALAERH, ALASVM, SGET04, STFTTR, SLACPY, 00304 + SLARHS, SLATB4, SLATMS, SPFTRI, SPFTRF, SPFTRS, 00305 + SPOT01, SPOT02, SPOT03, SPOTRI, SPOTRF, STRTTF 00306 * .. 00307 * .. Scalars in Common .. 00308 CHARACTER*32 SRNAMT 00309 * .. 00310 * .. Common blocks .. 00311 COMMON / SRNAMC / SRNAMT 00312 * .. 00313 * .. Data statements .. 00314 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00315 DATA UPLOS / 'U', 'L' / 00316 DATA FORMS / 'N', 'T' / 00317 * .. 00318 * .. Executable Statements .. 00319 * 00320 * Initialize constants and the random number seed. 00321 * 00322 NRUN = 0 00323 NFAIL = 0 00324 NERRS = 0 00325 DO 10 I = 1, 4 00326 ISEED( I ) = ISEEDY( I ) 00327 10 CONTINUE 00328 * 00329 DO 130 IIN = 1, NN 00330 * 00331 N = NVAL( IIN ) 00332 LDA = MAX( N, 1 ) 00333 LDB = MAX( N, 1 ) 00334 * 00335 DO 980 IIS = 1, NNS 00336 * 00337 NRHS = NSVAL( IIS ) 00338 * 00339 DO 120 IIT = 1, NNT 00340 * 00341 IMAT = NTVAL( IIT ) 00342 * 00343 * If N.EQ.0, only consider the first type 00344 * 00345 IF( N.EQ.0 .AND. IIT.GT.1 ) GO TO 120 00346 * 00347 * Skip types 3, 4, or 5 if the matrix size is too small. 00348 * 00349 IF( IMAT.EQ.4 .AND. N.LE.1 ) GO TO 120 00350 IF( IMAT.EQ.5 .AND. N.LE.2 ) GO TO 120 00351 * 00352 * Do first for UPLO = 'U', then for UPLO = 'L' 00353 * 00354 DO 110 IUPLO = 1, 2 00355 UPLO = UPLOS( IUPLO ) 00356 * 00357 * Do first for CFORM = 'N', then for CFORM = 'C' 00358 * 00359 DO 100 IFORM = 1, 2 00360 CFORM = FORMS( IFORM ) 00361 * 00362 * Set up parameters with SLATB4 and generate a test 00363 * matrix with SLATMS. 00364 * 00365 CALL SLATB4( 'SPO', IMAT, N, N, CTYPE, KL, KU, 00366 + ANORM, MODE, CNDNUM, DIST ) 00367 * 00368 SRNAMT = 'SLATMS' 00369 CALL SLATMS( N, N, DIST, ISEED, CTYPE, 00370 + S_WORK_SLATMS, 00371 + MODE, CNDNUM, ANORM, KL, KU, UPLO, A, 00372 + LDA, S_WORK_SLATMS, INFO ) 00373 * 00374 * Check error code from SLATMS. 00375 * 00376 IF( INFO.NE.0 ) THEN 00377 CALL ALAERH( 'SPF', 'SLATMS', INFO, 0, UPLO, N, 00378 + N, -1, -1, -1, IIT, NFAIL, NERRS, 00379 + NOUT ) 00380 GO TO 100 00381 END IF 00382 * 00383 * For types 3-5, zero one row and column of the matrix to 00384 * test that INFO is returned correctly. 00385 * 00386 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5 00387 IF( ZEROT ) THEN 00388 IF( IIT.EQ.3 ) THEN 00389 IZERO = 1 00390 ELSE IF( IIT.EQ.4 ) THEN 00391 IZERO = N 00392 ELSE 00393 IZERO = N / 2 + 1 00394 END IF 00395 IOFF = ( IZERO-1 )*LDA 00396 * 00397 * Set row and column IZERO of A to 0. 00398 * 00399 IF( IUPLO.EQ.1 ) THEN 00400 DO 20 I = 1, IZERO - 1 00401 A( IOFF+I ) = ZERO 00402 20 CONTINUE 00403 IOFF = IOFF + IZERO 00404 DO 30 I = IZERO, N 00405 A( IOFF ) = ZERO 00406 IOFF = IOFF + LDA 00407 30 CONTINUE 00408 ELSE 00409 IOFF = IZERO 00410 DO 40 I = 1, IZERO - 1 00411 A( IOFF ) = ZERO 00412 IOFF = IOFF + LDA 00413 40 CONTINUE 00414 IOFF = IOFF - IZERO 00415 DO 50 I = IZERO, N 00416 A( IOFF+I ) = ZERO 00417 50 CONTINUE 00418 END IF 00419 ELSE 00420 IZERO = 0 00421 END IF 00422 * 00423 * Save a copy of the matrix A in ASAV. 00424 * 00425 CALL SLACPY( UPLO, N, N, A, LDA, ASAV, LDA ) 00426 * 00427 * Compute the condition number of A (RCONDC). 00428 * 00429 IF( ZEROT ) THEN 00430 RCONDC = ZERO 00431 ELSE 00432 * 00433 * Compute the 1-norm of A. 00434 * 00435 ANORM = SLANSY( '1', UPLO, N, A, LDA, 00436 + S_WORK_SLANSY ) 00437 * 00438 * Factor the matrix A. 00439 * 00440 CALL SPOTRF( UPLO, N, A, LDA, INFO ) 00441 * 00442 * Form the inverse of A. 00443 * 00444 CALL SPOTRI( UPLO, N, A, LDA, INFO ) 00445 * 00446 * Compute the 1-norm condition number of A. 00447 * 00448 AINVNM = SLANSY( '1', UPLO, N, A, LDA, 00449 + S_WORK_SLANSY ) 00450 RCONDC = ( ONE / ANORM ) / AINVNM 00451 * 00452 * Restore the matrix A. 00453 * 00454 CALL SLACPY( UPLO, N, N, ASAV, LDA, A, LDA ) 00455 * 00456 END IF 00457 * 00458 * Form an exact solution and set the right hand side. 00459 * 00460 SRNAMT = 'SLARHS' 00461 CALL SLARHS( 'SPO', 'N', UPLO, ' ', N, N, KL, KU, 00462 + NRHS, A, LDA, XACT, LDA, B, LDA, 00463 + ISEED, INFO ) 00464 CALL SLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA ) 00465 * 00466 * Compute the L*L' or U'*U factorization of the 00467 * matrix and solve the system. 00468 * 00469 CALL SLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) 00470 CALL SLACPY( 'Full', N, NRHS, B, LDB, X, LDB ) 00471 * 00472 SRNAMT = 'STRTTF' 00473 CALL STRTTF( CFORM, UPLO, N, AFAC, LDA, ARF, INFO ) 00474 SRNAMT = 'SPFTRF' 00475 CALL SPFTRF( CFORM, UPLO, N, ARF, INFO ) 00476 * 00477 * Check error code from SPFTRF. 00478 * 00479 IF( INFO.NE.IZERO ) THEN 00480 * 00481 * LANGOU: there is a small hick here: IZERO should 00482 * always be INFO however if INFO is ZERO, ALAERH does not 00483 * complain. 00484 * 00485 CALL ALAERH( 'SPF', 'SPFSV ', INFO, IZERO, 00486 + UPLO, N, N, -1, -1, NRHS, IIT, 00487 + NFAIL, NERRS, NOUT ) 00488 GO TO 100 00489 END IF 00490 * 00491 * Skip the tests if INFO is not 0. 00492 * 00493 IF( INFO.NE.0 ) THEN 00494 GO TO 100 00495 END IF 00496 * 00497 SRNAMT = 'SPFTRS' 00498 CALL SPFTRS( CFORM, UPLO, N, NRHS, ARF, X, LDB, 00499 + INFO ) 00500 * 00501 SRNAMT = 'STFTTR' 00502 CALL STFTTR( CFORM, UPLO, N, ARF, AFAC, LDA, INFO ) 00503 * 00504 * Reconstruct matrix from factors and compute 00505 * residual. 00506 * 00507 CALL SLACPY( UPLO, N, N, AFAC, LDA, ASAV, LDA ) 00508 CALL SPOT01( UPLO, N, A, LDA, AFAC, LDA, 00509 + S_WORK_SPOT01, RESULT( 1 ) ) 00510 CALL SLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA ) 00511 * 00512 * Form the inverse and compute the residual. 00513 * 00514 IF(MOD(N,2).EQ.0)THEN 00515 CALL SLACPY( 'A', N+1, N/2, ARF, N+1, ARFINV, 00516 + N+1 ) 00517 ELSE 00518 CALL SLACPY( 'A', N, (N+1)/2, ARF, N, ARFINV, 00519 + N ) 00520 END IF 00521 * 00522 SRNAMT = 'SPFTRI' 00523 CALL SPFTRI( CFORM, UPLO, N, ARFINV , INFO ) 00524 * 00525 SRNAMT = 'STFTTR' 00526 CALL STFTTR( CFORM, UPLO, N, ARFINV, AINV, LDA, 00527 + INFO ) 00528 * 00529 * Check error code from SPFTRI. 00530 * 00531 IF( INFO.NE.0 ) 00532 + CALL ALAERH( 'SPO', 'SPFTRI', INFO, 0, UPLO, N, 00533 + N, -1, -1, -1, IMAT, NFAIL, NERRS, 00534 + NOUT ) 00535 * 00536 CALL SPOT03( UPLO, N, A, LDA, AINV, LDA, 00537 + S_TEMP_SPOT03, LDA, S_WORK_SPOT03, 00538 + RCONDC, RESULT( 2 ) ) 00539 * 00540 * Compute residual of the computed solution. 00541 * 00542 CALL SLACPY( 'Full', N, NRHS, B, LDA, 00543 + S_TEMP_SPOT02, LDA ) 00544 CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA, 00545 + S_TEMP_SPOT02, LDA, S_WORK_SPOT02, 00546 + RESULT( 3 ) ) 00547 * 00548 * Check solution from generated exact solution. 00549 00550 CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00551 + RESULT( 4 ) ) 00552 NT = 4 00553 * 00554 * Print information about the tests that did not 00555 * pass the threshold. 00556 * 00557 DO 60 K = 1, NT 00558 IF( RESULT( K ).GE.THRESH ) THEN 00559 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00560 + CALL ALADHD( NOUT, 'SPF' ) 00561 WRITE( NOUT, FMT = 9999 )'SPFSV ', UPLO, 00562 + N, IIT, K, RESULT( K ) 00563 NFAIL = NFAIL + 1 00564 END IF 00565 60 CONTINUE 00566 NRUN = NRUN + NT 00567 100 CONTINUE 00568 110 CONTINUE 00569 120 CONTINUE 00570 980 CONTINUE 00571 130 CONTINUE 00572 * 00573 * Print a summary of the results. 00574 * 00575 CALL ALASVM( 'SPF', NOUT, NFAIL, NRUN, NERRS ) 00576 * 00577 9999 FORMAT( 1X, A6, ', UPLO=''', A1, ''', N =', I5, ', type ', I1, 00578 + ', test(', I1, ')=', G12.5 ) 00579 * 00580 RETURN 00581 * 00582 * End of SDRVRFP 00583 * 00584 END