![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZCHKHP 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE ZCHKHP( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR, 00012 * NMAX, A, AFAC, AINV, B, X, XACT, WORK, RWORK, 00013 * IWORK, NOUT ) 00014 * 00015 * .. Scalar Arguments .. 00016 * LOGICAL TSTERR 00017 * INTEGER NMAX, NN, NNS, NOUT 00018 * DOUBLE PRECISION THRESH 00019 * .. 00020 * .. Array Arguments .. 00021 * LOGICAL DOTYPE( * ) 00022 * INTEGER IWORK( * ), NSVAL( * ), NVAL( * ) 00023 * DOUBLE PRECISION RWORK( * ) 00024 * COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), 00025 * $ WORK( * ), X( * ), XACT( * ) 00026 * .. 00027 * 00028 * 00029 *> \par Purpose: 00030 * ============= 00031 *> 00032 *> \verbatim 00033 *> 00034 *> ZCHKHP tests ZHPTRF, -TRI, -TRS, -RFS, and -CON 00035 *> \endverbatim 00036 * 00037 * Arguments: 00038 * ========== 00039 * 00040 *> \param[in] DOTYPE 00041 *> \verbatim 00042 *> DOTYPE is LOGICAL array, dimension (NTYPES) 00043 *> The matrix types to be used for testing. Matrices of type j 00044 *> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 00045 *> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 00046 *> \endverbatim 00047 *> 00048 *> \param[in] NN 00049 *> \verbatim 00050 *> NN is INTEGER 00051 *> The number of values of N contained in the vector NVAL. 00052 *> \endverbatim 00053 *> 00054 *> \param[in] NVAL 00055 *> \verbatim 00056 *> NVAL is INTEGER array, dimension (NN) 00057 *> The values of the matrix dimension N. 00058 *> \endverbatim 00059 *> 00060 *> \param[in] NNS 00061 *> \verbatim 00062 *> NNS is INTEGER 00063 *> The number of values of NRHS contained in the vector NSVAL. 00064 *> \endverbatim 00065 *> 00066 *> \param[in] NSVAL 00067 *> \verbatim 00068 *> NSVAL is INTEGER array, dimension (NNS) 00069 *> The values of the number of right hand sides NRHS. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] THRESH 00073 *> \verbatim 00074 *> THRESH is DOUBLE PRECISION 00075 *> The threshold value for the test ratios. A result is 00076 *> included in the output file if RESULT >= THRESH. To have 00077 *> every test ratio printed, use THRESH = 0. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] TSTERR 00081 *> \verbatim 00082 *> TSTERR is LOGICAL 00083 *> Flag that indicates whether error exits are to be tested. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] NMAX 00087 *> \verbatim 00088 *> NMAX is INTEGER 00089 *> The maximum value permitted for N, used in dimensioning the 00090 *> work arrays. 00091 *> \endverbatim 00092 *> 00093 *> \param[out] A 00094 *> \verbatim 00095 *> A is COMPLEX*16 array, dimension 00096 *> (NMAX*(NMAX+1)/2) 00097 *> \endverbatim 00098 *> 00099 *> \param[out] AFAC 00100 *> \verbatim 00101 *> AFAC is COMPLEX*16 array, dimension 00102 *> (NMAX*(NMAX+1)/2) 00103 *> \endverbatim 00104 *> 00105 *> \param[out] AINV 00106 *> \verbatim 00107 *> AINV is COMPLEX*16 array, dimension 00108 *> (NMAX*(NMAX+1)/2) 00109 *> \endverbatim 00110 *> 00111 *> \param[out] B 00112 *> \verbatim 00113 *> B is COMPLEX*16 array, dimension (NMAX*NSMAX) 00114 *> where NSMAX is the largest entry in NSVAL. 00115 *> \endverbatim 00116 *> 00117 *> \param[out] X 00118 *> \verbatim 00119 *> X is COMPLEX*16 array, dimension (NMAX*NSMAX) 00120 *> \endverbatim 00121 *> 00122 *> \param[out] XACT 00123 *> \verbatim 00124 *> XACT is COMPLEX*16 array, dimension (NMAX*NSMAX) 00125 *> \endverbatim 00126 *> 00127 *> \param[out] WORK 00128 *> \verbatim 00129 *> WORK is COMPLEX*16 array, dimension 00130 *> (NMAX*max(2,NSMAX)) 00131 *> \endverbatim 00132 *> 00133 *> \param[out] RWORK 00134 *> \verbatim 00135 *> RWORK is DOUBLE PRECISION array, 00136 *> dimension (NMAX+2*NSMAX) 00137 *> \endverbatim 00138 *> 00139 *> \param[out] IWORK 00140 *> \verbatim 00141 *> IWORK is INTEGER array, dimension (NMAX) 00142 *> \endverbatim 00143 *> 00144 *> \param[in] NOUT 00145 *> \verbatim 00146 *> NOUT is INTEGER 00147 *> The unit number for output. 00148 *> \endverbatim 00149 * 00150 * Authors: 00151 * ======== 00152 * 00153 *> \author Univ. of Tennessee 00154 *> \author Univ. of California Berkeley 00155 *> \author Univ. of Colorado Denver 00156 *> \author NAG Ltd. 00157 * 00158 *> \date November 2011 00159 * 00160 *> \ingroup complex16_lin 00161 * 00162 * ===================================================================== 00163 SUBROUTINE ZCHKHP( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR, 00164 $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, RWORK, 00165 $ IWORK, NOUT ) 00166 * 00167 * -- LAPACK test routine (version 3.4.0) -- 00168 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00169 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00170 * November 2011 00171 * 00172 * .. Scalar Arguments .. 00173 LOGICAL TSTERR 00174 INTEGER NMAX, NN, NNS, NOUT 00175 DOUBLE PRECISION THRESH 00176 * .. 00177 * .. Array Arguments .. 00178 LOGICAL DOTYPE( * ) 00179 INTEGER IWORK( * ), NSVAL( * ), NVAL( * ) 00180 DOUBLE PRECISION RWORK( * ) 00181 COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), 00182 $ WORK( * ), X( * ), XACT( * ) 00183 * .. 00184 * 00185 * ===================================================================== 00186 * 00187 * .. Parameters .. 00188 DOUBLE PRECISION ZERO 00189 PARAMETER ( ZERO = 0.0D+0 ) 00190 INTEGER NTYPES 00191 PARAMETER ( NTYPES = 10 ) 00192 INTEGER NTESTS 00193 PARAMETER ( NTESTS = 8 ) 00194 * .. 00195 * .. Local Scalars .. 00196 LOGICAL TRFCON, ZEROT 00197 CHARACTER DIST, PACKIT, TYPE, UPLO, XTYPE 00198 CHARACTER*3 PATH 00199 INTEGER I, I1, I2, IMAT, IN, INFO, IOFF, IRHS, IUPLO, 00200 $ IZERO, J, K, KL, KU, LDA, MODE, N, NERRS, 00201 $ NFAIL, NIMAT, NPP, NRHS, NRUN, NT 00202 DOUBLE PRECISION ANORM, CNDNUM, RCOND, RCONDC 00203 * .. 00204 * .. Local Arrays .. 00205 CHARACTER UPLOS( 2 ) 00206 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00207 DOUBLE PRECISION RESULT( NTESTS ) 00208 * .. 00209 * .. External Functions .. 00210 LOGICAL LSAME 00211 DOUBLE PRECISION DGET06, ZLANHP 00212 EXTERNAL LSAME, DGET06, ZLANHP 00213 * .. 00214 * .. External Subroutines .. 00215 EXTERNAL ALAERH, ALAHD, ALASUM, ZCOPY, ZERRSY, ZGET04, 00216 $ ZHPCON, ZHPRFS, ZHPT01, ZHPTRF, ZHPTRI, ZHPTRS, 00217 $ ZLACPY, ZLAIPD, ZLARHS, ZLATB4, ZLATMS, ZPPT02, 00218 $ ZPPT03, ZPPT05 00219 * .. 00220 * .. Intrinsic Functions .. 00221 INTRINSIC MAX, MIN 00222 * .. 00223 * .. Scalars in Common .. 00224 LOGICAL LERR, OK 00225 CHARACTER*32 SRNAMT 00226 INTEGER INFOT, NUNIT 00227 * .. 00228 * .. Common blocks .. 00229 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00230 COMMON / SRNAMC / SRNAMT 00231 * .. 00232 * .. Data statements .. 00233 DATA ISEEDY / 1988, 1989, 1990, 1991 / 00234 DATA UPLOS / 'U', 'L' / 00235 * .. 00236 * .. Executable Statements .. 00237 * 00238 * Initialize constants and the random number seed. 00239 * 00240 PATH( 1: 1 ) = 'Zomplex precision' 00241 PATH( 2: 3 ) = 'HP' 00242 NRUN = 0 00243 NFAIL = 0 00244 NERRS = 0 00245 DO 10 I = 1, 4 00246 ISEED( I ) = ISEEDY( I ) 00247 10 CONTINUE 00248 * 00249 * Test the error exits 00250 * 00251 IF( TSTERR ) 00252 $ CALL ZERRSY( PATH, NOUT ) 00253 INFOT = 0 00254 * 00255 * Do for each value of N in NVAL 00256 * 00257 DO 170 IN = 1, NN 00258 N = NVAL( IN ) 00259 LDA = MAX( N, 1 ) 00260 XTYPE = 'N' 00261 NIMAT = NTYPES 00262 IF( N.LE.0 ) 00263 $ NIMAT = 1 00264 * 00265 IZERO = 0 00266 DO 160 IMAT = 1, NIMAT 00267 * 00268 * Do the tests only if DOTYPE( IMAT ) is true. 00269 * 00270 IF( .NOT.DOTYPE( IMAT ) ) 00271 $ GO TO 160 00272 * 00273 * Skip types 3, 4, 5, or 6 if the matrix size is too small. 00274 * 00275 ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 00276 IF( ZEROT .AND. N.LT.IMAT-2 ) 00277 $ GO TO 160 00278 * 00279 * Do first for UPLO = 'U', then for UPLO = 'L' 00280 * 00281 DO 150 IUPLO = 1, 2 00282 UPLO = UPLOS( IUPLO ) 00283 IF( LSAME( UPLO, 'U' ) ) THEN 00284 PACKIT = 'C' 00285 ELSE 00286 PACKIT = 'R' 00287 END IF 00288 * 00289 * Set up parameters with ZLATB4 and generate a test matrix 00290 * with ZLATMS. 00291 * 00292 CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00293 $ CNDNUM, DIST ) 00294 * 00295 SRNAMT = 'ZLATMS' 00296 CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, 00297 $ CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK, 00298 $ INFO ) 00299 * 00300 * Check error code from ZLATMS. 00301 * 00302 IF( INFO.NE.0 ) THEN 00303 CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N, -1, 00304 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00305 GO TO 150 00306 END IF 00307 * 00308 * For types 3-6, zero one or more rows and columns of 00309 * the matrix to test that INFO is returned correctly. 00310 * 00311 IF( ZEROT ) THEN 00312 IF( IMAT.EQ.3 ) THEN 00313 IZERO = 1 00314 ELSE IF( IMAT.EQ.4 ) THEN 00315 IZERO = N 00316 ELSE 00317 IZERO = N / 2 + 1 00318 END IF 00319 * 00320 IF( IMAT.LT.6 ) THEN 00321 * 00322 * Set row and column IZERO to zero. 00323 * 00324 IF( IUPLO.EQ.1 ) THEN 00325 IOFF = ( IZERO-1 )*IZERO / 2 00326 DO 20 I = 1, IZERO - 1 00327 A( IOFF+I ) = ZERO 00328 20 CONTINUE 00329 IOFF = IOFF + IZERO 00330 DO 30 I = IZERO, N 00331 A( IOFF ) = ZERO 00332 IOFF = IOFF + I 00333 30 CONTINUE 00334 ELSE 00335 IOFF = IZERO 00336 DO 40 I = 1, IZERO - 1 00337 A( IOFF ) = ZERO 00338 IOFF = IOFF + N - I 00339 40 CONTINUE 00340 IOFF = IOFF - IZERO 00341 DO 50 I = IZERO, N 00342 A( IOFF+I ) = ZERO 00343 50 CONTINUE 00344 END IF 00345 ELSE 00346 IOFF = 0 00347 IF( IUPLO.EQ.1 ) THEN 00348 * 00349 * Set the first IZERO rows and columns to zero. 00350 * 00351 DO 70 J = 1, N 00352 I2 = MIN( J, IZERO ) 00353 DO 60 I = 1, I2 00354 A( IOFF+I ) = ZERO 00355 60 CONTINUE 00356 IOFF = IOFF + J 00357 70 CONTINUE 00358 ELSE 00359 * 00360 * Set the last IZERO rows and columns to zero. 00361 * 00362 DO 90 J = 1, N 00363 I1 = MAX( J, IZERO ) 00364 DO 80 I = I1, N 00365 A( IOFF+I ) = ZERO 00366 80 CONTINUE 00367 IOFF = IOFF + N - J 00368 90 CONTINUE 00369 END IF 00370 END IF 00371 ELSE 00372 IZERO = 0 00373 END IF 00374 * 00375 * Set the imaginary part of the diagonals. 00376 * 00377 IF( IUPLO.EQ.1 ) THEN 00378 CALL ZLAIPD( N, A, 2, 1 ) 00379 ELSE 00380 CALL ZLAIPD( N, A, N, -1 ) 00381 END IF 00382 * 00383 * Compute the L*D*L' or U*D*U' factorization of the matrix. 00384 * 00385 NPP = N*( N+1 ) / 2 00386 CALL ZCOPY( NPP, A, 1, AFAC, 1 ) 00387 SRNAMT = 'ZHPTRF' 00388 CALL ZHPTRF( UPLO, N, AFAC, IWORK, INFO ) 00389 * 00390 * Adjust the expected value of INFO to account for 00391 * pivoting. 00392 * 00393 K = IZERO 00394 IF( K.GT.0 ) THEN 00395 100 CONTINUE 00396 IF( IWORK( K ).LT.0 ) THEN 00397 IF( IWORK( K ).NE.-K ) THEN 00398 K = -IWORK( K ) 00399 GO TO 100 00400 END IF 00401 ELSE IF( IWORK( K ).NE.K ) THEN 00402 K = IWORK( K ) 00403 GO TO 100 00404 END IF 00405 END IF 00406 * 00407 * Check error code from ZHPTRF. 00408 * 00409 IF( INFO.NE.K ) 00410 $ CALL ALAERH( PATH, 'ZHPTRF', INFO, K, UPLO, N, N, -1, 00411 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00412 IF( INFO.NE.0 ) THEN 00413 TRFCON = .TRUE. 00414 ELSE 00415 TRFCON = .FALSE. 00416 END IF 00417 * 00418 *+ TEST 1 00419 * Reconstruct matrix from factors and compute residual. 00420 * 00421 CALL ZHPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA, RWORK, 00422 $ RESULT( 1 ) ) 00423 NT = 1 00424 * 00425 *+ TEST 2 00426 * Form the inverse and compute the residual. 00427 * 00428 IF( .NOT.TRFCON ) THEN 00429 CALL ZCOPY( NPP, AFAC, 1, AINV, 1 ) 00430 SRNAMT = 'ZHPTRI' 00431 CALL ZHPTRI( UPLO, N, AINV, IWORK, WORK, INFO ) 00432 * 00433 * Check error code from ZHPTRI. 00434 * 00435 IF( INFO.NE.0 ) 00436 $ CALL ALAERH( PATH, 'ZHPTRI', INFO, 0, UPLO, N, N, 00437 $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00438 * 00439 CALL ZPPT03( UPLO, N, A, AINV, WORK, LDA, RWORK, 00440 $ RCONDC, RESULT( 2 ) ) 00441 NT = 2 00442 END IF 00443 * 00444 * Print information about the tests that did not pass 00445 * the threshold. 00446 * 00447 DO 110 K = 1, NT 00448 IF( RESULT( K ).GE.THRESH ) THEN 00449 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00450 $ CALL ALAHD( NOUT, PATH ) 00451 WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, K, 00452 $ RESULT( K ) 00453 NFAIL = NFAIL + 1 00454 END IF 00455 110 CONTINUE 00456 NRUN = NRUN + NT 00457 * 00458 * Do only the condition estimate if INFO is not 0. 00459 * 00460 IF( TRFCON ) THEN 00461 RCONDC = ZERO 00462 GO TO 140 00463 END IF 00464 * 00465 DO 130 IRHS = 1, NNS 00466 NRHS = NSVAL( IRHS ) 00467 * 00468 *+ TEST 3 00469 * Solve and compute residual for A * X = B. 00470 * 00471 SRNAMT = 'ZLARHS' 00472 CALL ZLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, 00473 $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED, 00474 $ INFO ) 00475 XTYPE = 'C' 00476 CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00477 * 00478 SRNAMT = 'ZHPTRS' 00479 CALL ZHPTRS( UPLO, N, NRHS, AFAC, IWORK, X, LDA, 00480 $ INFO ) 00481 * 00482 * Check error code from ZHPTRS. 00483 * 00484 IF( INFO.NE.0 ) 00485 $ CALL ALAERH( PATH, 'ZHPTRS', INFO, 0, UPLO, N, N, 00486 $ -1, -1, NRHS, IMAT, NFAIL, NERRS, 00487 $ NOUT ) 00488 * 00489 CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) 00490 CALL ZPPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA, 00491 $ RWORK, RESULT( 3 ) ) 00492 * 00493 *+ TEST 4 00494 * Check solution from generated exact solution. 00495 * 00496 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00497 $ RESULT( 4 ) ) 00498 * 00499 *+ TESTS 5, 6, and 7 00500 * Use iterative refinement to improve the solution. 00501 * 00502 SRNAMT = 'ZHPRFS' 00503 CALL ZHPRFS( UPLO, N, NRHS, A, AFAC, IWORK, B, LDA, X, 00504 $ LDA, RWORK, RWORK( NRHS+1 ), WORK, 00505 $ RWORK( 2*NRHS+1 ), INFO ) 00506 * 00507 * Check error code from ZHPRFS. 00508 * 00509 IF( INFO.NE.0 ) 00510 $ CALL ALAERH( PATH, 'ZHPRFS', INFO, 0, UPLO, N, N, 00511 $ -1, -1, NRHS, IMAT, NFAIL, NERRS, 00512 $ NOUT ) 00513 * 00514 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00515 $ RESULT( 5 ) ) 00516 CALL ZPPT05( UPLO, N, NRHS, A, B, LDA, X, LDA, XACT, 00517 $ LDA, RWORK, RWORK( NRHS+1 ), 00518 $ RESULT( 6 ) ) 00519 * 00520 * Print information about the tests that did not pass 00521 * the threshold. 00522 * 00523 DO 120 K = 3, 7 00524 IF( RESULT( K ).GE.THRESH ) THEN 00525 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00526 $ CALL ALAHD( NOUT, PATH ) 00527 WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, IMAT, 00528 $ K, RESULT( K ) 00529 NFAIL = NFAIL + 1 00530 END IF 00531 120 CONTINUE 00532 NRUN = NRUN + 5 00533 130 CONTINUE 00534 * 00535 *+ TEST 8 00536 * Get an estimate of RCOND = 1/CNDNUM. 00537 * 00538 140 CONTINUE 00539 ANORM = ZLANHP( '1', UPLO, N, A, RWORK ) 00540 SRNAMT = 'ZHPCON' 00541 CALL ZHPCON( UPLO, N, AFAC, IWORK, ANORM, RCOND, WORK, 00542 $ INFO ) 00543 * 00544 * Check error code from ZHPCON. 00545 * 00546 IF( INFO.NE.0 ) 00547 $ CALL ALAERH( PATH, 'ZHPCON', INFO, 0, UPLO, N, N, -1, 00548 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00549 * 00550 RESULT( 8 ) = DGET06( RCOND, RCONDC ) 00551 * 00552 * Print the test ratio if it is .GE. THRESH. 00553 * 00554 IF( RESULT( 8 ).GE.THRESH ) THEN 00555 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00556 $ CALL ALAHD( NOUT, PATH ) 00557 WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, 8, 00558 $ RESULT( 8 ) 00559 NFAIL = NFAIL + 1 00560 END IF 00561 NRUN = NRUN + 1 00562 150 CONTINUE 00563 160 CONTINUE 00564 170 CONTINUE 00565 * 00566 * Print a summary of the results. 00567 * 00568 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00569 * 00570 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', type ', I2, ', test ', 00571 $ I2, ', ratio =', G12.5 ) 00572 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', 00573 $ I2, ', test(', I2, ') =', G12.5 ) 00574 RETURN 00575 * 00576 * End of ZCHKHP 00577 * 00578 END