![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CCHKPT 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 CCHKPT( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR, 00012 * A, D, E, B, X, XACT, WORK, RWORK, NOUT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * LOGICAL TSTERR 00016 * INTEGER NN, NNS, NOUT 00017 * REAL THRESH 00018 * .. 00019 * .. Array Arguments .. 00020 * LOGICAL DOTYPE( * ) 00021 * INTEGER NSVAL( * ), NVAL( * ) 00022 * REAL D( * ), RWORK( * ) 00023 * COMPLEX A( * ), B( * ), E( * ), WORK( * ), X( * ), 00024 * $ XACT( * ) 00025 * .. 00026 * 00027 * 00028 *> \par Purpose: 00029 * ============= 00030 *> 00031 *> \verbatim 00032 *> 00033 *> CCHKPT tests CPTTRF, -TRS, -RFS, and -CON 00034 *> \endverbatim 00035 * 00036 * Arguments: 00037 * ========== 00038 * 00039 *> \param[in] DOTYPE 00040 *> \verbatim 00041 *> DOTYPE is LOGICAL array, dimension (NTYPES) 00042 *> The matrix types to be used for testing. Matrices of type j 00043 *> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 00044 *> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 00045 *> \endverbatim 00046 *> 00047 *> \param[in] NN 00048 *> \verbatim 00049 *> NN is INTEGER 00050 *> The number of values of N contained in the vector NVAL. 00051 *> \endverbatim 00052 *> 00053 *> \param[in] NVAL 00054 *> \verbatim 00055 *> NVAL is INTEGER array, dimension (NN) 00056 *> The values of the matrix dimension N. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] NNS 00060 *> \verbatim 00061 *> NNS is INTEGER 00062 *> The number of values of NRHS contained in the vector NSVAL. 00063 *> \endverbatim 00064 *> 00065 *> \param[in] NSVAL 00066 *> \verbatim 00067 *> NSVAL is INTEGER array, dimension (NNS) 00068 *> The values of the number of right hand sides NRHS. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] THRESH 00072 *> \verbatim 00073 *> THRESH is REAL 00074 *> The threshold value for the test ratios. A result is 00075 *> included in the output file if RESULT >= THRESH. To have 00076 *> every test ratio printed, use THRESH = 0. 00077 *> \endverbatim 00078 *> 00079 *> \param[in] TSTERR 00080 *> \verbatim 00081 *> TSTERR is LOGICAL 00082 *> Flag that indicates whether error exits are to be tested. 00083 *> \endverbatim 00084 *> 00085 *> \param[out] A 00086 *> \verbatim 00087 *> A is COMPLEX array, dimension (NMAX*2) 00088 *> \endverbatim 00089 *> 00090 *> \param[out] D 00091 *> \verbatim 00092 *> D is REAL array, dimension (NMAX*2) 00093 *> \endverbatim 00094 *> 00095 *> \param[out] E 00096 *> \verbatim 00097 *> E is COMPLEX array, dimension (NMAX*2) 00098 *> \endverbatim 00099 *> 00100 *> \param[out] B 00101 *> \verbatim 00102 *> B is COMPLEX array, dimension (NMAX*NSMAX) 00103 *> where NSMAX is the largest entry in NSVAL. 00104 *> \endverbatim 00105 *> 00106 *> \param[out] X 00107 *> \verbatim 00108 *> X is COMPLEX array, dimension (NMAX*NSMAX) 00109 *> \endverbatim 00110 *> 00111 *> \param[out] XACT 00112 *> \verbatim 00113 *> XACT is COMPLEX array, dimension (NMAX*NSMAX) 00114 *> \endverbatim 00115 *> 00116 *> \param[out] WORK 00117 *> \verbatim 00118 *> WORK is COMPLEX array, dimension 00119 *> (NMAX*max(3,NSMAX)) 00120 *> \endverbatim 00121 *> 00122 *> \param[out] RWORK 00123 *> \verbatim 00124 *> RWORK is REAL array, dimension 00125 *> (max(NMAX,2*NSMAX)) 00126 *> \endverbatim 00127 *> 00128 *> \param[in] NOUT 00129 *> \verbatim 00130 *> NOUT is INTEGER 00131 *> The unit number for output. 00132 *> \endverbatim 00133 * 00134 * Authors: 00135 * ======== 00136 * 00137 *> \author Univ. of Tennessee 00138 *> \author Univ. of California Berkeley 00139 *> \author Univ. of Colorado Denver 00140 *> \author NAG Ltd. 00141 * 00142 *> \date November 2011 00143 * 00144 *> \ingroup complex_lin 00145 * 00146 * ===================================================================== 00147 SUBROUTINE CCHKPT( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR, 00148 $ A, D, E, B, X, XACT, WORK, RWORK, NOUT ) 00149 * 00150 * -- LAPACK test routine (version 3.4.0) -- 00151 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00152 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00153 * November 2011 00154 * 00155 * .. Scalar Arguments .. 00156 LOGICAL TSTERR 00157 INTEGER NN, NNS, NOUT 00158 REAL THRESH 00159 * .. 00160 * .. Array Arguments .. 00161 LOGICAL DOTYPE( * ) 00162 INTEGER NSVAL( * ), NVAL( * ) 00163 REAL D( * ), RWORK( * ) 00164 COMPLEX A( * ), B( * ), E( * ), WORK( * ), X( * ), 00165 $ XACT( * ) 00166 * .. 00167 * 00168 * ===================================================================== 00169 * 00170 * .. Parameters .. 00171 REAL ONE, ZERO 00172 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00173 INTEGER NTYPES 00174 PARAMETER ( NTYPES = 12 ) 00175 INTEGER NTESTS 00176 PARAMETER ( NTESTS = 7 ) 00177 * .. 00178 * .. Local Scalars .. 00179 LOGICAL ZEROT 00180 CHARACTER DIST, TYPE, UPLO 00181 CHARACTER*3 PATH 00182 INTEGER I, IA, IMAT, IN, INFO, IRHS, IUPLO, IX, IZERO, 00183 $ J, K, KL, KU, LDA, MODE, N, NERRS, NFAIL, 00184 $ NIMAT, NRHS, NRUN 00185 REAL AINVNM, ANORM, COND, DMAX, RCOND, RCONDC 00186 * .. 00187 * .. Local Arrays .. 00188 CHARACTER UPLOS( 2 ) 00189 INTEGER ISEED( 4 ), ISEEDY( 4 ) 00190 REAL RESULT( NTESTS ) 00191 COMPLEX Z( 3 ) 00192 * .. 00193 * .. External Functions .. 00194 INTEGER ISAMAX 00195 REAL CLANHT, SCASUM, SGET06 00196 EXTERNAL ISAMAX, CLANHT, SCASUM, SGET06 00197 * .. 00198 * .. External Subroutines .. 00199 EXTERNAL ALAERH, ALAHD, ALASUM, CCOPY, CERRGT, CGET04, 00200 $ CLACPY, CLAPTM, CLARNV, CLATB4, CLATMS, CPTCON, 00201 $ CPTRFS, CPTT01, CPTT02, CPTT05, CPTTRF, CPTTRS, 00202 $ CSSCAL, SCOPY, SLARNV, SSCAL 00203 * .. 00204 * .. Intrinsic Functions .. 00205 INTRINSIC ABS, MAX, REAL 00206 * .. 00207 * .. Scalars in Common .. 00208 LOGICAL LERR, OK 00209 CHARACTER*32 SRNAMT 00210 INTEGER INFOT, NUNIT 00211 * .. 00212 * .. Common blocks .. 00213 COMMON / INFOC / INFOT, NUNIT, OK, LERR 00214 COMMON / SRNAMC / SRNAMT 00215 * .. 00216 * .. Data statements .. 00217 DATA ISEEDY / 0, 0, 0, 1 / , UPLOS / 'U', 'L' / 00218 * .. 00219 * .. Executable Statements .. 00220 * 00221 PATH( 1: 1 ) = 'Complex precision' 00222 PATH( 2: 3 ) = 'PT' 00223 NRUN = 0 00224 NFAIL = 0 00225 NERRS = 0 00226 DO 10 I = 1, 4 00227 ISEED( I ) = ISEEDY( I ) 00228 10 CONTINUE 00229 * 00230 * Test the error exits 00231 * 00232 IF( TSTERR ) 00233 $ CALL CERRGT( PATH, NOUT ) 00234 INFOT = 0 00235 * 00236 DO 120 IN = 1, NN 00237 * 00238 * Do for each value of N in NVAL. 00239 * 00240 N = NVAL( IN ) 00241 LDA = MAX( 1, N ) 00242 NIMAT = NTYPES 00243 IF( N.LE.0 ) 00244 $ NIMAT = 1 00245 * 00246 DO 110 IMAT = 1, NIMAT 00247 * 00248 * Do the tests only if DOTYPE( IMAT ) is true. 00249 * 00250 IF( N.GT.0 .AND. .NOT.DOTYPE( IMAT ) ) 00251 $ GO TO 110 00252 * 00253 * Set up parameters with CLATB4. 00254 * 00255 CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00256 $ COND, DIST ) 00257 * 00258 ZEROT = IMAT.GE.8 .AND. IMAT.LE.10 00259 IF( IMAT.LE.6 ) THEN 00260 * 00261 * Type 1-6: generate a Hermitian tridiagonal matrix of 00262 * known condition number in lower triangular band storage. 00263 * 00264 SRNAMT = 'CLATMS' 00265 CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, COND, 00266 $ ANORM, KL, KU, 'B', A, 2, WORK, INFO ) 00267 * 00268 * Check the error code from CLATMS. 00269 * 00270 IF( INFO.NE.0 ) THEN 00271 CALL ALAERH( PATH, 'CLATMS', INFO, 0, ' ', N, N, KL, 00272 $ KU, -1, IMAT, NFAIL, NERRS, NOUT ) 00273 GO TO 110 00274 END IF 00275 IZERO = 0 00276 * 00277 * Copy the matrix to D and E. 00278 * 00279 IA = 1 00280 DO 20 I = 1, N - 1 00281 D( I ) = REAL( A( IA ) ) 00282 E( I ) = A( IA+1 ) 00283 IA = IA + 2 00284 20 CONTINUE 00285 IF( N.GT.0 ) 00286 $ D( N ) = REAL( A( IA ) ) 00287 ELSE 00288 * 00289 * Type 7-12: generate a diagonally dominant matrix with 00290 * unknown condition number in the vectors D and E. 00291 * 00292 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 7 ) ) THEN 00293 * 00294 * Let E be complex, D real, with values from [-1,1]. 00295 * 00296 CALL SLARNV( 2, ISEED, N, D ) 00297 CALL CLARNV( 2, ISEED, N-1, E ) 00298 * 00299 * Make the tridiagonal matrix diagonally dominant. 00300 * 00301 IF( N.EQ.1 ) THEN 00302 D( 1 ) = ABS( D( 1 ) ) 00303 ELSE 00304 D( 1 ) = ABS( D( 1 ) ) + ABS( E( 1 ) ) 00305 D( N ) = ABS( D( N ) ) + ABS( E( N-1 ) ) 00306 DO 30 I = 2, N - 1 00307 D( I ) = ABS( D( I ) ) + ABS( E( I ) ) + 00308 $ ABS( E( I-1 ) ) 00309 30 CONTINUE 00310 END IF 00311 * 00312 * Scale D and E so the maximum element is ANORM. 00313 * 00314 IX = ISAMAX( N, D, 1 ) 00315 DMAX = D( IX ) 00316 CALL SSCAL( N, ANORM / DMAX, D, 1 ) 00317 CALL CSSCAL( N-1, ANORM / DMAX, E, 1 ) 00318 * 00319 ELSE IF( IZERO.GT.0 ) THEN 00320 * 00321 * Reuse the last matrix by copying back the zeroed out 00322 * elements. 00323 * 00324 IF( IZERO.EQ.1 ) THEN 00325 D( 1 ) = Z( 2 ) 00326 IF( N.GT.1 ) 00327 $ E( 1 ) = Z( 3 ) 00328 ELSE IF( IZERO.EQ.N ) THEN 00329 E( N-1 ) = Z( 1 ) 00330 D( N ) = Z( 2 ) 00331 ELSE 00332 E( IZERO-1 ) = Z( 1 ) 00333 D( IZERO ) = Z( 2 ) 00334 E( IZERO ) = Z( 3 ) 00335 END IF 00336 END IF 00337 * 00338 * For types 8-10, set one row and column of the matrix to 00339 * zero. 00340 * 00341 IZERO = 0 00342 IF( IMAT.EQ.8 ) THEN 00343 IZERO = 1 00344 Z( 2 ) = D( 1 ) 00345 D( 1 ) = ZERO 00346 IF( N.GT.1 ) THEN 00347 Z( 3 ) = E( 1 ) 00348 E( 1 ) = ZERO 00349 END IF 00350 ELSE IF( IMAT.EQ.9 ) THEN 00351 IZERO = N 00352 IF( N.GT.1 ) THEN 00353 Z( 1 ) = E( N-1 ) 00354 E( N-1 ) = ZERO 00355 END IF 00356 Z( 2 ) = D( N ) 00357 D( N ) = ZERO 00358 ELSE IF( IMAT.EQ.10 ) THEN 00359 IZERO = ( N+1 ) / 2 00360 IF( IZERO.GT.1 ) THEN 00361 Z( 1 ) = E( IZERO-1 ) 00362 Z( 3 ) = E( IZERO ) 00363 E( IZERO-1 ) = ZERO 00364 E( IZERO ) = ZERO 00365 END IF 00366 Z( 2 ) = D( IZERO ) 00367 D( IZERO ) = ZERO 00368 END IF 00369 END IF 00370 * 00371 CALL SCOPY( N, D, 1, D( N+1 ), 1 ) 00372 IF( N.GT.1 ) 00373 $ CALL CCOPY( N-1, E, 1, E( N+1 ), 1 ) 00374 * 00375 *+ TEST 1 00376 * Factor A as L*D*L' and compute the ratio 00377 * norm(L*D*L' - A) / (n * norm(A) * EPS ) 00378 * 00379 CALL CPTTRF( N, D( N+1 ), E( N+1 ), INFO ) 00380 * 00381 * Check error code from CPTTRF. 00382 * 00383 IF( INFO.NE.IZERO ) THEN 00384 CALL ALAERH( PATH, 'CPTTRF', INFO, IZERO, ' ', N, N, -1, 00385 $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) 00386 GO TO 110 00387 END IF 00388 * 00389 IF( INFO.GT.0 ) THEN 00390 RCONDC = ZERO 00391 GO TO 100 00392 END IF 00393 * 00394 CALL CPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK, 00395 $ RESULT( 1 ) ) 00396 * 00397 * Print the test ratio if greater than or equal to THRESH. 00398 * 00399 IF( RESULT( 1 ).GE.THRESH ) THEN 00400 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00401 $ CALL ALAHD( NOUT, PATH ) 00402 WRITE( NOUT, FMT = 9999 )N, IMAT, 1, RESULT( 1 ) 00403 NFAIL = NFAIL + 1 00404 END IF 00405 NRUN = NRUN + 1 00406 * 00407 * Compute RCONDC = 1 / (norm(A) * norm(inv(A)) 00408 * 00409 * Compute norm(A). 00410 * 00411 ANORM = CLANHT( '1', N, D, E ) 00412 * 00413 * Use CPTTRS to solve for one column at a time of inv(A), 00414 * computing the maximum column sum as we go. 00415 * 00416 AINVNM = ZERO 00417 DO 50 I = 1, N 00418 DO 40 J = 1, N 00419 X( J ) = ZERO 00420 40 CONTINUE 00421 X( I ) = ONE 00422 CALL CPTTRS( 'Lower', N, 1, D( N+1 ), E( N+1 ), X, LDA, 00423 $ INFO ) 00424 AINVNM = MAX( AINVNM, SCASUM( N, X, 1 ) ) 00425 50 CONTINUE 00426 RCONDC = ONE / MAX( ONE, ANORM*AINVNM ) 00427 * 00428 DO 90 IRHS = 1, NNS 00429 NRHS = NSVAL( IRHS ) 00430 * 00431 * Generate NRHS random solution vectors. 00432 * 00433 IX = 1 00434 DO 60 J = 1, NRHS 00435 CALL CLARNV( 2, ISEED, N, XACT( IX ) ) 00436 IX = IX + LDA 00437 60 CONTINUE 00438 * 00439 DO 80 IUPLO = 1, 2 00440 * 00441 * Do first for UPLO = 'U', then for UPLO = 'L'. 00442 * 00443 UPLO = UPLOS( IUPLO ) 00444 * 00445 * Set the right hand side. 00446 * 00447 CALL CLAPTM( UPLO, N, NRHS, ONE, D, E, XACT, LDA, 00448 $ ZERO, B, LDA ) 00449 * 00450 *+ TEST 2 00451 * Solve A*x = b and compute the residual. 00452 * 00453 CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 00454 CALL CPTTRS( UPLO, N, NRHS, D( N+1 ), E( N+1 ), X, 00455 $ LDA, INFO ) 00456 * 00457 * Check error code from CPTTRS. 00458 * 00459 IF( INFO.NE.0 ) 00460 $ CALL ALAERH( PATH, 'CPTTRS', INFO, 0, UPLO, N, N, 00461 $ -1, -1, NRHS, IMAT, NFAIL, NERRS, 00462 $ NOUT ) 00463 * 00464 CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) 00465 CALL CPTT02( UPLO, N, NRHS, D, E, X, LDA, WORK, LDA, 00466 $ RESULT( 2 ) ) 00467 * 00468 *+ TEST 3 00469 * Check solution from generated exact solution. 00470 * 00471 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00472 $ RESULT( 3 ) ) 00473 * 00474 *+ TESTS 4, 5, and 6 00475 * Use iterative refinement to improve the solution. 00476 * 00477 SRNAMT = 'CPTRFS' 00478 CALL CPTRFS( UPLO, N, NRHS, D, E, D( N+1 ), E( N+1 ), 00479 $ B, LDA, X, LDA, RWORK, RWORK( NRHS+1 ), 00480 $ WORK, RWORK( 2*NRHS+1 ), INFO ) 00481 * 00482 * Check error code from CPTRFS. 00483 * 00484 IF( INFO.NE.0 ) 00485 $ CALL ALAERH( PATH, 'CPTRFS', INFO, 0, UPLO, N, N, 00486 $ -1, -1, NRHS, IMAT, NFAIL, NERRS, 00487 $ NOUT ) 00488 * 00489 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 00490 $ RESULT( 4 ) ) 00491 CALL CPTT05( N, NRHS, D, E, B, LDA, X, LDA, XACT, LDA, 00492 $ RWORK, RWORK( NRHS+1 ), RESULT( 5 ) ) 00493 * 00494 * Print information about the tests that did not pass the 00495 * threshold. 00496 * 00497 DO 70 K = 2, 6 00498 IF( RESULT( K ).GE.THRESH ) THEN 00499 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00500 $ CALL ALAHD( NOUT, PATH ) 00501 WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, IMAT, 00502 $ K, RESULT( K ) 00503 NFAIL = NFAIL + 1 00504 END IF 00505 70 CONTINUE 00506 NRUN = NRUN + 5 00507 * 00508 80 CONTINUE 00509 90 CONTINUE 00510 * 00511 *+ TEST 7 00512 * Estimate the reciprocal of the condition number of the 00513 * matrix. 00514 * 00515 100 CONTINUE 00516 SRNAMT = 'CPTCON' 00517 CALL CPTCON( N, D( N+1 ), E( N+1 ), ANORM, RCOND, RWORK, 00518 $ INFO ) 00519 * 00520 * Check error code from CPTCON. 00521 * 00522 IF( INFO.NE.0 ) 00523 $ CALL ALAERH( PATH, 'CPTCON', INFO, 0, ' ', N, N, -1, -1, 00524 $ -1, IMAT, NFAIL, NERRS, NOUT ) 00525 * 00526 RESULT( 7 ) = SGET06( RCOND, RCONDC ) 00527 * 00528 * Print the test ratio if greater than or equal to THRESH. 00529 * 00530 IF( RESULT( 7 ).GE.THRESH ) THEN 00531 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 00532 $ CALL ALAHD( NOUT, PATH ) 00533 WRITE( NOUT, FMT = 9999 )N, IMAT, 7, RESULT( 7 ) 00534 NFAIL = NFAIL + 1 00535 END IF 00536 NRUN = NRUN + 1 00537 110 CONTINUE 00538 120 CONTINUE 00539 * 00540 * Print a summary of the results. 00541 * 00542 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) 00543 * 00544 9999 FORMAT( ' N =', I5, ', type ', I2, ', test ', I2, ', ratio = ', 00545 $ G12.5 ) 00546 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS =', I3, 00547 $ ', type ', I2, ', test ', I2, ', ratio = ', G12.5 ) 00548 RETURN 00549 * 00550 * End of CCHKPT 00551 * 00552 END