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