 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DCHKPB( 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 * DCHKPB tests DPBTRF, -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 ONE, ZERO 00092 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00093 INTEGER NTYPES, NTESTS 00094 PARAMETER ( NTYPES = 8, NTESTS = 7 ) 00095 INTEGER NBW 00096 PARAMETER ( NBW = 4 ) 00097 * .. 00098 * .. Local Scalars .. 00099 LOGICAL ZEROT 00100 CHARACTER DIST, PACKIT, TYPE, UPLO, XTYPE 00101 CHARACTER*3 PATH 00102 INTEGER I, I1, I2, IKD, IMAT, IN, INB, INFO, IOFF, 00103 $ IRHS, IUPLO, IW, IZERO, K, KD, KL, KOFF, KU, 00104 $ LDA, LDAB, MODE, N, NB, NERRS, NFAIL, NIMAT, 00105 $ NKD, NRHS, NRUN 00106 DOUBLE PRECISION AINVNM, ANORM, CNDNUM, RCOND, RCONDC 00107 * .. 00108 * .. Local Arrays .. 00109 INTEGER ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW ) 00110 DOUBLE PRECISION RESULT( NTESTS ) 00111 * .. 00112 * .. External Functions .. 00113 DOUBLE PRECISION DGET06, DLANGE, DLANSB 00114 EXTERNAL DGET06, DLANGE, DLANSB 00115 * .. 00116 * .. External Subroutines .. 00117 EXTERNAL ALAERH, ALAHD, ALASUM, DCOPY, DERRPO, DGET04, 00118 $ DLACPY, DLARHS, DLASET, DLATB4, DLATMS, DPBCON, 00119 $ DPBRFS, DPBT01, DPBT02, DPBT05, DPBTRF, DPBTRS, 00120 $ DSWAP, XLAENV 00121 * .. 00122 * .. Intrinsic Functions .. 00123 INTRINSIC MAX, MIN 00124 * .. 00125 * .. Scalars in Common .. 00126 LOGICAL LERR, OK 00127 CHARACTER*32 SRNAMT 00128 INTEGER INFOT, NUNIT 00129 * .. 00130 * .. Common blocks .. 00131 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00132 COMMON / SRNAMC / SRNAMT 00133 * .. 00134 * .. Data statements .. 00135 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00136 * .. 00137 * .. Executable Statements .. 00138 * 00139 * Initialize constants and the random number seed. 00140 * 00141 PATH( 1: 1 ) = 'Double precision' 00142 PATH( 2: 3 ) = 'PB' 00143 NRUN = 0 00144 NFAIL = 0 00145 NERRS = 0 00146 DO 10 I = 1, 4 00147 ISEED( I ) = ISEEDY( I ) 00148 10 CONTINUE 00149 * 00150 * Test the error exits 00151 * 00152 IF( TSTERR ) 00153 $ CALL DERRPO( PATH, NOUT ) 00154 INFOT = 0 00155 CALL XLAENV( 2, 2 ) 00156 KDVAL( 1 ) = 0 00157 * 00158 * Do for each value of N in NVAL 00159 * 00160 DO 90 IN = 1, NN 00161 N = NVAL( IN ) 00162 LDA = MAX( N, 1 ) 00163 XTYPE = 'N' 00164 * 00165 * Set limits on the number of loop iterations. 00166 * 00167 NKD = MAX( 1, MIN( N, 4 ) ) 00168 NIMAT = NTYPES 00169 IF( N.EQ.0 ) 00170 $ NIMAT = 1 00171 * 00172 KDVAL( 2 ) = N + ( N+1 ) / 4 00173 KDVAL( 3 ) = ( 3*N-1 ) / 4 00174 KDVAL( 4 ) = ( N+1 ) / 4 00175 * 00176 DO 80 IKD = 1, NKD 00177 * 00178 * Do for KD = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This order 00179 * makes it easier to skip redundant values for small values 00180 * of N. 00181 * 00182 KD = KDVAL( IKD ) 00183 LDAB = KD + 1 00184 * 00185 * Do first for UPLO = 'U', then for UPLO = 'L' 00186 * 00187 DO 70 IUPLO = 1, 2 00188 KOFF = 1 00189 IF( IUPLO.EQ.1 ) THEN 00190 UPLO = 'U' 00191 KOFF = MAX( 1, KD+2-N ) 00192 PACKIT = 'Q' 00193 ELSE 00194 UPLO = 'L' 00195 PACKIT = 'B' 00196 END IF 00197 * 00198 DO 60 IMAT = 1, NIMAT 00199 * 00200 * Do the tests only if DOTYPE( IMAT ) is true. 00201 * 00202 IF( .NOT.DOTYPE( IMAT ) ) 00203 $ GO TO 60 00204 * 00205 * Skip types 2, 3, or 4 if the matrix size is too small. 00206 * 00207 ZEROT = IMAT.GE.2 .AND. IMAT.LE.4 00208 IF( ZEROT .AND. N.LT.IMAT-1 ) 00209 $ GO TO 60 00210 * 00211 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN 00212 * 00213 * Set up parameters with DLATB4 and generate a test 00214 * matrix with DLATMS. 00215 * 00216 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, 00217 $ MODE, CNDNUM, DIST ) 00218 * 00219 SRNAMT = 'DLATMS' 00220 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 00221 $ CNDNUM, ANORM, KD, KD, PACKIT, 00222 $ A( KOFF ), LDAB, WORK, INFO ) 00223 * 00224 * Check error code from DLATMS. 00225 * 00226 IF( INFO.NE.0 ) THEN 00227 CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, 00228 $ N, KD, KD, -1, IMAT, NFAIL, NERRS, 00229 $ NOUT ) 00230 GO TO 60 00231 END IF 00232 ELSE IF( IZERO.GT.0 ) THEN 00233 * 00234 * Use the same matrix for types 3 and 4 as for type 00235 * 2 by copying back the zeroed out column, 00236 * 00237 IW = 2*LDA + 1 00238 IF( IUPLO.EQ.1 ) THEN 00239 IOFF = ( IZERO-1 )*LDAB + KD + 1 00240 CALL DCOPY( IZERO-I1, WORK( IW ), 1, 00241 $ A( IOFF-IZERO+I1 ), 1 ) 00242 IW = IW + IZERO - I1 00243 CALL DCOPY( I2-IZERO+1, WORK( IW ), 1, 00244 $ A( IOFF ), MAX( LDAB-1, 1 ) ) 00245 ELSE 00246 IOFF = ( I1-1 )*LDAB + 1 00247 CALL DCOPY( IZERO-I1, WORK( IW ), 1, 00248 $ A( IOFF+IZERO-I1 ), 00249 $ MAX( LDAB-1, 1 ) ) 00250 IOFF = ( IZERO-1 )*LDAB + 1 00251 IW = IW + IZERO - I1 00252 CALL DCOPY( I2-IZERO+1, WORK( IW ), 1, 00253 $ A( IOFF ), 1 ) 00254 END IF 00255 END IF 00256 * 00257 * For types 2-4, zero one row and column of the matrix 00258 * to test that INFO is returned correctly. 00259 * 00260 IZERO = 0 00261 IF( ZEROT ) THEN 00262 IF( IMAT.EQ.2 ) THEN 00263 IZERO = 1 00264 ELSE IF( IMAT.EQ.3 ) THEN 00265 IZERO = N 00266 ELSE 00267 IZERO = N / 2 + 1 00268 END IF 00269 * 00270 * Save the zeroed out row and column in WORK(*,3) 00271 * 00272 IW = 2*LDA 00273 DO 20 I = 1, MIN( 2*KD+1, N ) 00274 WORK( IW+I ) = ZERO 00275 20 CONTINUE 00276 IW = IW + 1 00277 I1 = MAX( IZERO-KD, 1 ) 00278 I2 = MIN( IZERO+KD, N ) 00279 * 00280 IF( IUPLO.EQ.1 ) THEN 00281 IOFF = ( IZERO-1 )*LDAB + KD + 1 00282 CALL DSWAP( IZERO-I1, A( IOFF-IZERO+I1 ), 1, 00283 $ WORK( IW ), 1 ) 00284 IW = IW + IZERO - I1 00285 CALL DSWAP( I2-IZERO+1, A( IOFF ), 00286 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 ) 00287 ELSE 00288 IOFF = ( I1-1 )*LDAB + 1 00289 CALL DSWAP( IZERO-I1, A( IOFF+IZERO-I1 ), 00290 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 ) 00291 IOFF = ( IZERO-1 )*LDAB + 1 00292 IW = IW + IZERO - I1 00293 CALL DSWAP( I2-IZERO+1, A( IOFF ), 1, 00294 $ WORK( IW ), 1 ) 00295 END IF 00296 END IF 00297 * 00298 * Do for each value of NB in NBVAL 00299 * 00300 DO 50 INB = 1, NNB 00301 NB = NBVAL( INB ) 00302 CALL XLAENV( 1, NB ) 00303 * 00304 * Compute the L*L' or U'*U factorization of the band 00305 * matrix. 00306 * 00307 CALL DLACPY( 'Full', KD+1, N, A, LDAB, AFAC, LDAB ) 00308 SRNAMT = 'DPBTRF' 00309 CALL DPBTRF( UPLO, N, KD, AFAC, LDAB, INFO ) 00310 * 00311 * Check error code from DPBTRF. 00312 * 00313 IF( INFO.NE.IZERO ) THEN 00314 CALL ALAERH( PATH, 'DPBTRF', INFO, IZERO, UPLO, 00315 $ N, N, KD, KD, NB, IMAT, NFAIL, 00316 $ NERRS, NOUT ) 00317 GO TO 50 00318 END IF 00319 * 00320 * Skip the tests if INFO is not 0. 00321 * 00322 IF( INFO.NE.0 ) 00323 $ GO TO 50 00324 * 00325 *+ TEST 1 00326 * Reconstruct matrix from factors and compute 00327 * residual. 00328 * 00329 CALL DLACPY( 'Full', KD+1, N, AFAC, LDAB, AINV, 00330 $ LDAB ) 00331 CALL DPBT01( UPLO, N, KD, A, LDAB, AINV, LDAB, 00332 $ RWORK, RESULT( 1 ) ) 00333 * 00334 * Print the test ratio if it is .GE. THRESH. 00335 * 00336 IF( RESULT( 1 ).GE.THRESH ) THEN 00337 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00338 $ CALL ALAHD( NOUT, PATH ) 00339 WRITE( NOUT, FMT = 9999 )UPLO, N, KD, NB, IMAT, 00340 $ 1, RESULT( 1 ) 00341 NFAIL = NFAIL + 1 00342 END IF 00343 NRUN = NRUN + 1 00344 * 00345 * Only do other tests if this is the first blocksize. 00346 * 00347 IF( INB.GT.1 ) 00348 $ GO TO 50 00349 * 00350 * Form the inverse of A so we can get a good estimate 00351 * of RCONDC = 1/(norm(A) * norm(inv(A))). 00352 * 00353 CALL DLASET( 'Full', N, N, ZERO, ONE, AINV, LDA ) 00354 SRNAMT = 'DPBTRS' 00355 CALL DPBTRS( UPLO, N, KD, N, AFAC, LDAB, AINV, LDA, 00356 $ INFO ) 00357 * 00358 * Compute RCONDC = 1/(norm(A) * norm(inv(A))). 00359 * 00360 ANORM = DLANSB( '1', UPLO, N, KD, A, LDAB, RWORK ) 00361 AINVNM = DLANGE( '1', N, N, AINV, LDA, RWORK ) 00362 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00363 RCONDC = ONE 00364 ELSE 00365 RCONDC = ( ONE / ANORM ) / AINVNM 00366 END IF 00367 * 00368 DO 40 IRHS = 1, NNS 00369 NRHS = NSVAL( IRHS ) 00370 * 00371 *+ TEST 2 00372 * Solve and compute residual for A * X = B. 00373 * 00374 SRNAMT = 'DLARHS' 00375 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KD, 00376 $ KD, NRHS, A, LDAB, XACT, LDA, B, 00377 $ LDA, ISEED, INFO ) 00378 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00379 * 00380 SRNAMT = 'DPBTRS' 00381 CALL DPBTRS( UPLO, N, KD, NRHS, AFAC, LDAB, X, 00382 $ LDA, INFO ) 00383 * 00384 * Check error code from DPBTRS. 00385 * 00386 IF( INFO.NE.0 ) 00387 $ CALL ALAERH( PATH, 'DPBTRS', INFO, 0, UPLO, 00388 $ N, N, KD, KD, NRHS, IMAT, NFAIL, 00389 $ NERRS, NOUT ) 00390 * 00391 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, 00392 $ LDA ) 00393 CALL DPBT02( UPLO, N, KD, NRHS, A, LDAB, X, LDA, 00394 $ WORK, LDA, RWORK, RESULT( 2 ) ) 00395 * 00396 *+ TEST 3 00397 * Check solution from generated exact solution. 00398 * 00399 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00400 $ RESULT( 3 ) ) 00401 * 00402 *+ TESTS 4, 5, and 6 00403 * Use iterative refinement to improve the solution. 00404 * 00405 SRNAMT = 'DPBRFS' 00406 CALL DPBRFS( UPLO, N, KD, NRHS, A, LDAB, AFAC, 00407 $ LDAB, B, LDA, X, LDA, RWORK, 00408 $ RWORK( NRHS+1 ), WORK, IWORK, 00409 $ INFO ) 00410 * 00411 * Check error code from DPBRFS. 00412 * 00413 IF( INFO.NE.0 ) 00414 $ CALL ALAERH( PATH, 'DPBRFS', INFO, 0, UPLO, 00415 $ N, N, KD, KD, NRHS, IMAT, NFAIL, 00416 $ NERRS, NOUT ) 00417 * 00418 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00419 $ RESULT( 4 ) ) 00420 CALL DPBT05( UPLO, N, KD, NRHS, A, LDAB, B, LDA, 00421 $ X, LDA, XACT, LDA, RWORK, 00422 $ RWORK( NRHS+1 ), RESULT( 5 ) ) 00423 * 00424 * Print information about the tests that did not 00425 * pass the threshold. 00426 * 00427 DO 30 K = 2, 6 00428 IF( RESULT( K ).GE.THRESH ) THEN 00429 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00430 $ CALL ALAHD( NOUT, PATH ) 00431 WRITE( NOUT, FMT = 9998 )UPLO, N, KD, 00432 $ NRHS, IMAT, K, RESULT( K ) 00433 NFAIL = NFAIL + 1 00434 END IF 00435 30 CONTINUE 00436 NRUN = NRUN + 5 00437 40 CONTINUE 00438 * 00439 *+ TEST 7 00440 * Get an estimate of RCOND = 1/CNDNUM. 00441 * 00442 SRNAMT = 'DPBCON' 00443 CALL DPBCON( UPLO, N, KD, AFAC, LDAB, ANORM, RCOND, 00444 $ WORK, IWORK, INFO ) 00445 * 00446 * Check error code from DPBCON. 00447 * 00448 IF( INFO.NE.0 ) 00449 $ CALL ALAERH( PATH, 'DPBCON', INFO, 0, UPLO, N, 00450 $ N, KD, KD, -1, IMAT, NFAIL, NERRS, 00451 $ NOUT ) 00452 * 00453 RESULT( 7 ) = DGET06( RCOND, RCONDC ) 00454 * 00455 * Print the test ratio if it is .GE. THRESH. 00456 * 00457 IF( RESULT( 7 ).GE.THRESH ) THEN 00458 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00459 $ CALL ALAHD( NOUT, PATH ) 00460 WRITE( NOUT, FMT = 9997 )UPLO, N, KD, IMAT, 7, 00461 $ RESULT( 7 ) 00462 NFAIL = NFAIL + 1 00463 END IF 00464 NRUN = NRUN + 1 00465 50 CONTINUE 00466 60 CONTINUE 00467 70 CONTINUE 00468 80 CONTINUE 00469 90 CONTINUE 00470 * 00471 * Print a summary of the results. 00472 * 00473 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00474 * 00475 9999 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NB=', I4, 00476 $ ', type ', I2, ', test ', I2, ', ratio= ', G12.5 ) 00477 9998 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NRHS=', I3, 00478 $ ', type ', I2, ', test(', I2, ') = ', G12.5 ) 00479 9997 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ',', 10X, 00480 $ ' type ', I2, ', test(', I2, ') = ', G12.5 ) 00481 RETURN 00482 * 00483 * End of DCHKPB 00484 * 00485 END
1.7.3 
