 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CCHKHP( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR, 00002 $ NMAX, A, AFAC, AINV, 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 RWORK( * ) 00018 COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ), 00019 $ WORK( * ), X( * ), XACT( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * CCHKHP tests CHPTRF, -TRI, -TRS, -RFS, and -CON 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 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) REAL 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 maximum value permitted for N, used in dimensioning the 00057 * work arrays. 00058 * 00059 * A (workspace) COMPLEX array, dimension 00060 * (NMAX*(NMAX+1)/2) 00061 * 00062 * AFAC (workspace) COMPLEX array, dimension 00063 * (NMAX*(NMAX+1)/2) 00064 * 00065 * AINV (workspace) COMPLEX array, dimension 00066 * (NMAX*(NMAX+1)/2) 00067 * 00068 * B (workspace) COMPLEX array, dimension (NMAX*NSMAX) 00069 * where NSMAX is the largest entry in NSVAL. 00070 * 00071 * X (workspace) COMPLEX array, dimension (NMAX*NSMAX) 00072 * 00073 * XACT (workspace) COMPLEX array, dimension (NMAX*NSMAX) 00074 * 00075 * WORK (workspace) COMPLEX array, dimension 00076 * (NMAX*max(2,NSMAX)) 00077 * 00078 * RWORK (workspace) REAL array, 00079 * dimension (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 REAL ZERO 00090 PARAMETER ( ZERO = 0.0E+0 ) 00091 INTEGER NTYPES 00092 PARAMETER ( NTYPES = 10 ) 00093 INTEGER NTESTS 00094 PARAMETER ( NTESTS = 8 ) 00095 * .. 00096 * .. Local Scalars .. 00097 LOGICAL TRFCON, ZEROT 00098 CHARACTER DIST, PACKIT, TYPE, UPLO, XTYPE 00099 CHARACTER*3 PATH 00100 INTEGER I, I1, I2, IMAT, IN, INFO, IOFF, IRHS, IUPLO, 00101 $ IZERO, J, K, KL, KU, LDA, MODE, N, NERRS, 00102 $ NFAIL, NIMAT, NPP, NRHS, NRUN, NT 00103 REAL ANORM, CNDNUM, RCOND, RCONDC 00104 * .. 00105 * .. Local Arrays .. 00106 CHARACTER UPLOS( 2 ) 00107 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00108 REAL RESULT( NTESTS ) 00109 * .. 00110 * .. External Functions .. 00111 LOGICAL LSAME 00112 REAL CLANHP, SGET06 00113 EXTERNAL LSAME, CLANHP, SGET06 00114 * .. 00115 * .. External Subroutines .. 00116 EXTERNAL ALAERH, ALAHD, ALASUM, CCOPY, CERRSY, CGET04, 00117 $ CHPCON, CHPRFS, CHPT01, CHPTRF, CHPTRI, CHPTRS, 00118 $ CLACPY, CLAIPD, CLARHS, CLATB4, CLATMS, CPPT02, 00119 $ CPPT03, CPPT05 00120 * .. 00121 * .. Intrinsic Functions .. 00122 INTRINSIC MAX, MIN 00123 * .. 00124 * .. Scalars in Common .. 00125 LOGICAL LERR, OK 00126 CHARACTER*32 SRNAMT 00127 INTEGER INFOT, NUNIT 00128 * .. 00129 * .. Common blocks .. 00130 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00131 COMMON / SRNAMC / SRNAMT 00132 * .. 00133 * .. Data statements .. 00134 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00135 DATA UPLOS / 'U', 'L' / 00136 * .. 00137 * .. Executable Statements .. 00138 * 00139 * Initialize constants and the random number seed. 00140 * 00141 PATH( 1: 1 ) = 'Complex precision' 00142 PATH( 2: 3 ) = 'HP' 00143 NRUN = 0 00144 NFAIL = 0 00145 NERRS = 0 00146 DO 10 I = 1, 4 00147 ISEED( I ) = ISEEDY( I ) 00148 10 CONTINUE 00149 * 00150 * Test the error exits 00151 * 00152 IF( TSTERR ) 00153 $ CALL CERRSY( PATH, NOUT ) 00154 INFOT = 0 00155 * 00156 * Do for each value of N in NVAL 00157 * 00158 DO 170 IN = 1, NN 00159 N = NVAL( IN ) 00160 LDA = MAX( N, 1 ) 00161 XTYPE = 'N' 00162 NIMAT = NTYPES 00163 IF( N.LE.0 ) 00164 $ NIMAT = 1 00165 * 00166 IZERO = 0 00167 DO 160 IMAT = 1, NIMAT 00168 * 00169 * Do the tests only if DOTYPE( IMAT ) is true. 00170 * 00171 IF( .NOT.DOTYPE( IMAT ) ) 00172 $ GO TO 160 00173 * 00174 * Skip types 3, 4, 5, or 6 if the matrix size is too small. 00175 * 00176 ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 00177 IF( ZEROT .AND. N.LT.IMAT-2 ) 00178 $ GO TO 160 00179 * 00180 * Do first for UPLO = 'U', then for UPLO = 'L' 00181 * 00182 DO 150 IUPLO = 1, 2 00183 UPLO = UPLOS( IUPLO ) 00184 IF( LSAME( UPLO, 'U' ) ) THEN 00185 PACKIT = 'C' 00186 ELSE 00187 PACKIT = 'R' 00188 END IF 00189 * 00190 * Set up parameters with CLATB4 and generate a test matrix 00191 * with CLATMS. 00192 * 00193 CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00194 $ CNDNUM, DIST ) 00195 * 00196 SRNAMT = 'CLATMS' 00197 CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 00198 $ CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK, 00199 $ INFO ) 00200 * 00201 * Check error code from CLATMS. 00202 * 00203 IF( INFO.NE.0 ) THEN 00204 CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N, -1, 00205 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00206 GO TO 150 00207 END IF 00208 * 00209 * For types 3-6, zero one or more rows and columns of 00210 * the matrix to test that INFO is returned correctly. 00211 * 00212 IF( ZEROT ) THEN 00213 IF( IMAT.EQ.3 ) THEN 00214 IZERO = 1 00215 ELSE IF( IMAT.EQ.4 ) THEN 00216 IZERO = N 00217 ELSE 00218 IZERO = N / 2 + 1 00219 END IF 00220 * 00221 IF( IMAT.LT.6 ) THEN 00222 * 00223 * Set row and column IZERO to zero. 00224 * 00225 IF( IUPLO.EQ.1 ) THEN 00226 IOFF = ( IZERO-1 )*IZERO / 2 00227 DO 20 I = 1, IZERO - 1 00228 A( IOFF+I ) = ZERO 00229 20 CONTINUE 00230 IOFF = IOFF + IZERO 00231 DO 30 I = IZERO, N 00232 A( IOFF ) = ZERO 00233 IOFF = IOFF + I 00234 30 CONTINUE 00235 ELSE 00236 IOFF = IZERO 00237 DO 40 I = 1, IZERO - 1 00238 A( IOFF ) = ZERO 00239 IOFF = IOFF + N - I 00240 40 CONTINUE 00241 IOFF = IOFF - IZERO 00242 DO 50 I = IZERO, N 00243 A( IOFF+I ) = ZERO 00244 50 CONTINUE 00245 END IF 00246 ELSE 00247 IOFF = 0 00248 IF( IUPLO.EQ.1 ) THEN 00249 * 00250 * Set the first IZERO rows and columns to zero. 00251 * 00252 DO 70 J = 1, N 00253 I2 = MIN( J, IZERO ) 00254 DO 60 I = 1, I2 00255 A( IOFF+I ) = ZERO 00256 60 CONTINUE 00257 IOFF = IOFF + J 00258 70 CONTINUE 00259 ELSE 00260 * 00261 * Set the last IZERO rows and columns to zero. 00262 * 00263 DO 90 J = 1, N 00264 I1 = MAX( J, IZERO ) 00265 DO 80 I = I1, N 00266 A( IOFF+I ) = ZERO 00267 80 CONTINUE 00268 IOFF = IOFF + N - J 00269 90 CONTINUE 00270 END IF 00271 END IF 00272 ELSE 00273 IZERO = 0 00274 END IF 00275 * 00276 * Set the imaginary part of the diagonals. 00277 * 00278 IF( IUPLO.EQ.1 ) THEN 00279 CALL CLAIPD( N, A, 2, 1 ) 00280 ELSE 00281 CALL CLAIPD( N, A, N, -1 ) 00282 END IF 00283 * 00284 * Compute the L*D*L' or U*D*U' factorization of the matrix. 00285 * 00286 NPP = N*( N+1 ) / 2 00287 CALL CCOPY( NPP, A, 1, AFAC, 1 ) 00288 SRNAMT = 'CHPTRF' 00289 CALL CHPTRF( UPLO, N, AFAC, IWORK, INFO ) 00290 * 00291 * Adjust the expected value of INFO to account for 00292 * pivoting. 00293 * 00294 K = IZERO 00295 IF( K.GT.0 ) THEN 00296 100 CONTINUE 00297 IF( IWORK( K ).LT.0 ) THEN 00298 IF( IWORK( K ).NE.-K ) THEN 00299 K = -IWORK( K ) 00300 GO TO 100 00301 END IF 00302 ELSE IF( IWORK( K ).NE.K ) THEN 00303 K = IWORK( K ) 00304 GO TO 100 00305 END IF 00306 END IF 00307 * 00308 * Check error code from CHPTRF. 00309 * 00310 IF( INFO.NE.K ) 00311 $ CALL ALAERH( PATH, 'CHPTRF', INFO, K, UPLO, N, N, -1, 00312 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00313 IF( INFO.NE.0 ) THEN 00314 TRFCON = .TRUE. 00315 ELSE 00316 TRFCON = .FALSE. 00317 END IF 00318 * 00319 *+ TEST 1 00320 * Reconstruct matrix from factors and compute residual. 00321 * 00322 CALL CHPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA, RWORK, 00323 $ RESULT( 1 ) ) 00324 NT = 1 00325 * 00326 *+ TEST 2 00327 * Form the inverse and compute the residual. 00328 * 00329 IF( .NOT.TRFCON ) THEN 00330 CALL CCOPY( NPP, AFAC, 1, AINV, 1 ) 00331 SRNAMT = 'CHPTRI' 00332 CALL CHPTRI( UPLO, N, AINV, IWORK, WORK, INFO ) 00333 * 00334 * Check error code from CHPTRI. 00335 * 00336 IF( INFO.NE.0 ) 00337 $ CALL ALAERH( PATH, 'CHPTRI', INFO, 0, UPLO, N, N, 00338 $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00339 * 00340 CALL CPPT03( UPLO, N, A, AINV, WORK, LDA, RWORK, 00341 $ RCONDC, RESULT( 2 ) ) 00342 NT = 2 00343 END IF 00344 * 00345 * Print information about the tests that did not pass 00346 * the threshold. 00347 * 00348 DO 110 K = 1, NT 00349 IF( RESULT( K ).GE.THRESH ) THEN 00350 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00351 $ CALL ALAHD( NOUT, PATH ) 00352 WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, K, 00353 $ RESULT( K ) 00354 NFAIL = NFAIL + 1 00355 END IF 00356 110 CONTINUE 00357 NRUN = NRUN + NT 00358 * 00359 * Do only the condition estimate if INFO is not 0. 00360 * 00361 IF( TRFCON ) THEN 00362 RCONDC = ZERO 00363 GO TO 140 00364 END IF 00365 * 00366 DO 130 IRHS = 1, NNS 00367 NRHS = NSVAL( IRHS ) 00368 * 00369 *+ TEST 3 00370 * Solve and compute residual for A * X = B. 00371 * 00372 SRNAMT = 'CLARHS' 00373 CALL CLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, 00374 $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED, 00375 $ INFO ) 00376 XTYPE = 'C' 00377 CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00378 * 00379 SRNAMT = 'CHPTRS' 00380 CALL CHPTRS( UPLO, N, NRHS, AFAC, IWORK, X, LDA, 00381 $ INFO ) 00382 * 00383 * Check error code from CHPTRS. 00384 * 00385 IF( INFO.NE.0 ) 00386 $ CALL ALAERH( PATH, 'CHPTRS', INFO, 0, UPLO, N, N, 00387 $ -1, -1, NRHS, IMAT, NFAIL, NERRS, 00388 $ NOUT ) 00389 * 00390 CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) 00391 CALL CPPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA, 00392 $ RWORK, RESULT( 3 ) ) 00393 * 00394 *+ TEST 4 00395 * Check solution from generated exact solution. 00396 * 00397 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00398 $ RESULT( 4 ) ) 00399 * 00400 *+ TESTS 5, 6, and 7 00401 * Use iterative refinement to improve the solution. 00402 * 00403 SRNAMT = 'CHPRFS' 00404 CALL CHPRFS( UPLO, N, NRHS, A, AFAC, IWORK, B, LDA, X, 00405 $ LDA, RWORK, RWORK( NRHS+1 ), WORK, 00406 $ RWORK( 2*NRHS+1 ), INFO ) 00407 * 00408 * Check error code from CHPRFS. 00409 * 00410 IF( INFO.NE.0 ) 00411 $ CALL ALAERH( PATH, 'CHPRFS', INFO, 0, UPLO, N, N, 00412 $ -1, -1, NRHS, IMAT, NFAIL, NERRS, 00413 $ NOUT ) 00414 * 00415 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00416 $ RESULT( 5 ) ) 00417 CALL CPPT05( UPLO, N, NRHS, A, B, LDA, X, LDA, XACT, 00418 $ LDA, RWORK, RWORK( NRHS+1 ), 00419 $ RESULT( 6 ) ) 00420 * 00421 * Print information about the tests that did not pass 00422 * the threshold. 00423 * 00424 DO 120 K = 3, 7 00425 IF( RESULT( K ).GE.THRESH ) THEN 00426 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00427 $ CALL ALAHD( NOUT, PATH ) 00428 WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, IMAT, 00429 $ K, RESULT( K ) 00430 NFAIL = NFAIL + 1 00431 END IF 00432 120 CONTINUE 00433 NRUN = NRUN + 5 00434 130 CONTINUE 00435 * 00436 *+ TEST 8 00437 * Get an estimate of RCOND = 1/CNDNUM. 00438 * 00439 140 CONTINUE 00440 ANORM = CLANHP( '1', UPLO, N, A, RWORK ) 00441 SRNAMT = 'CHPCON' 00442 CALL CHPCON( UPLO, N, AFAC, IWORK, ANORM, RCOND, WORK, 00443 $ INFO ) 00444 * 00445 * Check error code from CHPCON. 00446 * 00447 IF( INFO.NE.0 ) 00448 $ CALL ALAERH( PATH, 'CHPCON', INFO, 0, UPLO, N, N, -1, 00449 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00450 * 00451 RESULT( 8 ) = SGET06( RCOND, RCONDC ) 00452 * 00453 * Print the test ratio if it is .GE. THRESH. 00454 * 00455 IF( RESULT( 8 ).GE.THRESH ) THEN 00456 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00457 $ CALL ALAHD( NOUT, PATH ) 00458 WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, 8, 00459 $ RESULT( 8 ) 00460 NFAIL = NFAIL + 1 00461 END IF 00462 NRUN = NRUN + 1 00463 150 CONTINUE 00464 160 CONTINUE 00465 170 CONTINUE 00466 * 00467 * Print a summary of the results. 00468 * 00469 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00470 * 00471 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', type ', I2, ', test ', 00472 $ I2, ', ratio =', G12.5 ) 00473 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', 00474 $ I2, ', test(', I2, ') =', G12.5 ) 00475 RETURN 00476 * 00477 * End of CCHKHP 00478 * 00479 END
1.7.3 
