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