 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DCHKTB( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR, 00002 $ NMAX, AB, AINV, B, X, XACT, WORK, RWORK, IWORK, 00003 $ 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 AB( * ), AINV( * ), B( * ), RWORK( * ), 00018 $ WORK( * ), X( * ), XACT( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * DCHKTB tests DTBTRS, -RFS, and -CON, and DLATBS. 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 column 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 leading dimension of the work arrays. 00056 * NMAX >= the maximum value of N in NVAL. 00057 * 00058 * AB (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX) 00059 * 00060 * AINV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX) 00061 * 00062 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX) 00063 * where NSMAX is the largest entry in NSVAL. 00064 * 00065 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX) 00066 * 00067 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX) 00068 * 00069 * WORK (workspace) DOUBLE PRECISION array, dimension 00070 * (NMAX*max(3,NSMAX)) 00071 * 00072 * RWORK (workspace) DOUBLE PRECISION array, dimension 00073 * (max(NMAX,2*NSMAX)) 00074 * 00075 * IWORK (workspace) INTEGER array, dimension (NMAX) 00076 * 00077 * NOUT (input) INTEGER 00078 * The unit number for output. 00079 * 00080 * ===================================================================== 00081 * 00082 * .. Parameters .. 00083 INTEGER NTYPE1, NTYPES 00084 PARAMETER ( NTYPE1 = 9, NTYPES = 17 ) 00085 INTEGER NTESTS 00086 PARAMETER ( NTESTS = 8 ) 00087 INTEGER NTRAN 00088 PARAMETER ( NTRAN = 3 ) 00089 DOUBLE PRECISION ONE, ZERO 00090 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00091 * .. 00092 * .. Local Scalars .. 00093 CHARACTER DIAG, NORM, TRANS, UPLO, XTYPE 00094 CHARACTER*3 PATH 00095 INTEGER I, IDIAG, IK, IMAT, IN, INFO, IRHS, ITRAN, 00096 $ IUPLO, J, K, KD, LDA, LDAB, N, NERRS, NFAIL, 00097 $ NIMAT, NIMAT2, NK, NRHS, NRUN 00098 DOUBLE PRECISION AINVNM, ANORM, RCOND, RCONDC, RCONDI, RCONDO, 00099 $ SCALE 00100 * .. 00101 * .. Local Arrays .. 00102 CHARACTER TRANSS( NTRAN ), UPLOS( 2 ) 00103 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00104 DOUBLE PRECISION RESULT( NTESTS ) 00105 * .. 00106 * .. External Functions .. 00107 LOGICAL LSAME 00108 DOUBLE PRECISION DLANTB, DLANTR 00109 EXTERNAL LSAME, DLANTB, DLANTR 00110 * .. 00111 * .. External Subroutines .. 00112 EXTERNAL ALAERH, ALAHD, ALASUM, DCOPY, DERRTR, DGET04, 00113 $ DLACPY, DLARHS, DLASET, DLATBS, DLATTB, DTBCON, 00114 $ DTBRFS, DTBSV, DTBT02, DTBT03, DTBT05, DTBT06, 00115 $ DTBTRS 00116 * .. 00117 * .. Scalars in Common .. 00118 LOGICAL LERR, OK 00119 CHARACTER*32 SRNAMT 00120 INTEGER INFOT, IOUNIT 00121 * .. 00122 * .. Common blocks .. 00123 COMMON / INFOC / INFOT, IOUNIT, OK, LERR 00124 COMMON / SRNAMC / SRNAMT 00125 * .. 00126 * .. Intrinsic Functions .. 00127 INTRINSIC MAX, MIN 00128 * .. 00129 * .. Data statements .. 00130 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00131 DATA UPLOS / 'U', 'L' / , TRANSS / 'N', 'T', 'C' / 00132 * .. 00133 * .. Executable Statements .. 00134 * 00135 * Initialize constants and the random number seed. 00136 * 00137 PATH( 1: 1 ) = 'Double precision' 00138 PATH( 2: 3 ) = 'TB' 00139 NRUN = 0 00140 NFAIL = 0 00141 NERRS = 0 00142 DO 10 I = 1, 4 00143 ISEED( I ) = ISEEDY( I ) 00144 10 CONTINUE 00145 * 00146 * Test the error exits 00147 * 00148 IF( TSTERR ) 00149 $ CALL DERRTR( PATH, NOUT ) 00150 INFOT = 0 00151 * 00152 DO 140 IN = 1, NN 00153 * 00154 * Do for each value of N in NVAL 00155 * 00156 N = NVAL( IN ) 00157 LDA = MAX( 1, N ) 00158 XTYPE = 'N' 00159 NIMAT = NTYPE1 00160 NIMAT2 = NTYPES 00161 IF( N.LE.0 ) THEN 00162 NIMAT = 1 00163 NIMAT2 = NTYPE1 + 1 00164 END IF 00165 * 00166 NK = MIN( N+1, 4 ) 00167 DO 130 IK = 1, NK 00168 * 00169 * Do for KD = 0, N, (3N-1)/4, and (N+1)/4. This order makes 00170 * it easier to skip redundant values for small values of N. 00171 * 00172 IF( IK.EQ.1 ) THEN 00173 KD = 0 00174 ELSE IF( IK.EQ.2 ) THEN 00175 KD = MAX( N, 0 ) 00176 ELSE IF( IK.EQ.3 ) THEN 00177 KD = ( 3*N-1 ) / 4 00178 ELSE IF( IK.EQ.4 ) THEN 00179 KD = ( N+1 ) / 4 00180 END IF 00181 LDAB = KD + 1 00182 * 00183 DO 90 IMAT = 1, NIMAT 00184 * 00185 * Do the tests only if DOTYPE( IMAT ) is true. 00186 * 00187 IF( .NOT.DOTYPE( IMAT ) ) 00188 $ GO TO 90 00189 * 00190 DO 80 IUPLO = 1, 2 00191 * 00192 * Do first for UPLO = 'U', then for UPLO = 'L' 00193 * 00194 UPLO = UPLOS( IUPLO ) 00195 * 00196 * Call DLATTB to generate a triangular test matrix. 00197 * 00198 SRNAMT = 'DLATTB' 00199 CALL DLATTB( IMAT, UPLO, 'No transpose', DIAG, ISEED, 00200 $ N, KD, AB, LDAB, X, WORK, INFO ) 00201 * 00202 * Set IDIAG = 1 for non-unit matrices, 2 for unit. 00203 * 00204 IF( LSAME( DIAG, 'N' ) ) THEN 00205 IDIAG = 1 00206 ELSE 00207 IDIAG = 2 00208 END IF 00209 * 00210 * Form the inverse of A so we can get a good estimate 00211 * of RCONDC = 1/(norm(A) * norm(inv(A))). 00212 * 00213 CALL DLASET( 'Full', N, N, ZERO, ONE, AINV, LDA ) 00214 IF( LSAME( UPLO, 'U' ) ) THEN 00215 DO 20 J = 1, N 00216 CALL DTBSV( UPLO, 'No transpose', DIAG, J, KD, 00217 $ AB, LDAB, AINV( ( J-1 )*LDA+1 ), 1 ) 00218 20 CONTINUE 00219 ELSE 00220 DO 30 J = 1, N 00221 CALL DTBSV( UPLO, 'No transpose', DIAG, N-J+1, 00222 $ KD, AB( ( J-1 )*LDAB+1 ), LDAB, 00223 $ AINV( ( J-1 )*LDA+J ), 1 ) 00224 30 CONTINUE 00225 END IF 00226 * 00227 * Compute the 1-norm condition number of A. 00228 * 00229 ANORM = DLANTB( '1', UPLO, DIAG, N, KD, AB, LDAB, 00230 $ RWORK ) 00231 AINVNM = DLANTR( '1', UPLO, DIAG, N, N, AINV, LDA, 00232 $ RWORK ) 00233 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00234 RCONDO = ONE 00235 ELSE 00236 RCONDO = ( ONE / ANORM ) / AINVNM 00237 END IF 00238 * 00239 * Compute the infinity-norm condition number of A. 00240 * 00241 ANORM = DLANTB( 'I', UPLO, DIAG, N, KD, AB, LDAB, 00242 $ RWORK ) 00243 AINVNM = DLANTR( 'I', UPLO, DIAG, N, N, AINV, LDA, 00244 $ RWORK ) 00245 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 00246 RCONDI = ONE 00247 ELSE 00248 RCONDI = ( ONE / ANORM ) / AINVNM 00249 END IF 00250 * 00251 DO 60 IRHS = 1, NNS 00252 NRHS = NSVAL( IRHS ) 00253 XTYPE = 'N' 00254 * 00255 DO 50 ITRAN = 1, NTRAN 00256 * 00257 * Do for op(A) = A, A**T, or A**H. 00258 * 00259 TRANS = TRANSS( ITRAN ) 00260 IF( ITRAN.EQ.1 ) THEN 00261 NORM = 'O' 00262 RCONDC = RCONDO 00263 ELSE 00264 NORM = 'I' 00265 RCONDC = RCONDI 00266 END IF 00267 * 00268 *+ TEST 1 00269 * Solve and compute residual for op(A)*x = b. 00270 * 00271 SRNAMT = 'DLARHS' 00272 CALL DLARHS( PATH, XTYPE, UPLO, TRANS, N, N, KD, 00273 $ IDIAG, NRHS, AB, LDAB, XACT, LDA, 00274 $ B, LDA, ISEED, INFO ) 00275 XTYPE = 'C' 00276 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00277 * 00278 SRNAMT = 'DTBTRS' 00279 CALL DTBTRS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, 00280 $ LDAB, X, LDA, INFO ) 00281 * 00282 * Check error code from DTBTRS. 00283 * 00284 IF( INFO.NE.0 ) 00285 $ CALL ALAERH( PATH, 'DTBTRS', INFO, 0, 00286 $ UPLO // TRANS // DIAG, N, N, KD, 00287 $ KD, NRHS, IMAT, NFAIL, NERRS, 00288 $ NOUT ) 00289 * 00290 CALL DTBT02( UPLO, TRANS, DIAG, N, KD, NRHS, AB, 00291 $ LDAB, X, LDA, B, LDA, WORK, 00292 $ RESULT( 1 ) ) 00293 * 00294 *+ TEST 2 00295 * Check solution from generated exact solution. 00296 * 00297 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00298 $ RESULT( 2 ) ) 00299 * 00300 *+ TESTS 3, 4, and 5 00301 * Use iterative refinement to improve the solution 00302 * and compute error bounds. 00303 * 00304 SRNAMT = 'DTBRFS' 00305 CALL DTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, 00306 $ LDAB, B, LDA, X, LDA, RWORK, 00307 $ RWORK( NRHS+1 ), WORK, IWORK, 00308 $ INFO ) 00309 * 00310 * Check error code from DTBRFS. 00311 * 00312 IF( INFO.NE.0 ) 00313 $ CALL ALAERH( PATH, 'DTBRFS', INFO, 0, 00314 $ UPLO // TRANS // DIAG, N, N, KD, 00315 $ KD, NRHS, IMAT, NFAIL, NERRS, 00316 $ NOUT ) 00317 * 00318 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00319 $ RESULT( 3 ) ) 00320 CALL DTBT05( UPLO, TRANS, DIAG, N, KD, NRHS, AB, 00321 $ LDAB, B, LDA, X, LDA, XACT, LDA, 00322 $ RWORK, RWORK( NRHS+1 ), 00323 $ RESULT( 4 ) ) 00324 * 00325 * Print information about the tests that did not 00326 * pass the threshold. 00327 * 00328 DO 40 K = 1, 5 00329 IF( RESULT( K ).GE.THRESH ) THEN 00330 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00331 $ CALL ALAHD( NOUT, PATH ) 00332 WRITE( NOUT, FMT = 9999 )UPLO, TRANS, 00333 $ DIAG, N, KD, NRHS, IMAT, K, RESULT( K ) 00334 NFAIL = NFAIL + 1 00335 END IF 00336 40 CONTINUE 00337 NRUN = NRUN + 5 00338 50 CONTINUE 00339 60 CONTINUE 00340 * 00341 *+ TEST 6 00342 * Get an estimate of RCOND = 1/CNDNUM. 00343 * 00344 DO 70 ITRAN = 1, 2 00345 IF( ITRAN.EQ.1 ) THEN 00346 NORM = 'O' 00347 RCONDC = RCONDO 00348 ELSE 00349 NORM = 'I' 00350 RCONDC = RCONDI 00351 END IF 00352 SRNAMT = 'DTBCON' 00353 CALL DTBCON( NORM, UPLO, DIAG, N, KD, AB, LDAB, 00354 $ RCOND, WORK, IWORK, INFO ) 00355 * 00356 * Check error code from DTBCON. 00357 * 00358 IF( INFO.NE.0 ) 00359 $ CALL ALAERH( PATH, 'DTBCON', INFO, 0, 00360 $ NORM // UPLO // DIAG, N, N, KD, KD, 00361 $ -1, IMAT, NFAIL, NERRS, NOUT ) 00362 * 00363 CALL DTBT06( RCOND, RCONDC, UPLO, DIAG, N, KD, AB, 00364 $ LDAB, RWORK, RESULT( 6 ) ) 00365 * 00366 * Print information about the tests that did not pass 00367 * the threshold. 00368 * 00369 IF( RESULT( 6 ).GE.THRESH ) THEN 00370 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00371 $ CALL ALAHD( NOUT, PATH ) 00372 WRITE( NOUT, FMT = 9998 ) 'DTBCON', NORM, UPLO, 00373 $ DIAG, N, KD, IMAT, 6, RESULT( 6 ) 00374 NFAIL = NFAIL + 1 00375 END IF 00376 NRUN = NRUN + 1 00377 70 CONTINUE 00378 80 CONTINUE 00379 90 CONTINUE 00380 * 00381 * Use pathological test matrices to test DLATBS. 00382 * 00383 DO 120 IMAT = NTYPE1 + 1, NIMAT2 00384 * 00385 * Do the tests only if DOTYPE( IMAT ) is true. 00386 * 00387 IF( .NOT.DOTYPE( IMAT ) ) 00388 $ GO TO 120 00389 * 00390 DO 110 IUPLO = 1, 2 00391 * 00392 * Do first for UPLO = 'U', then for UPLO = 'L' 00393 * 00394 UPLO = UPLOS( IUPLO ) 00395 DO 100 ITRAN = 1, NTRAN 00396 * 00397 * Do for op(A) = A, A**T, and A**H. 00398 * 00399 TRANS = TRANSS( ITRAN ) 00400 * 00401 * Call DLATTB to generate a triangular test matrix. 00402 * 00403 SRNAMT = 'DLATTB' 00404 CALL DLATTB( IMAT, UPLO, TRANS, DIAG, ISEED, N, KD, 00405 $ AB, LDAB, X, WORK, INFO ) 00406 * 00407 *+ TEST 7 00408 * Solve the system op(A)*x = b 00409 * 00410 SRNAMT = 'DLATBS' 00411 CALL DCOPY( N, X, 1, B, 1 ) 00412 CALL DLATBS( UPLO, TRANS, DIAG, 'N', N, KD, AB, 00413 $ LDAB, B, SCALE, RWORK, INFO ) 00414 * 00415 * Check error code from DLATBS. 00416 * 00417 IF( INFO.NE.0 ) 00418 $ CALL ALAERH( PATH, 'DLATBS', INFO, 0, 00419 $ UPLO // TRANS // DIAG // 'N', N, N, 00420 $ KD, KD, -1, IMAT, NFAIL, NERRS, 00421 $ NOUT ) 00422 * 00423 CALL DTBT03( UPLO, TRANS, DIAG, N, KD, 1, AB, LDAB, 00424 $ SCALE, RWORK, ONE, B, LDA, X, LDA, 00425 $ WORK, RESULT( 7 ) ) 00426 * 00427 *+ TEST 8 00428 * Solve op(A)*x = b again with NORMIN = 'Y'. 00429 * 00430 CALL DCOPY( N, X, 1, B, 1 ) 00431 CALL DLATBS( UPLO, TRANS, DIAG, 'Y', N, KD, AB, 00432 $ LDAB, B, SCALE, RWORK, INFO ) 00433 * 00434 * Check error code from DLATBS. 00435 * 00436 IF( INFO.NE.0 ) 00437 $ CALL ALAERH( PATH, 'DLATBS', INFO, 0, 00438 $ UPLO // TRANS // DIAG // 'Y', N, N, 00439 $ KD, KD, -1, IMAT, NFAIL, NERRS, 00440 $ NOUT ) 00441 * 00442 CALL DTBT03( UPLO, TRANS, DIAG, N, KD, 1, AB, LDAB, 00443 $ SCALE, RWORK, ONE, B, LDA, X, LDA, 00444 $ WORK, RESULT( 8 ) ) 00445 * 00446 * Print information about the tests that did not pass 00447 * the threshold. 00448 * 00449 IF( RESULT( 7 ).GE.THRESH ) THEN 00450 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00451 $ CALL ALAHD( NOUT, PATH ) 00452 WRITE( NOUT, FMT = 9997 )'DLATBS', UPLO, TRANS, 00453 $ DIAG, 'N', N, KD, IMAT, 7, RESULT( 7 ) 00454 NFAIL = NFAIL + 1 00455 END IF 00456 IF( RESULT( 8 ).GE.THRESH ) THEN 00457 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00458 $ CALL ALAHD( NOUT, PATH ) 00459 WRITE( NOUT, FMT = 9997 )'DLATBS', UPLO, TRANS, 00460 $ DIAG, 'Y', N, KD, IMAT, 8, RESULT( 8 ) 00461 NFAIL = NFAIL + 1 00462 END IF 00463 NRUN = NRUN + 2 00464 100 CONTINUE 00465 110 CONTINUE 00466 120 CONTINUE 00467 130 CONTINUE 00468 140 CONTINUE 00469 * 00470 * Print a summary of the results. 00471 * 00472 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00473 * 00474 9999 FORMAT( ' UPLO=''', A1, ''', TRANS=''', A1, '' 00475 ', $ DIAG=''', A1, ''', N=', I5, ', KD=', I5, ', NRHS=', I5, 00476 $ ', type ', I2, ', test(', I2, ')=', G12.5 ) 00477 9998 FORMAT( 1X, A, '( ''', A1, ''', ''', A1, ''', ''', A1, ''',', 00478 $ I5, ',', I5, ', ... ), type ', I2, ', test(', I2, ')=', 00479 $ G12.5 ) 00480 9997 FORMAT( 1X, A, '( ''', A1, ''', ''', A1, ''', ''', A1, ''', ''', 00481 $ A1, ''',', I5, ',', I5, ', ... ), type ', I2, ', test(', 00482 $ I1, ')=', G12.5 ) 00483 RETURN 00484 * 00485 * End of DCHKTB 00486 * 00487 END
1.7.3 
