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