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