 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DCHKPO( 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.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, 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 * DCHKPO tests DPOTRF, -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 * 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 (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 = 9 ) 00095 INTEGER NTESTS 00096 PARAMETER ( NTESTS = 8 ) 00097 * .. 00098 * .. Local Scalars .. 00099 LOGICAL ZEROT 00100 CHARACTER DIST, TYPE, UPLO, XTYPE 00101 CHARACTER*3 PATH 00102 INTEGER I, IMAT, IN, INB, INFO, IOFF, IRHS, IUPLO, 00103 $ IZERO, K, KL, KU, LDA, MODE, N, NB, NERRS, 00104 $ NFAIL, NIMAT, NRHS, NRUN 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 * .. 00112 * .. External Functions .. 00113 DOUBLE PRECISION DGET06, DLANSY 00114 EXTERNAL DGET06, DLANSY 00115 * .. 00116 * .. External Subroutines .. 00117 EXTERNAL ALAERH, ALAHD, ALASUM, DERRPO, DGET04, DLACPY, 00118 $ DLARHS, DLATB4, DLATMS, DPOCON, DPORFS, DPOT01, 00119 $ DPOT02, DPOT03, DPOT05, DPOTRF, DPOTRI, DPOTRS, 00120 $ XLAENV 00121 * .. 00122 * .. Scalars in Common .. 00123 LOGICAL LERR, OK 00124 CHARACTER*32 SRNAMT 00125 INTEGER INFOT, NUNIT 00126 * .. 00127 * .. Common blocks .. 00128 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00129 COMMON / SRNAMC / SRNAMT 00130 * .. 00131 * .. Intrinsic Functions .. 00132 INTRINSIC MAX 00133 * .. 00134 * .. Data statements .. 00135 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00136 DATA UPLOS / 'U', 'L' / 00137 * .. 00138 * .. Executable Statements .. 00139 * 00140 * Initialize constants and the random number seed. 00141 * 00142 PATH( 1: 1 ) = 'Double precision' 00143 PATH( 2: 3 ) = 'PO' 00144 NRUN = 0 00145 NFAIL = 0 00146 NERRS = 0 00147 DO 10 I = 1, 4 00148 ISEED( I ) = ISEEDY( I ) 00149 10 CONTINUE 00150 * 00151 * Test the error exits 00152 * 00153 IF( TSTERR ) 00154 $ CALL DERRPO( PATH, NOUT ) 00155 INFOT = 0 00156 CALL XLAENV( 2, 2 ) 00157 * 00158 * Do for each value of N in NVAL 00159 * 00160 DO 120 IN = 1, NN 00161 N = NVAL( IN ) 00162 LDA = MAX( N, 1 ) 00163 XTYPE = 'N' 00164 NIMAT = NTYPES 00165 IF( N.LE.0 ) 00166 $ NIMAT = 1 00167 * 00168 IZERO = 0 00169 DO 110 IMAT = 1, NIMAT 00170 * 00171 * Do the tests only if DOTYPE( IMAT ) is true. 00172 * 00173 IF( .NOT.DOTYPE( IMAT ) ) 00174 $ GO TO 110 00175 * 00176 * Skip types 3, 4, or 5 if the matrix size is too small. 00177 * 00178 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5 00179 IF( ZEROT .AND. N.LT.IMAT-2 ) 00180 $ GO TO 110 00181 * 00182 * Do first for UPLO = 'U', then for UPLO = 'L' 00183 * 00184 DO 100 IUPLO = 1, 2 00185 UPLO = UPLOS( IUPLO ) 00186 * 00187 * Set up parameters with DLATB4 and generate a test matrix 00188 * with DLATMS. 00189 * 00190 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00191 $ CNDNUM, DIST ) 00192 * 00193 SRNAMT = 'DLATMS' 00194 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 00195 $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, 00196 $ INFO ) 00197 * 00198 * Check error code from DLATMS. 00199 * 00200 IF( INFO.NE.0 ) THEN 00201 CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1, 00202 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00203 GO TO 100 00204 END IF 00205 * 00206 * For types 3-5, zero one row and column of the matrix to 00207 * test that INFO is returned correctly. 00208 * 00209 IF( ZEROT ) THEN 00210 IF( IMAT.EQ.3 ) THEN 00211 IZERO = 1 00212 ELSE IF( IMAT.EQ.4 ) THEN 00213 IZERO = N 00214 ELSE 00215 IZERO = N / 2 + 1 00216 END IF 00217 IOFF = ( IZERO-1 )*LDA 00218 * 00219 * Set row and column IZERO of A to 0. 00220 * 00221 IF( IUPLO.EQ.1 ) THEN 00222 DO 20 I = 1, IZERO - 1 00223 A( IOFF+I ) = ZERO 00224 20 CONTINUE 00225 IOFF = IOFF + IZERO 00226 DO 30 I = IZERO, N 00227 A( IOFF ) = ZERO 00228 IOFF = IOFF + LDA 00229 30 CONTINUE 00230 ELSE 00231 IOFF = IZERO 00232 DO 40 I = 1, IZERO - 1 00233 A( IOFF ) = ZERO 00234 IOFF = IOFF + LDA 00235 40 CONTINUE 00236 IOFF = IOFF - IZERO 00237 DO 50 I = IZERO, N 00238 A( IOFF+I ) = ZERO 00239 50 CONTINUE 00240 END IF 00241 ELSE 00242 IZERO = 0 00243 END IF 00244 * 00245 * Do for each value of NB in NBVAL 00246 * 00247 DO 90 INB = 1, NNB 00248 NB = NBVAL( INB ) 00249 CALL XLAENV( 1, NB ) 00250 * 00251 * Compute the L*L' or U'*U factorization of the matrix. 00252 * 00253 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) 00254 SRNAMT = 'DPOTRF' 00255 CALL DPOTRF( UPLO, N, AFAC, LDA, INFO ) 00256 * 00257 * Check error code from DPOTRF. 00258 * 00259 IF( INFO.NE.IZERO ) THEN 00260 CALL ALAERH( PATH, 'DPOTRF', INFO, IZERO, UPLO, N, 00261 $ N, -1, -1, NB, IMAT, NFAIL, NERRS, 00262 $ NOUT ) 00263 GO TO 90 00264 END IF 00265 * 00266 * Skip the tests if INFO is not 0. 00267 * 00268 IF( INFO.NE.0 ) 00269 $ GO TO 90 00270 * 00271 *+ TEST 1 00272 * Reconstruct matrix from factors and compute residual. 00273 * 00274 CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) 00275 CALL DPOT01( UPLO, N, A, LDA, AINV, LDA, RWORK, 00276 $ RESULT( 1 ) ) 00277 * 00278 *+ TEST 2 00279 * Form the inverse and compute the residual. 00280 * 00281 CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) 00282 SRNAMT = 'DPOTRI' 00283 CALL DPOTRI( UPLO, N, AINV, LDA, INFO ) 00284 * 00285 * Check error code from DPOTRI. 00286 * 00287 IF( INFO.NE.0 ) 00288 $ CALL ALAERH( PATH, 'DPOTRI', INFO, 0, UPLO, N, N, 00289 $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00290 * 00291 CALL DPOT03( UPLO, N, A, LDA, AINV, LDA, WORK, LDA, 00292 $ RWORK, RCONDC, RESULT( 2 ) ) 00293 * 00294 * Print information about the tests that did not pass 00295 * the threshold. 00296 * 00297 DO 60 K = 1, 2 00298 IF( RESULT( K ).GE.THRESH ) THEN 00299 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00300 $ CALL ALAHD( NOUT, PATH ) 00301 WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K, 00302 $ RESULT( K ) 00303 NFAIL = NFAIL + 1 00304 END IF 00305 60 CONTINUE 00306 NRUN = NRUN + 2 00307 * 00308 * Skip the rest of the tests unless this is the first 00309 * blocksize. 00310 * 00311 IF( INB.NE.1 ) 00312 $ GO TO 90 00313 * 00314 DO 80 IRHS = 1, NNS 00315 NRHS = NSVAL( IRHS ) 00316 * 00317 *+ TEST 3 00318 * Solve and compute residual for A * X = B . 00319 * 00320 SRNAMT = 'DLARHS' 00321 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, 00322 $ NRHS, A, LDA, XACT, LDA, B, LDA, 00323 $ ISEED, INFO ) 00324 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00325 * 00326 SRNAMT = 'DPOTRS' 00327 CALL DPOTRS( UPLO, N, NRHS, AFAC, LDA, X, LDA, 00328 $ INFO ) 00329 * 00330 * Check error code from DPOTRS. 00331 * 00332 IF( INFO.NE.0 ) 00333 $ CALL ALAERH( PATH, 'DPOTRS', INFO, 0, UPLO, N, 00334 $ N, -1, -1, NRHS, IMAT, NFAIL, 00335 $ NERRS, NOUT ) 00336 * 00337 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) 00338 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, 00339 $ LDA, RWORK, RESULT( 3 ) ) 00340 * 00341 *+ TEST 4 00342 * Check solution from generated exact solution. 00343 * 00344 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00345 $ RESULT( 4 ) ) 00346 * 00347 *+ TESTS 5, 6, and 7 00348 * Use iterative refinement to improve the solution. 00349 * 00350 SRNAMT = 'DPORFS' 00351 CALL DPORFS( UPLO, N, NRHS, A, LDA, AFAC, LDA, B, 00352 $ LDA, X, LDA, RWORK, RWORK( NRHS+1 ), 00353 $ WORK, IWORK, INFO ) 00354 * 00355 * Check error code from DPORFS. 00356 * 00357 IF( INFO.NE.0 ) 00358 $ CALL ALAERH( PATH, 'DPORFS', INFO, 0, UPLO, N, 00359 $ N, -1, -1, NRHS, IMAT, NFAIL, 00360 $ NERRS, NOUT ) 00361 * 00362 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00363 $ RESULT( 5 ) ) 00364 CALL DPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA, 00365 $ XACT, LDA, RWORK, RWORK( NRHS+1 ), 00366 $ RESULT( 6 ) ) 00367 * 00368 * Print information about the tests that did not pass 00369 * the threshold. 00370 * 00371 DO 70 K = 3, 7 00372 IF( RESULT( K ).GE.THRESH ) THEN 00373 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00374 $ CALL ALAHD( NOUT, PATH ) 00375 WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, 00376 $ IMAT, K, RESULT( K ) 00377 NFAIL = NFAIL + 1 00378 END IF 00379 70 CONTINUE 00380 NRUN = NRUN + 5 00381 80 CONTINUE 00382 * 00383 *+ TEST 8 00384 * Get an estimate of RCOND = 1/CNDNUM. 00385 * 00386 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK ) 00387 SRNAMT = 'DPOCON' 00388 CALL DPOCON( UPLO, N, AFAC, LDA, ANORM, RCOND, WORK, 00389 $ IWORK, INFO ) 00390 * 00391 * Check error code from DPOCON. 00392 * 00393 IF( INFO.NE.0 ) 00394 $ CALL ALAERH( PATH, 'DPOCON', INFO, 0, UPLO, N, N, 00395 $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00396 * 00397 RESULT( 8 ) = DGET06( RCOND, RCONDC ) 00398 * 00399 * Print the test ratio if it is .GE. THRESH. 00400 * 00401 IF( RESULT( 8 ).GE.THRESH ) THEN 00402 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00403 $ CALL ALAHD( NOUT, PATH ) 00404 WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 8, 00405 $ RESULT( 8 ) 00406 NFAIL = NFAIL + 1 00407 END IF 00408 NRUN = NRUN + 1 00409 90 CONTINUE 00410 100 CONTINUE 00411 110 CONTINUE 00412 120 CONTINUE 00413 * 00414 * Print a summary of the results. 00415 * 00416 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00417 * 00418 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ', 00419 $ I2, ', test ', I2, ', ratio =', G12.5 ) 00420 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', 00421 $ I2, ', test(', I2, ') =', G12.5 ) 00422 9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2, 00423 $ ', test(', I2, ') =', G12.5 ) 00424 RETURN 00425 * 00426 * End of DCHKPO 00427 * 00428 END
1.7.3 
