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