 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZCHKPB( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, 00002 $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, 00003 $ XACT, WORK, RWORK, 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 NBVAL( * ), NSVAL( * ), NVAL( * ) 00017 DOUBLE PRECISION RWORK( * ) 00018 COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), 00019 $ WORK( * ), X( * ), XACT( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * ZCHKPB tests ZPBTRF, -TRS, -RFS, and -CON. 00026 * 00027 * Arguments 00028 * ========= 00029 * 00030 * DOTYPE (input) LOGICAL array, dimension (NTYPES) 00031 * The matrix types to be used for testing. Matrices of type j 00032 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 00033 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 00034 * 00035 * NN (input) INTEGER 00036 * The number of values of N contained in the vector NVAL. 00037 * 00038 * NVAL (input) INTEGER array, dimension (NN) 00039 * The values of the matrix dimension N. 00040 * 00041 * NNB (input) INTEGER 00042 * The number of values of NB contained in the vector NBVAL. 00043 * 00044 * NBVAL (input) INTEGER array, dimension (NBVAL) 00045 * The values of the blocksize NB. 00046 * 00047 * NNS (input) INTEGER 00048 * The number of values of NRHS contained in the vector NSVAL. 00049 * 00050 * NSVAL (input) INTEGER array, dimension (NNS) 00051 * The values of the number of right hand sides NRHS. 00052 * 00053 * THRESH (input) DOUBLE PRECISION 00054 * The threshold value for the test ratios. A result is 00055 * included in the output file if RESULT >= THRESH. To have 00056 * every test ratio printed, use THRESH = 0. 00057 * 00058 * TSTERR (input) LOGICAL 00059 * Flag that indicates whether error exits are to be tested. 00060 * 00061 * NMAX (input) INTEGER 00062 * The maximum value permitted for N, used in dimensioning the 00063 * work arrays. 00064 * 00065 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX) 00066 * 00067 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX) 00068 * 00069 * AINV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX) 00070 * 00071 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX) 00072 * where NSMAX is the largest entry in NSVAL. 00073 * 00074 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX) 00075 * 00076 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX) 00077 * 00078 * WORK (workspace) DOUBLE PRECISION array, dimension 00079 * (NMAX*max(3,NSMAX)) 00080 * 00081 * RWORK (workspace) DOUBLE PRECISION array, dimension 00082 * (max(NMAX,2*NSMAX)) 00083 * 00084 * NOUT (input) INTEGER 00085 * The unit number for output. 00086 * 00087 * ===================================================================== 00088 * 00089 * .. Parameters .. 00090 DOUBLE PRECISION ONE, ZERO 00091 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00092 INTEGER NTYPES, NTESTS 00093 PARAMETER ( NTYPES = 8, NTESTS = 7 ) 00094 INTEGER NBW 00095 PARAMETER ( NBW = 4 ) 00096 * .. 00097 * .. Local Scalars .. 00098 LOGICAL ZEROT 00099 CHARACTER DIST, PACKIT, TYPE, UPLO, XTYPE 00100 CHARACTER*3 PATH 00101 INTEGER I, I1, I2, IKD, IMAT, IN, INB, INFO, IOFF, 00102 $ IRHS, IUPLO, IW, IZERO, K, KD, KL, KOFF, KU, 00103 $ LDA, LDAB, MODE, N, NB, NERRS, NFAIL, NIMAT, 00104 $ NKD, NRHS, NRUN 00105 DOUBLE PRECISION AINVNM, ANORM, CNDNUM, RCOND, RCONDC 00106 * .. 00107 * .. Local Arrays .. 00108 INTEGER ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW ) 00109 DOUBLE PRECISION RESULT( NTESTS ) 00110 * .. 00111 * .. External Functions .. 00112 DOUBLE PRECISION DGET06, ZLANGE, ZLANHB 00113 EXTERNAL DGET06, ZLANGE, ZLANHB 00114 * .. 00115 * .. External Subroutines .. 00116 EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZCOPY, ZERRPO, 00117 $ ZGET04, ZLACPY, ZLAIPD, ZLARHS, ZLASET, ZLATB4, 00118 $ ZLATMS, ZPBCON, ZPBRFS, ZPBT01, ZPBT02, ZPBT05, 00119 $ ZPBTRF, ZPBTRS, ZSWAP 00120 * .. 00121 * .. Intrinsic Functions .. 00122 INTRINSIC DCMPLX, MAX, MIN 00123 * .. 00124 * .. Scalars in Common .. 00125 LOGICAL LERR, OK 00126 CHARACTER*32 SRNAMT 00127 INTEGER INFOT, NUNIT 00128 * .. 00129 * .. Common blocks .. 00130 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00131 COMMON / SRNAMC / SRNAMT 00132 * .. 00133 * .. Data statements .. 00134 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00135 * .. 00136 * .. Executable Statements .. 00137 * 00138 * Initialize constants and the random number seed. 00139 * 00140 PATH( 1: 1 ) = 'Zomplex precision' 00141 PATH( 2: 3 ) = 'PB' 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 ZERRPO( PATH, NOUT ) 00153 INFOT = 0 00154 KDVAL( 1 ) = 0 00155 * 00156 * Do for each value of N in NVAL 00157 * 00158 DO 90 IN = 1, NN 00159 N = NVAL( IN ) 00160 LDA = MAX( N, 1 ) 00161 XTYPE = 'N' 00162 * 00163 * Set limits on the number of loop iterations. 00164 * 00165 NKD = MAX( 1, MIN( N, 4 ) ) 00166 NIMAT = NTYPES 00167 IF( N.EQ.0 ) 00168 $ NIMAT = 1 00169 * 00170 KDVAL( 2 ) = N + ( N+1 ) / 4 00171 KDVAL( 3 ) = ( 3*N-1 ) / 4 00172 KDVAL( 4 ) = ( N+1 ) / 4 00173 * 00174 DO 80 IKD = 1, NKD 00175 * 00176 * Do for KD = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This order 00177 * makes it easier to skip redundant values for small values 00178 * of N. 00179 * 00180 KD = KDVAL( IKD ) 00181 LDAB = KD + 1 00182 * 00183 * Do first for UPLO = 'U', then for UPLO = 'L' 00184 * 00185 DO 70 IUPLO = 1, 2 00186 KOFF = 1 00187 IF( IUPLO.EQ.1 ) THEN 00188 UPLO = 'U' 00189 KOFF = MAX( 1, KD+2-N ) 00190 PACKIT = 'Q' 00191 ELSE 00192 UPLO = 'L' 00193 PACKIT = 'B' 00194 END IF 00195 * 00196 DO 60 IMAT = 1, NIMAT 00197 * 00198 * Do the tests only if DOTYPE( IMAT ) is true. 00199 * 00200 IF( .NOT.DOTYPE( IMAT ) ) 00201 $ GO TO 60 00202 * 00203 * Skip types 2, 3, or 4 if the matrix size is too small. 00204 * 00205 ZEROT = IMAT.GE.2 .AND. IMAT.LE.4 00206 IF( ZEROT .AND. N.LT.IMAT-1 ) 00207 $ GO TO 60 00208 * 00209 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN 00210 * 00211 * Set up parameters with ZLATB4 and generate a test 00212 * matrix with ZLATMS. 00213 * 00214 CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, 00215 $ MODE, CNDNUM, DIST ) 00216 * 00217 SRNAMT = 'ZLATMS' 00218 CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 00219 $ CNDNUM, ANORM, KD, KD, PACKIT, 00220 $ A( KOFF ), LDAB, WORK, INFO ) 00221 * 00222 * Check error code from ZLATMS. 00223 * 00224 IF( INFO.NE.0 ) THEN 00225 CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, 00226 $ N, KD, KD, -1, IMAT, NFAIL, NERRS, 00227 $ NOUT ) 00228 GO TO 60 00229 END IF 00230 ELSE IF( IZERO.GT.0 ) THEN 00231 * 00232 * Use the same matrix for types 3 and 4 as for type 00233 * 2 by copying back the zeroed out column, 00234 * 00235 IW = 2*LDA + 1 00236 IF( IUPLO.EQ.1 ) THEN 00237 IOFF = ( IZERO-1 )*LDAB + KD + 1 00238 CALL ZCOPY( IZERO-I1, WORK( IW ), 1, 00239 $ A( IOFF-IZERO+I1 ), 1 ) 00240 IW = IW + IZERO - I1 00241 CALL ZCOPY( I2-IZERO+1, WORK( IW ), 1, 00242 $ A( IOFF ), MAX( LDAB-1, 1 ) ) 00243 ELSE 00244 IOFF = ( I1-1 )*LDAB + 1 00245 CALL ZCOPY( IZERO-I1, WORK( IW ), 1, 00246 $ A( IOFF+IZERO-I1 ), 00247 $ MAX( LDAB-1, 1 ) ) 00248 IOFF = ( IZERO-1 )*LDAB + 1 00249 IW = IW + IZERO - I1 00250 CALL ZCOPY( I2-IZERO+1, WORK( IW ), 1, 00251 $ A( IOFF ), 1 ) 00252 END IF 00253 END IF 00254 * 00255 * For types 2-4, zero one row and column of the matrix 00256 * to test that INFO is returned correctly. 00257 * 00258 IZERO = 0 00259 IF( ZEROT ) THEN 00260 IF( IMAT.EQ.2 ) THEN 00261 IZERO = 1 00262 ELSE IF( IMAT.EQ.3 ) THEN 00263 IZERO = N 00264 ELSE 00265 IZERO = N / 2 + 1 00266 END IF 00267 * 00268 * Save the zeroed out row and column in WORK(*,3) 00269 * 00270 IW = 2*LDA 00271 DO 20 I = 1, MIN( 2*KD+1, N ) 00272 WORK( IW+I ) = ZERO 00273 20 CONTINUE 00274 IW = IW + 1 00275 I1 = MAX( IZERO-KD, 1 ) 00276 I2 = MIN( IZERO+KD, N ) 00277 * 00278 IF( IUPLO.EQ.1 ) THEN 00279 IOFF = ( IZERO-1 )*LDAB + KD + 1 00280 CALL ZSWAP( IZERO-I1, A( IOFF-IZERO+I1 ), 1, 00281 $ WORK( IW ), 1 ) 00282 IW = IW + IZERO - I1 00283 CALL ZSWAP( I2-IZERO+1, A( IOFF ), 00284 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 ) 00285 ELSE 00286 IOFF = ( I1-1 )*LDAB + 1 00287 CALL ZSWAP( IZERO-I1, A( IOFF+IZERO-I1 ), 00288 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 ) 00289 IOFF = ( IZERO-1 )*LDAB + 1 00290 IW = IW + IZERO - I1 00291 CALL ZSWAP( I2-IZERO+1, A( IOFF ), 1, 00292 $ WORK( IW ), 1 ) 00293 END IF 00294 END IF 00295 * 00296 * Set the imaginary part of the diagonals. 00297 * 00298 IF( IUPLO.EQ.1 ) THEN 00299 CALL ZLAIPD( N, A( KD+1 ), LDAB, 0 ) 00300 ELSE 00301 CALL ZLAIPD( N, A( 1 ), LDAB, 0 ) 00302 END IF 00303 * 00304 * Do for each value of NB in NBVAL 00305 * 00306 DO 50 INB = 1, NNB 00307 NB = NBVAL( INB ) 00308 CALL XLAENV( 1, NB ) 00309 * 00310 * Compute the L*L' or U'*U factorization of the band 00311 * matrix. 00312 * 00313 CALL ZLACPY( 'Full', KD+1, N, A, LDAB, AFAC, LDAB ) 00314 SRNAMT = 'ZPBTRF' 00315 CALL ZPBTRF( UPLO, N, KD, AFAC, LDAB, INFO ) 00316 * 00317 * Check error code from ZPBTRF. 00318 * 00319 IF( INFO.NE.IZERO ) THEN 00320 CALL ALAERH( PATH, 'ZPBTRF', INFO, IZERO, UPLO, 00321 $ N, N, KD, KD, NB, IMAT, NFAIL, 00322 $ NERRS, NOUT ) 00323 GO TO 50 00324 END IF 00325 * 00326 * Skip the tests if INFO is not 0. 00327 * 00328 IF( INFO.NE.0 ) 00329 $ GO TO 50 00330 * 00331 *+ TEST 1 00332 * Reconstruct matrix from factors and compute 00333 * residual. 00334 * 00335 CALL ZLACPY( 'Full', KD+1, N, AFAC, LDAB, AINV, 00336 $ LDAB ) 00337 CALL ZPBT01( UPLO, N, KD, A, LDAB, AINV, LDAB, 00338 $ RWORK, RESULT( 1 ) ) 00339 * 00340 * Print the test ratio if it is .GE. THRESH. 00341 * 00342 IF( RESULT( 1 ).GE.THRESH ) THEN 00343 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00344 $ CALL ALAHD( NOUT, PATH ) 00345 WRITE( NOUT, FMT = 9999 )UPLO, N, KD, NB, IMAT, 00346 $ 1, RESULT( 1 ) 00347 NFAIL = NFAIL + 1 00348 END IF 00349 NRUN = NRUN + 1 00350 * 00351 * Only do other tests if this is the first blocksize. 00352 * 00353 IF( INB.GT.1 ) 00354 $ GO TO 50 00355 * 00356 * Form the inverse of A so we can get a good estimate 00357 * of RCONDC = 1/(norm(A) * norm(inv(A))). 00358 * 00359 CALL ZLASET( 'Full', N, N, DCMPLX( ZERO ), 00360 $ DCMPLX( ONE ), AINV, LDA ) 00361 SRNAMT = 'ZPBTRS' 00362 CALL ZPBTRS( UPLO, N, KD, N, AFAC, LDAB, AINV, LDA, 00363 $ INFO ) 00364 * 00365 * Compute RCONDC = 1/(norm(A) * norm(inv(A))). 00366 * 00367 ANORM = ZLANHB( '1', UPLO, N, KD, A, LDAB, RWORK ) 00368 AINVNM = ZLANGE( '1', N, N, AINV, LDA, RWORK ) 00369 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00370 RCONDC = ONE 00371 ELSE 00372 RCONDC = ( ONE / ANORM ) / AINVNM 00373 END IF 00374 * 00375 DO 40 IRHS = 1, NNS 00376 NRHS = NSVAL( IRHS ) 00377 * 00378 *+ TEST 2 00379 * Solve and compute residual for A * X = B. 00380 * 00381 SRNAMT = 'ZLARHS' 00382 CALL ZLARHS( PATH, XTYPE, UPLO, ' ', N, N, KD, 00383 $ KD, NRHS, A, LDAB, XACT, LDA, B, 00384 $ LDA, ISEED, INFO ) 00385 CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00386 * 00387 SRNAMT = 'ZPBTRS' 00388 CALL ZPBTRS( UPLO, N, KD, NRHS, AFAC, LDAB, X, 00389 $ LDA, INFO ) 00390 * 00391 * Check error code from ZPBTRS. 00392 * 00393 IF( INFO.NE.0 ) 00394 $ CALL ALAERH( PATH, 'ZPBTRS', INFO, 0, UPLO, 00395 $ N, N, KD, KD, NRHS, IMAT, NFAIL, 00396 $ NERRS, NOUT ) 00397 * 00398 CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, 00399 $ LDA ) 00400 CALL ZPBT02( UPLO, N, KD, NRHS, A, LDAB, X, LDA, 00401 $ WORK, LDA, RWORK, RESULT( 2 ) ) 00402 * 00403 *+ TEST 3 00404 * Check solution from generated exact solution. 00405 * 00406 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00407 $ RESULT( 3 ) ) 00408 * 00409 *+ TESTS 4, 5, and 6 00410 * Use iterative refinement to improve the solution. 00411 * 00412 SRNAMT = 'ZPBRFS' 00413 CALL ZPBRFS( UPLO, N, KD, NRHS, A, LDAB, AFAC, 00414 $ LDAB, B, LDA, X, LDA, RWORK, 00415 $ RWORK( NRHS+1 ), WORK, 00416 $ RWORK( 2*NRHS+1 ), INFO ) 00417 * 00418 * Check error code from ZPBRFS. 00419 * 00420 IF( INFO.NE.0 ) 00421 $ CALL ALAERH( PATH, 'ZPBRFS', INFO, 0, UPLO, 00422 $ N, N, KD, KD, NRHS, IMAT, NFAIL, 00423 $ NERRS, NOUT ) 00424 * 00425 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00426 $ RESULT( 4 ) ) 00427 CALL ZPBT05( UPLO, N, KD, NRHS, A, LDAB, B, LDA, 00428 $ X, LDA, XACT, LDA, RWORK, 00429 $ RWORK( NRHS+1 ), RESULT( 5 ) ) 00430 * 00431 * Print information about the tests that did not 00432 * pass the threshold. 00433 * 00434 DO 30 K = 2, 6 00435 IF( RESULT( K ).GE.THRESH ) THEN 00436 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00437 $ CALL ALAHD( NOUT, PATH ) 00438 WRITE( NOUT, FMT = 9998 )UPLO, N, KD, 00439 $ NRHS, IMAT, K, RESULT( K ) 00440 NFAIL = NFAIL + 1 00441 END IF 00442 30 CONTINUE 00443 NRUN = NRUN + 5 00444 40 CONTINUE 00445 * 00446 *+ TEST 7 00447 * Get an estimate of RCOND = 1/CNDNUM. 00448 * 00449 SRNAMT = 'ZPBCON' 00450 CALL ZPBCON( UPLO, N, KD, AFAC, LDAB, ANORM, RCOND, 00451 $ WORK, RWORK, INFO ) 00452 * 00453 * Check error code from ZPBCON. 00454 * 00455 IF( INFO.NE.0 ) 00456 $ CALL ALAERH( PATH, 'ZPBCON', INFO, 0, UPLO, N, 00457 $ N, KD, KD, -1, IMAT, NFAIL, NERRS, 00458 $ NOUT ) 00459 * 00460 RESULT( 7 ) = DGET06( RCOND, RCONDC ) 00461 * 00462 * Print the test ratio if it is .GE. THRESH. 00463 * 00464 IF( RESULT( 7 ).GE.THRESH ) THEN 00465 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00466 $ CALL ALAHD( NOUT, PATH ) 00467 WRITE( NOUT, FMT = 9997 )UPLO, N, KD, IMAT, 7, 00468 $ RESULT( 7 ) 00469 NFAIL = NFAIL + 1 00470 END IF 00471 NRUN = NRUN + 1 00472 50 CONTINUE 00473 60 CONTINUE 00474 70 CONTINUE 00475 80 CONTINUE 00476 90 CONTINUE 00477 * 00478 * Print a summary of the results. 00479 * 00480 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00481 * 00482 9999 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NB=', I4, 00483 $ ', type ', I2, ', test ', I2, ', ratio= ', G12.5 ) 00484 9998 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NRHS=', I3, 00485 $ ', type ', I2, ', test(', I2, ') = ', G12.5 ) 00486 9997 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ',', 10X, 00487 $ ' type ', I2, ', test(', I2, ') = ', G12.5 ) 00488 RETURN 00489 * 00490 * End of ZCHKPB 00491 * 00492 END
1.7.3 
