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