 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZCHKPT( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR, 00002 $ A, D, E, B, X, XACT, WORK, RWORK, NOUT ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 LOGICAL TSTERR 00010 INTEGER NN, NNS, NOUT 00011 DOUBLE PRECISION THRESH 00012 * .. 00013 * .. Array Arguments .. 00014 LOGICAL DOTYPE( * ) 00015 INTEGER NSVAL( * ), NVAL( * ) 00016 DOUBLE PRECISION D( * ), RWORK( * ) 00017 COMPLEX*16 A( * ), B( * ), E( * ), WORK( * ), X( * ), 00018 $ XACT( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * ZCHKPT tests ZPTTRF, -TRS, -RFS, and -CON 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 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 * A (workspace) COMPLEX*16 array, dimension (NMAX*2) 00055 * 00056 * D (workspace) DOUBLE PRECISION array, dimension (NMAX*2) 00057 * 00058 * E (workspace) COMPLEX*16 array, dimension (NMAX*2) 00059 * 00060 * B (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX) 00061 * where NSMAX is the largest entry in NSVAL. 00062 * 00063 * X (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX) 00064 * 00065 * XACT (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX) 00066 * 00067 * WORK (workspace) COMPLEX*16 array, dimension 00068 * (NMAX*max(3,NSMAX)) 00069 * 00070 * RWORK (workspace) DOUBLE PRECISION array, dimension 00071 * (max(NMAX,2*NSMAX)) 00072 * 00073 * NOUT (input) INTEGER 00074 * The unit number for output. 00075 * 00076 * ===================================================================== 00077 * 00078 * .. Parameters .. 00079 DOUBLE PRECISION ONE, ZERO 00080 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00081 INTEGER NTYPES 00082 PARAMETER ( NTYPES = 12 ) 00083 INTEGER NTESTS 00084 PARAMETER ( NTESTS = 7 ) 00085 * .. 00086 * .. Local Scalars .. 00087 LOGICAL ZEROT 00088 CHARACTER DIST, TYPE, UPLO 00089 CHARACTER*3 PATH 00090 INTEGER I, IA, IMAT, IN, INFO, IRHS, IUPLO, IX, IZERO, 00091 $ J, K, KL, KU, LDA, MODE, N, NERRS, NFAIL, 00092 $ NIMAT, NRHS, NRUN 00093 DOUBLE PRECISION AINVNM, ANORM, COND, DMAX, RCOND, RCONDC 00094 * .. 00095 * .. Local Arrays .. 00096 CHARACTER UPLOS( 2 ) 00097 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00098 DOUBLE PRECISION RESULT( NTESTS ) 00099 COMPLEX*16 Z( 3 ) 00100 * .. 00101 * .. External Functions .. 00102 INTEGER IDAMAX 00103 DOUBLE PRECISION DGET06, DZASUM, ZLANHT 00104 EXTERNAL IDAMAX, DGET06, DZASUM, ZLANHT 00105 * .. 00106 * .. External Subroutines .. 00107 EXTERNAL ALAERH, ALAHD, ALASUM, DCOPY, DLARNV, DSCAL, 00108 $ ZCOPY, ZDSCAL, ZERRGT, ZGET04, ZLACPY, ZLAPTM, 00109 $ ZLARNV, ZLATB4, ZLATMS, ZPTCON, ZPTRFS, ZPTT01, 00110 $ ZPTT02, ZPTT05, ZPTTRF, ZPTTRS 00111 * .. 00112 * .. Intrinsic Functions .. 00113 INTRINSIC ABS, DBLE, MAX 00114 * .. 00115 * .. Scalars in Common .. 00116 LOGICAL LERR, OK 00117 CHARACTER*32 SRNAMT 00118 INTEGER INFOT, NUNIT 00119 * .. 00120 * .. Common blocks .. 00121 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00122 COMMON / SRNAMC / SRNAMT 00123 * .. 00124 * .. Data statements .. 00125 DATA ISEEDY / 0, 0, 0, 1 / , UPLOS / 'U', 'L' / 00126 * .. 00127 * .. Executable Statements .. 00128 * 00129 PATH( 1: 1 ) = 'Zomplex precision' 00130 PATH( 2: 3 ) = 'PT' 00131 NRUN = 0 00132 NFAIL = 0 00133 NERRS = 0 00134 DO 10 I = 1, 4 00135 ISEED( I ) = ISEEDY( I ) 00136 10 CONTINUE 00137 * 00138 * Test the error exits 00139 * 00140 IF( TSTERR ) 00141 $ CALL ZERRGT( PATH, NOUT ) 00142 INFOT = 0 00143 * 00144 DO 120 IN = 1, NN 00145 * 00146 * Do for each value of N in NVAL. 00147 * 00148 N = NVAL( IN ) 00149 LDA = MAX( 1, N ) 00150 NIMAT = NTYPES 00151 IF( N.LE.0 ) 00152 $ NIMAT = 1 00153 * 00154 DO 110 IMAT = 1, NIMAT 00155 * 00156 * Do the tests only if DOTYPE( IMAT ) is true. 00157 * 00158 IF( N.GT.0 .AND. .NOT.DOTYPE( IMAT ) ) 00159 $ GO TO 110 00160 * 00161 * Set up parameters with ZLATB4. 00162 * 00163 CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00164 $ COND, DIST ) 00165 * 00166 ZEROT = IMAT.GE.8 .AND. IMAT.LE.10 00167 IF( IMAT.LE.6 ) THEN 00168 * 00169 * Type 1-6: generate a Hermitian tridiagonal matrix of 00170 * known condition number in lower triangular band storage. 00171 * 00172 SRNAMT = 'ZLATMS' 00173 CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, COND, 00174 $ ANORM, KL, KU, 'B', A, 2, WORK, INFO ) 00175 * 00176 * Check the error code from ZLATMS. 00177 * 00178 IF( INFO.NE.0 ) THEN 00179 CALL ALAERH( PATH, 'ZLATMS', INFO, 0, ' ', N, N, KL, 00180 $ KU, -1, IMAT, NFAIL, NERRS, NOUT ) 00181 GO TO 110 00182 END IF 00183 IZERO = 0 00184 * 00185 * Copy the matrix to D and E. 00186 * 00187 IA = 1 00188 DO 20 I = 1, N - 1 00189 D( I ) = DBLE( A( IA ) ) 00190 E( I ) = A( IA+1 ) 00191 IA = IA + 2 00192 20 CONTINUE 00193 IF( N.GT.0 ) 00194 $ D( N ) = DBLE( A( IA ) ) 00195 ELSE 00196 * 00197 * Type 7-12: generate a diagonally dominant matrix with 00198 * unknown condition number in the vectors D and E. 00199 * 00200 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 7 ) ) THEN 00201 * 00202 * Let E be complex, D real, with values from [-1,1]. 00203 * 00204 CALL DLARNV( 2, ISEED, N, D ) 00205 CALL ZLARNV( 2, ISEED, N-1, E ) 00206 * 00207 * Make the tridiagonal matrix diagonally dominant. 00208 * 00209 IF( N.EQ.1 ) THEN 00210 D( 1 ) = ABS( D( 1 ) ) 00211 ELSE 00212 D( 1 ) = ABS( D( 1 ) ) + ABS( E( 1 ) ) 00213 D( N ) = ABS( D( N ) ) + ABS( E( N-1 ) ) 00214 DO 30 I = 2, N - 1 00215 D( I ) = ABS( D( I ) ) + ABS( E( I ) ) + 00216 $ ABS( E( I-1 ) ) 00217 30 CONTINUE 00218 END IF 00219 * 00220 * Scale D and E so the maximum element is ANORM. 00221 * 00222 IX = IDAMAX( N, D, 1 ) 00223 DMAX = D( IX ) 00224 CALL DSCAL( N, ANORM / DMAX, D, 1 ) 00225 CALL ZDSCAL( N-1, ANORM / DMAX, E, 1 ) 00226 * 00227 ELSE IF( IZERO.GT.0 ) THEN 00228 * 00229 * Reuse the last matrix by copying back the zeroed out 00230 * elements. 00231 * 00232 IF( IZERO.EQ.1 ) THEN 00233 D( 1 ) = Z( 2 ) 00234 IF( N.GT.1 ) 00235 $ E( 1 ) = Z( 3 ) 00236 ELSE IF( IZERO.EQ.N ) THEN 00237 E( N-1 ) = Z( 1 ) 00238 D( N ) = Z( 2 ) 00239 ELSE 00240 E( IZERO-1 ) = Z( 1 ) 00241 D( IZERO ) = Z( 2 ) 00242 E( IZERO ) = Z( 3 ) 00243 END IF 00244 END IF 00245 * 00246 * For types 8-10, set one row and column of the matrix to 00247 * zero. 00248 * 00249 IZERO = 0 00250 IF( IMAT.EQ.8 ) THEN 00251 IZERO = 1 00252 Z( 2 ) = D( 1 ) 00253 D( 1 ) = ZERO 00254 IF( N.GT.1 ) THEN 00255 Z( 3 ) = E( 1 ) 00256 E( 1 ) = ZERO 00257 END IF 00258 ELSE IF( IMAT.EQ.9 ) THEN 00259 IZERO = N 00260 IF( N.GT.1 ) THEN 00261 Z( 1 ) = E( N-1 ) 00262 E( N-1 ) = ZERO 00263 END IF 00264 Z( 2 ) = D( N ) 00265 D( N ) = ZERO 00266 ELSE IF( IMAT.EQ.10 ) THEN 00267 IZERO = ( N+1 ) / 2 00268 IF( IZERO.GT.1 ) THEN 00269 Z( 1 ) = E( IZERO-1 ) 00270 Z( 3 ) = E( IZERO ) 00271 E( IZERO-1 ) = ZERO 00272 E( IZERO ) = ZERO 00273 END IF 00274 Z( 2 ) = D( IZERO ) 00275 D( IZERO ) = ZERO 00276 END IF 00277 END IF 00278 * 00279 CALL DCOPY( N, D, 1, D( N+1 ), 1 ) 00280 IF( N.GT.1 ) 00281 $ CALL ZCOPY( N-1, E, 1, E( N+1 ), 1 ) 00282 * 00283 *+ TEST 1 00284 * Factor A as L*D*L' and compute the ratio 00285 * norm(L*D*L' - A) / (n * norm(A) * EPS ) 00286 * 00287 CALL ZPTTRF( N, D( N+1 ), E( N+1 ), INFO ) 00288 * 00289 * Check error code from ZPTTRF. 00290 * 00291 IF( INFO.NE.IZERO ) THEN 00292 CALL ALAERH( PATH, 'ZPTTRF', INFO, IZERO, ' ', N, N, -1, 00293 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00294 GO TO 110 00295 END IF 00296 * 00297 IF( INFO.GT.0 ) THEN 00298 RCONDC = ZERO 00299 GO TO 100 00300 END IF 00301 * 00302 CALL ZPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK, 00303 $ RESULT( 1 ) ) 00304 * 00305 * Print the test ratio if greater than or equal to THRESH. 00306 * 00307 IF( RESULT( 1 ).GE.THRESH ) THEN 00308 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00309 $ CALL ALAHD( NOUT, PATH ) 00310 WRITE( NOUT, FMT = 9999 )N, IMAT, 1, RESULT( 1 ) 00311 NFAIL = NFAIL + 1 00312 END IF 00313 NRUN = NRUN + 1 00314 * 00315 * Compute RCONDC = 1 / (norm(A) * norm(inv(A)) 00316 * 00317 * Compute norm(A). 00318 * 00319 ANORM = ZLANHT( '1', N, D, E ) 00320 * 00321 * Use ZPTTRS to solve for one column at a time of inv(A), 00322 * computing the maximum column sum as we go. 00323 * 00324 AINVNM = ZERO 00325 DO 50 I = 1, N 00326 DO 40 J = 1, N 00327 X( J ) = ZERO 00328 40 CONTINUE 00329 X( I ) = ONE 00330 CALL ZPTTRS( 'Lower', N, 1, D( N+1 ), E( N+1 ), X, LDA, 00331 $ INFO ) 00332 AINVNM = MAX( AINVNM, DZASUM( N, X, 1 ) ) 00333 50 CONTINUE 00334 RCONDC = ONE / MAX( ONE, ANORM*AINVNM ) 00335 * 00336 DO 90 IRHS = 1, NNS 00337 NRHS = NSVAL( IRHS ) 00338 * 00339 * Generate NRHS random solution vectors. 00340 * 00341 IX = 1 00342 DO 60 J = 1, NRHS 00343 CALL ZLARNV( 2, ISEED, N, XACT( IX ) ) 00344 IX = IX + LDA 00345 60 CONTINUE 00346 * 00347 DO 80 IUPLO = 1, 2 00348 * 00349 * Do first for UPLO = 'U', then for UPLO = 'L'. 00350 * 00351 UPLO = UPLOS( IUPLO ) 00352 * 00353 * Set the right hand side. 00354 * 00355 CALL ZLAPTM( UPLO, N, NRHS, ONE, D, E, XACT, LDA, 00356 $ ZERO, B, LDA ) 00357 * 00358 *+ TEST 2 00359 * Solve A*x = b and compute the residual. 00360 * 00361 CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00362 CALL ZPTTRS( UPLO, N, NRHS, D( N+1 ), E( N+1 ), X, 00363 $ LDA, INFO ) 00364 * 00365 * Check error code from ZPTTRS. 00366 * 00367 IF( INFO.NE.0 ) 00368 $ CALL ALAERH( PATH, 'ZPTTRS', INFO, 0, UPLO, N, N, 00369 $ -1, -1, NRHS, IMAT, NFAIL, NERRS, 00370 $ NOUT ) 00371 * 00372 CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) 00373 CALL ZPTT02( UPLO, N, NRHS, D, E, X, LDA, WORK, LDA, 00374 $ RESULT( 2 ) ) 00375 * 00376 *+ TEST 3 00377 * Check solution from generated exact solution. 00378 * 00379 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00380 $ RESULT( 3 ) ) 00381 * 00382 *+ TESTS 4, 5, and 6 00383 * Use iterative refinement to improve the solution. 00384 * 00385 SRNAMT = 'ZPTRFS' 00386 CALL ZPTRFS( UPLO, N, NRHS, D, E, D( N+1 ), E( N+1 ), 00387 $ B, LDA, X, LDA, RWORK, RWORK( NRHS+1 ), 00388 $ WORK, RWORK( 2*NRHS+1 ), INFO ) 00389 * 00390 * Check error code from ZPTRFS. 00391 * 00392 IF( INFO.NE.0 ) 00393 $ CALL ALAERH( PATH, 'ZPTRFS', INFO, 0, UPLO, N, N, 00394 $ -1, -1, NRHS, IMAT, NFAIL, NERRS, 00395 $ NOUT ) 00396 * 00397 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00398 $ RESULT( 4 ) ) 00399 CALL ZPTT05( N, NRHS, D, E, B, LDA, X, LDA, XACT, LDA, 00400 $ RWORK, RWORK( NRHS+1 ), RESULT( 5 ) ) 00401 * 00402 * Print information about the tests that did not pass the 00403 * threshold. 00404 * 00405 DO 70 K = 2, 6 00406 IF( RESULT( K ).GE.THRESH ) THEN 00407 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00408 $ CALL ALAHD( NOUT, PATH ) 00409 WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, IMAT, 00410 $ K, RESULT( K ) 00411 NFAIL = NFAIL + 1 00412 END IF 00413 70 CONTINUE 00414 NRUN = NRUN + 5 00415 * 00416 80 CONTINUE 00417 90 CONTINUE 00418 * 00419 *+ TEST 7 00420 * Estimate the reciprocal of the condition number of the 00421 * matrix. 00422 * 00423 100 CONTINUE 00424 SRNAMT = 'ZPTCON' 00425 CALL ZPTCON( N, D( N+1 ), E( N+1 ), ANORM, RCOND, RWORK, 00426 $ INFO ) 00427 * 00428 * Check error code from ZPTCON. 00429 * 00430 IF( INFO.NE.0 ) 00431 $ CALL ALAERH( PATH, 'ZPTCON', INFO, 0, ' ', N, N, -1, -1, 00432 $ -1, IMAT, NFAIL, NERRS, NOUT ) 00433 * 00434 RESULT( 7 ) = DGET06( RCOND, RCONDC ) 00435 * 00436 * Print the test ratio if greater than or equal to THRESH. 00437 * 00438 IF( RESULT( 7 ).GE.THRESH ) THEN 00439 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00440 $ CALL ALAHD( NOUT, PATH ) 00441 WRITE( NOUT, FMT = 9999 )N, IMAT, 7, RESULT( 7 ) 00442 NFAIL = NFAIL + 1 00443 END IF 00444 NRUN = NRUN + 1 00445 110 CONTINUE 00446 120 CONTINUE 00447 * 00448 * Print a summary of the results. 00449 * 00450 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00451 * 00452 9999 FORMAT( ' N =', I5, ', type ', I2, ', test ', I2, ', ratio = ', 00453 $ G12.5 ) 00454 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS =', I3, 00455 $ ', type ', I2, ', test ', I2, ', ratio = ', G12.5 ) 00456 RETURN 00457 * 00458 * End of ZCHKPT 00459 * 00460 END
1.7.3 
