|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DCHKSP( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR, 00002 $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, RWORK, 00003 $ IWORK, NOUT ) 00004 * 00005 * -- LAPACK test routine (version 3.1) -- 00006 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 LOGICAL TSTERR 00011 INTEGER NMAX, NN, NNS, NOUT 00012 DOUBLE PRECISION THRESH 00013 * .. 00014 * .. Array Arguments .. 00015 LOGICAL DOTYPE( * ) 00016 INTEGER IWORK( * ), NSVAL( * ), NVAL( * ) 00017 DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ), 00018 $ RWORK( * ), WORK( * ), X( * ), XACT( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * DCHKSP tests DSPTRF, -TRI, -TRS, -RFS, and -CON 00025 * 00026 * Arguments 00027 * ========= 00028 * 00029 * DOTYPE (input) LOGICAL array, dimension (NTYPES) 00030 * The matrix types to be used for testing. Matrices of type j 00031 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 00032 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 00033 * 00034 * NN (input) INTEGER 00035 * The number of values of N contained in the vector NVAL. 00036 * 00037 * NVAL (input) INTEGER array, dimension (NN) 00038 * The values of the matrix dimension N. 00039 * 00040 * NNS (input) INTEGER 00041 * The number of values of NRHS contained in the vector NSVAL. 00042 * 00043 * NSVAL (input) INTEGER array, dimension (NNS) 00044 * The values of the number of right hand sides NRHS. 00045 * 00046 * THRESH (input) DOUBLE PRECISION 00047 * The threshold value for the test ratios. A result is 00048 * included in the output file if RESULT >= THRESH. To have 00049 * every test ratio printed, use THRESH = 0. 00050 * 00051 * TSTERR (input) LOGICAL 00052 * Flag that indicates whether error exits are to be tested. 00053 * 00054 * NMAX (input) INTEGER 00055 * The maximum value permitted for N, used in dimensioning the 00056 * work arrays. 00057 * 00058 * A (workspace) DOUBLE PRECISION array, dimension 00059 * (NMAX*(NMAX+1)/2) 00060 * 00061 * AFAC (workspace) DOUBLE PRECISION array, dimension 00062 * (NMAX*(NMAX+1)/2) 00063 * 00064 * AINV (workspace) DOUBLE PRECISION array, dimension 00065 * (NMAX*(NMAX+1)/2) 00066 * 00067 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX) 00068 * where NSMAX is the largest entry in NSVAL. 00069 * 00070 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX) 00071 * 00072 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX) 00073 * 00074 * WORK (workspace) DOUBLE PRECISION array, dimension 00075 * (NMAX*max(2,NSMAX)) 00076 * 00077 * RWORK (workspace) DOUBLE PRECISION array, 00078 * dimension (NMAX+2*NSMAX) 00079 * 00080 * IWORK (workspace) INTEGER array, dimension (2*NMAX) 00081 * 00082 * NOUT (input) INTEGER 00083 * The unit number for output. 00084 * 00085 * ===================================================================== 00086 * 00087 * .. Parameters .. 00088 DOUBLE PRECISION ZERO 00089 PARAMETER ( ZERO = 0.0D+0 ) 00090 INTEGER NTYPES 00091 PARAMETER ( NTYPES = 10 ) 00092 INTEGER NTESTS 00093 PARAMETER ( NTESTS = 8 ) 00094 * .. 00095 * .. Local Scalars .. 00096 LOGICAL TRFCON, ZEROT 00097 CHARACTER DIST, PACKIT, TYPE, UPLO, XTYPE 00098 CHARACTER*3 PATH 00099 INTEGER I, I1, I2, IMAT, IN, INFO, IOFF, IRHS, IUPLO, 00100 $ IZERO, J, K, KL, KU, LDA, MODE, N, NERRS, 00101 $ NFAIL, NIMAT, NPP, NRHS, NRUN, NT 00102 DOUBLE PRECISION ANORM, CNDNUM, RCOND, RCONDC 00103 * .. 00104 * .. Local Arrays .. 00105 CHARACTER UPLOS( 2 ) 00106 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00107 DOUBLE PRECISION RESULT( NTESTS ) 00108 * .. 00109 * .. External Functions .. 00110 LOGICAL LSAME 00111 DOUBLE PRECISION DGET06, DLANSP 00112 EXTERNAL LSAME, DGET06, DLANSP 00113 * .. 00114 * .. External Subroutines .. 00115 EXTERNAL ALAERH, ALAHD, ALASUM, DCOPY, DERRSY, DGET04, 00116 $ DLACPY, DLARHS, DLATB4, DLATMS, DPPT02, DPPT03, 00117 $ DPPT05, DSPCON, DSPRFS, DSPT01, DSPTRF, DSPTRI, 00118 $ DSPTRS 00119 * .. 00120 * .. Intrinsic Functions .. 00121 INTRINSIC MAX, MIN 00122 * .. 00123 * .. Scalars in Common .. 00124 LOGICAL LERR, OK 00125 CHARACTER*32 SRNAMT 00126 INTEGER INFOT, NUNIT 00127 * .. 00128 * .. Common blocks .. 00129 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00130 COMMON / SRNAMC / SRNAMT 00131 * .. 00132 * .. Data statements .. 00133 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00134 DATA UPLOS / 'U', 'L' / 00135 * .. 00136 * .. Executable Statements .. 00137 * 00138 * Initialize constants and the random number seed. 00139 * 00140 PATH( 1: 1 ) = 'Double precision' 00141 PATH( 2: 3 ) = 'SP' 00142 NRUN = 0 00143 NFAIL = 0 00144 NERRS = 0 00145 DO 10 I = 1, 4 00146 ISEED( I ) = ISEEDY( I ) 00147 10 CONTINUE 00148 * 00149 * Test the error exits 00150 * 00151 IF( TSTERR ) 00152 $ CALL DERRSY( PATH, NOUT ) 00153 INFOT = 0 00154 * 00155 * Do for each value of N in NVAL 00156 * 00157 DO 170 IN = 1, NN 00158 N = NVAL( IN ) 00159 LDA = MAX( N, 1 ) 00160 XTYPE = 'N' 00161 NIMAT = NTYPES 00162 IF( N.LE.0 ) 00163 $ NIMAT = 1 00164 * 00165 IZERO = 0 00166 DO 160 IMAT = 1, NIMAT 00167 * 00168 * Do the tests only if DOTYPE( IMAT ) is true. 00169 * 00170 IF( .NOT.DOTYPE( IMAT ) ) 00171 $ GO TO 160 00172 * 00173 * Skip types 3, 4, 5, or 6 if the matrix size is too small. 00174 * 00175 ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 00176 IF( ZEROT .AND. N.LT.IMAT-2 ) 00177 $ GO TO 160 00178 * 00179 * Do first for UPLO = 'U', then for UPLO = 'L' 00180 * 00181 DO 150 IUPLO = 1, 2 00182 UPLO = UPLOS( IUPLO ) 00183 IF( LSAME( UPLO, 'U' ) ) THEN 00184 PACKIT = 'C' 00185 ELSE 00186 PACKIT = 'R' 00187 END IF 00188 * 00189 * Set up parameters with DLATB4 and generate a test matrix 00190 * with DLATMS. 00191 * 00192 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00193 $ CNDNUM, DIST ) 00194 * 00195 SRNAMT = 'DLATMS' 00196 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 00197 $ CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK, 00198 $ INFO ) 00199 * 00200 * Check error code from DLATMS. 00201 * 00202 IF( INFO.NE.0 ) THEN 00203 CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1, 00204 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00205 GO TO 150 00206 END IF 00207 * 00208 * For types 3-6, zero one or more rows and columns of 00209 * the matrix to test that INFO is returned correctly. 00210 * 00211 IF( ZEROT ) THEN 00212 IF( IMAT.EQ.3 ) THEN 00213 IZERO = 1 00214 ELSE IF( IMAT.EQ.4 ) THEN 00215 IZERO = N 00216 ELSE 00217 IZERO = N / 2 + 1 00218 END IF 00219 * 00220 IF( IMAT.LT.6 ) THEN 00221 * 00222 * Set row and column IZERO to zero. 00223 * 00224 IF( IUPLO.EQ.1 ) THEN 00225 IOFF = ( IZERO-1 )*IZERO / 2 00226 DO 20 I = 1, IZERO - 1 00227 A( IOFF+I ) = ZERO 00228 20 CONTINUE 00229 IOFF = IOFF + IZERO 00230 DO 30 I = IZERO, N 00231 A( IOFF ) = ZERO 00232 IOFF = IOFF + I 00233 30 CONTINUE 00234 ELSE 00235 IOFF = IZERO 00236 DO 40 I = 1, IZERO - 1 00237 A( IOFF ) = ZERO 00238 IOFF = IOFF + N - I 00239 40 CONTINUE 00240 IOFF = IOFF - IZERO 00241 DO 50 I = IZERO, N 00242 A( IOFF+I ) = ZERO 00243 50 CONTINUE 00244 END IF 00245 ELSE 00246 IOFF = 0 00247 IF( IUPLO.EQ.1 ) THEN 00248 * 00249 * Set the first IZERO rows and columns to zero. 00250 * 00251 DO 70 J = 1, N 00252 I2 = MIN( J, IZERO ) 00253 DO 60 I = 1, I2 00254 A( IOFF+I ) = ZERO 00255 60 CONTINUE 00256 IOFF = IOFF + J 00257 70 CONTINUE 00258 ELSE 00259 * 00260 * Set the last IZERO rows and columns to zero. 00261 * 00262 DO 90 J = 1, N 00263 I1 = MAX( J, IZERO ) 00264 DO 80 I = I1, N 00265 A( IOFF+I ) = ZERO 00266 80 CONTINUE 00267 IOFF = IOFF + N - J 00268 90 CONTINUE 00269 END IF 00270 END IF 00271 ELSE 00272 IZERO = 0 00273 END IF 00274 * 00275 * Compute the L*D*L' or U*D*U' factorization of the matrix. 00276 * 00277 NPP = N*( N+1 ) / 2 00278 CALL DCOPY( NPP, A, 1, AFAC, 1 ) 00279 SRNAMT = 'DSPTRF' 00280 CALL DSPTRF( UPLO, N, AFAC, IWORK, INFO ) 00281 * 00282 * Adjust the expected value of INFO to account for 00283 * pivoting. 00284 * 00285 K = IZERO 00286 IF( K.GT.0 ) THEN 00287 100 CONTINUE 00288 IF( IWORK( K ).LT.0 ) THEN 00289 IF( IWORK( K ).NE.-K ) THEN 00290 K = -IWORK( K ) 00291 GO TO 100 00292 END IF 00293 ELSE IF( IWORK( K ).NE.K ) THEN 00294 K = IWORK( K ) 00295 GO TO 100 00296 END IF 00297 END IF 00298 * 00299 * Check error code from DSPTRF. 00300 * 00301 IF( INFO.NE.K ) 00302 $ CALL ALAERH( PATH, 'DSPTRF', INFO, K, UPLO, N, N, -1, 00303 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00304 IF( INFO.NE.0 ) THEN 00305 TRFCON = .TRUE. 00306 ELSE 00307 TRFCON = .FALSE. 00308 END IF 00309 * 00310 *+ TEST 1 00311 * Reconstruct matrix from factors and compute residual. 00312 * 00313 CALL DSPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA, RWORK, 00314 $ RESULT( 1 ) ) 00315 NT = 1 00316 * 00317 *+ TEST 2 00318 * Form the inverse and compute the residual. 00319 * 00320 IF( .NOT.TRFCON ) THEN 00321 CALL DCOPY( NPP, AFAC, 1, AINV, 1 ) 00322 SRNAMT = 'DSPTRI' 00323 CALL DSPTRI( UPLO, N, AINV, IWORK, WORK, INFO ) 00324 * 00325 * Check error code from DSPTRI. 00326 * 00327 IF( INFO.NE.0 ) 00328 $ CALL ALAERH( PATH, 'DSPTRI', INFO, 0, UPLO, N, N, 00329 $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00330 * 00331 CALL DPPT03( UPLO, N, A, AINV, WORK, LDA, RWORK, 00332 $ RCONDC, RESULT( 2 ) ) 00333 NT = 2 00334 END IF 00335 * 00336 * Print information about the tests that did not pass 00337 * the threshold. 00338 * 00339 DO 110 K = 1, NT 00340 IF( RESULT( K ).GE.THRESH ) THEN 00341 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00342 $ CALL ALAHD( NOUT, PATH ) 00343 WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, K, 00344 $ RESULT( K ) 00345 NFAIL = NFAIL + 1 00346 END IF 00347 110 CONTINUE 00348 NRUN = NRUN + NT 00349 * 00350 * Do only the condition estimate if INFO is not 0. 00351 * 00352 IF( TRFCON ) THEN 00353 RCONDC = ZERO 00354 GO TO 140 00355 END IF 00356 * 00357 DO 130 IRHS = 1, NNS 00358 NRHS = NSVAL( IRHS ) 00359 * 00360 *+ TEST 3 00361 * Solve and compute residual for A * X = B. 00362 * 00363 SRNAMT = 'DLARHS' 00364 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, 00365 $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED, 00366 $ INFO ) 00367 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00368 * 00369 SRNAMT = 'DSPTRS' 00370 CALL DSPTRS( UPLO, N, NRHS, AFAC, IWORK, X, LDA, 00371 $ INFO ) 00372 * 00373 * Check error code from DSPTRS. 00374 * 00375 IF( INFO.NE.0 ) 00376 $ CALL ALAERH( PATH, 'DSPTRS', INFO, 0, UPLO, N, N, 00377 $ -1, -1, NRHS, IMAT, NFAIL, NERRS, 00378 $ NOUT ) 00379 * 00380 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) 00381 CALL DPPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA, 00382 $ RWORK, RESULT( 3 ) ) 00383 * 00384 *+ TEST 4 00385 * Check solution from generated exact solution. 00386 * 00387 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00388 $ RESULT( 4 ) ) 00389 * 00390 *+ TESTS 5, 6, and 7 00391 * Use iterative refinement to improve the solution. 00392 * 00393 SRNAMT = 'DSPRFS' 00394 CALL DSPRFS( UPLO, N, NRHS, A, AFAC, IWORK, B, LDA, X, 00395 $ LDA, RWORK, RWORK( NRHS+1 ), WORK, 00396 $ IWORK( N+1 ), INFO ) 00397 * 00398 * Check error code from DSPRFS. 00399 * 00400 IF( INFO.NE.0 ) 00401 $ CALL ALAERH( PATH, 'DSPRFS', INFO, 0, UPLO, N, N, 00402 $ -1, -1, NRHS, IMAT, NFAIL, NERRS, 00403 $ NOUT ) 00404 * 00405 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00406 $ RESULT( 5 ) ) 00407 CALL DPPT05( UPLO, N, NRHS, A, B, LDA, X, LDA, XACT, 00408 $ LDA, RWORK, RWORK( NRHS+1 ), 00409 $ RESULT( 6 ) ) 00410 * 00411 * Print information about the tests that did not pass 00412 * the threshold. 00413 * 00414 DO 120 K = 3, 7 00415 IF( RESULT( K ).GE.THRESH ) THEN 00416 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00417 $ CALL ALAHD( NOUT, PATH ) 00418 WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, IMAT, 00419 $ K, RESULT( K ) 00420 NFAIL = NFAIL + 1 00421 END IF 00422 120 CONTINUE 00423 NRUN = NRUN + 5 00424 130 CONTINUE 00425 * 00426 *+ TEST 8 00427 * Get an estimate of RCOND = 1/CNDNUM. 00428 * 00429 140 CONTINUE 00430 ANORM = DLANSP( '1', UPLO, N, A, RWORK ) 00431 SRNAMT = 'DSPCON' 00432 CALL DSPCON( UPLO, N, AFAC, IWORK, ANORM, RCOND, WORK, 00433 $ IWORK( N+1 ), INFO ) 00434 * 00435 * Check error code from DSPCON. 00436 * 00437 IF( INFO.NE.0 ) 00438 $ CALL ALAERH( PATH, 'DSPCON', INFO, 0, UPLO, N, N, -1, 00439 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00440 * 00441 RESULT( 8 ) = DGET06( RCOND, RCONDC ) 00442 * 00443 * Print the test ratio if it is .GE. THRESH. 00444 * 00445 IF( RESULT( 8 ).GE.THRESH ) THEN 00446 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00447 $ CALL ALAHD( NOUT, PATH ) 00448 WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, 8, 00449 $ RESULT( 8 ) 00450 NFAIL = NFAIL + 1 00451 END IF 00452 NRUN = NRUN + 1 00453 150 CONTINUE 00454 160 CONTINUE 00455 170 CONTINUE 00456 * 00457 * Print a summary of the results. 00458 * 00459 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00460 * 00461 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', type ', I2, ', test ', 00462 $ I2, ', ratio =', G12.5 ) 00463 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', 00464 $ I2, ', test(', I2, ') =', G12.5 ) 00465 RETURN 00466 * 00467 * End of DCHKSP 00468 * 00469 END
1.7.3 
