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