1*> \brief \b DDRVPT 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE DDRVPT( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, D, 12* E, B, X, XACT, WORK, RWORK, NOUT ) 13* 14* .. Scalar Arguments .. 15* LOGICAL TSTERR 16* INTEGER NN, NOUT, NRHS 17* DOUBLE PRECISION THRESH 18* .. 19* .. Array Arguments .. 20* LOGICAL DOTYPE( * ) 21* INTEGER NVAL( * ) 22* DOUBLE PRECISION A( * ), B( * ), D( * ), E( * ), RWORK( * ), 23* $ WORK( * ), X( * ), XACT( * ) 24* .. 25* 26* 27*> \par Purpose: 28* ============= 29*> 30*> \verbatim 31*> 32*> DDRVPT tests DPTSV and -SVX. 33*> \endverbatim 34* 35* Arguments: 36* ========== 37* 38*> \param[in] DOTYPE 39*> \verbatim 40*> DOTYPE is LOGICAL array, dimension (NTYPES) 41*> The matrix types to be used for testing. Matrices of type j 42*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = 43*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. 44*> \endverbatim 45*> 46*> \param[in] NN 47*> \verbatim 48*> NN is INTEGER 49*> The number of values of N contained in the vector NVAL. 50*> \endverbatim 51*> 52*> \param[in] NVAL 53*> \verbatim 54*> NVAL is INTEGER array, dimension (NN) 55*> The values of the matrix dimension N. 56*> \endverbatim 57*> 58*> \param[in] NRHS 59*> \verbatim 60*> NRHS is INTEGER 61*> The number of right hand side vectors to be generated for 62*> each linear system. 63*> \endverbatim 64*> 65*> \param[in] THRESH 66*> \verbatim 67*> THRESH is DOUBLE PRECISION 68*> The threshold value for the test ratios. A result is 69*> included in the output file if RESULT >= THRESH. To have 70*> every test ratio printed, use THRESH = 0. 71*> \endverbatim 72*> 73*> \param[in] TSTERR 74*> \verbatim 75*> TSTERR is LOGICAL 76*> Flag that indicates whether error exits are to be tested. 77*> \endverbatim 78*> 79*> \param[out] A 80*> \verbatim 81*> A is DOUBLE PRECISION array, dimension (NMAX*2) 82*> \endverbatim 83*> 84*> \param[out] D 85*> \verbatim 86*> D is DOUBLE PRECISION array, dimension (NMAX*2) 87*> \endverbatim 88*> 89*> \param[out] E 90*> \verbatim 91*> E is DOUBLE PRECISION array, dimension (NMAX*2) 92*> \endverbatim 93*> 94*> \param[out] B 95*> \verbatim 96*> B is DOUBLE PRECISION array, dimension (NMAX*NRHS) 97*> \endverbatim 98*> 99*> \param[out] X 100*> \verbatim 101*> X is DOUBLE PRECISION array, dimension (NMAX*NRHS) 102*> \endverbatim 103*> 104*> \param[out] XACT 105*> \verbatim 106*> XACT is DOUBLE PRECISION array, dimension (NMAX*NRHS) 107*> \endverbatim 108*> 109*> \param[out] WORK 110*> \verbatim 111*> WORK is DOUBLE PRECISION array, dimension 112*> (NMAX*max(3,NRHS)) 113*> \endverbatim 114*> 115*> \param[out] RWORK 116*> \verbatim 117*> RWORK is DOUBLE PRECISION array, dimension 118*> (max(NMAX,2*NRHS)) 119*> \endverbatim 120*> 121*> \param[in] NOUT 122*> \verbatim 123*> NOUT is INTEGER 124*> The unit number for output. 125*> \endverbatim 126* 127* Authors: 128* ======== 129* 130*> \author Univ. of Tennessee 131*> \author Univ. of California Berkeley 132*> \author Univ. of Colorado Denver 133*> \author NAG Ltd. 134* 135*> \ingroup double_lin 136* 137* ===================================================================== 138 SUBROUTINE DDRVPT( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, D, 139 $ E, B, X, XACT, WORK, RWORK, NOUT ) 140* 141* -- LAPACK test routine -- 142* -- LAPACK is a software package provided by Univ. of Tennessee, -- 143* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 144* 145* .. Scalar Arguments .. 146 LOGICAL TSTERR 147 INTEGER NN, NOUT, NRHS 148 DOUBLE PRECISION THRESH 149* .. 150* .. Array Arguments .. 151 LOGICAL DOTYPE( * ) 152 INTEGER NVAL( * ) 153 DOUBLE PRECISION A( * ), B( * ), D( * ), E( * ), RWORK( * ), 154 $ WORK( * ), X( * ), XACT( * ) 155* .. 156* 157* ===================================================================== 158* 159* .. Parameters .. 160 DOUBLE PRECISION ONE, ZERO 161 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 162 INTEGER NTYPES 163 PARAMETER ( NTYPES = 12 ) 164 INTEGER NTESTS 165 PARAMETER ( NTESTS = 6 ) 166* .. 167* .. Local Scalars .. 168 LOGICAL ZEROT 169 CHARACTER DIST, FACT, TYPE 170 CHARACTER*3 PATH 171 INTEGER I, IA, IFACT, IMAT, IN, INFO, IX, IZERO, J, K, 172 $ K1, KL, KU, LDA, MODE, N, NERRS, NFAIL, NIMAT, 173 $ NRUN, NT 174 DOUBLE PRECISION AINVNM, ANORM, COND, DMAX, RCOND, RCONDC 175* .. 176* .. Local Arrays .. 177 INTEGER ISEED( 4 ), ISEEDY( 4 ) 178 DOUBLE PRECISION RESULT( NTESTS ), Z( 3 ) 179* .. 180* .. External Functions .. 181 INTEGER IDAMAX 182 DOUBLE PRECISION DASUM, DGET06, DLANST 183 EXTERNAL IDAMAX, DASUM, DGET06, DLANST 184* .. 185* .. External Subroutines .. 186 EXTERNAL ALADHD, ALAERH, ALASVM, DCOPY, DERRVX, DGET04, 187 $ DLACPY, DLAPTM, DLARNV, DLASET, DLATB4, DLATMS, 188 $ DPTSV, DPTSVX, DPTT01, DPTT02, DPTT05, DPTTRF, 189 $ DPTTRS, DSCAL 190* .. 191* .. Intrinsic Functions .. 192 INTRINSIC ABS, MAX 193* .. 194* .. Scalars in Common .. 195 LOGICAL LERR, OK 196 CHARACTER*32 SRNAMT 197 INTEGER INFOT, NUNIT 198* .. 199* .. Common blocks .. 200 COMMON / INFOC / INFOT, NUNIT, OK, LERR 201 COMMON / SRNAMC / SRNAMT 202* .. 203* .. Data statements .. 204 DATA ISEEDY / 0, 0, 0, 1 / 205* .. 206* .. Executable Statements .. 207* 208 PATH( 1: 1 ) = 'Double precision' 209 PATH( 2: 3 ) = 'PT' 210 NRUN = 0 211 NFAIL = 0 212 NERRS = 0 213 DO 10 I = 1, 4 214 ISEED( I ) = ISEEDY( I ) 215 10 CONTINUE 216* 217* Test the error exits 218* 219 IF( TSTERR ) 220 $ CALL DERRVX( PATH, NOUT ) 221 INFOT = 0 222* 223 DO 120 IN = 1, NN 224* 225* Do for each value of N in NVAL. 226* 227 N = NVAL( IN ) 228 LDA = MAX( 1, N ) 229 NIMAT = NTYPES 230 IF( N.LE.0 ) 231 $ NIMAT = 1 232* 233 DO 110 IMAT = 1, NIMAT 234* 235* Do the tests only if DOTYPE( IMAT ) is true. 236* 237 IF( N.GT.0 .AND. .NOT.DOTYPE( IMAT ) ) 238 $ GO TO 110 239* 240* Set up parameters with DLATB4. 241* 242 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 243 $ COND, DIST ) 244* 245 ZEROT = IMAT.GE.8 .AND. IMAT.LE.10 246 IF( IMAT.LE.6 ) THEN 247* 248* Type 1-6: generate a symmetric tridiagonal matrix of 249* known condition number in lower triangular band storage. 250* 251 SRNAMT = 'DLATMS' 252 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, COND, 253 $ ANORM, KL, KU, 'B', A, 2, WORK, INFO ) 254* 255* Check the error code from DLATMS. 256* 257 IF( INFO.NE.0 ) THEN 258 CALL ALAERH( PATH, 'DLATMS', INFO, 0, ' ', N, N, KL, 259 $ KU, -1, IMAT, NFAIL, NERRS, NOUT ) 260 GO TO 110 261 END IF 262 IZERO = 0 263* 264* Copy the matrix to D and E. 265* 266 IA = 1 267 DO 20 I = 1, N - 1 268 D( I ) = A( IA ) 269 E( I ) = A( IA+1 ) 270 IA = IA + 2 271 20 CONTINUE 272 IF( N.GT.0 ) 273 $ D( N ) = A( IA ) 274 ELSE 275* 276* Type 7-12: generate a diagonally dominant matrix with 277* unknown condition number in the vectors D and E. 278* 279 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 7 ) ) THEN 280* 281* Let D and E have values from [-1,1]. 282* 283 CALL DLARNV( 2, ISEED, N, D ) 284 CALL DLARNV( 2, ISEED, N-1, E ) 285* 286* Make the tridiagonal matrix diagonally dominant. 287* 288 IF( N.EQ.1 ) THEN 289 D( 1 ) = ABS( D( 1 ) ) 290 ELSE 291 D( 1 ) = ABS( D( 1 ) ) + ABS( E( 1 ) ) 292 D( N ) = ABS( D( N ) ) + ABS( E( N-1 ) ) 293 DO 30 I = 2, N - 1 294 D( I ) = ABS( D( I ) ) + ABS( E( I ) ) + 295 $ ABS( E( I-1 ) ) 296 30 CONTINUE 297 END IF 298* 299* Scale D and E so the maximum element is ANORM. 300* 301 IX = IDAMAX( N, D, 1 ) 302 DMAX = D( IX ) 303 CALL DSCAL( N, ANORM / DMAX, D, 1 ) 304 IF( N.GT.1 ) 305 $ CALL DSCAL( N-1, ANORM / DMAX, E, 1 ) 306* 307 ELSE IF( IZERO.GT.0 ) THEN 308* 309* Reuse the last matrix by copying back the zeroed out 310* elements. 311* 312 IF( IZERO.EQ.1 ) THEN 313 D( 1 ) = Z( 2 ) 314 IF( N.GT.1 ) 315 $ E( 1 ) = Z( 3 ) 316 ELSE IF( IZERO.EQ.N ) THEN 317 E( N-1 ) = Z( 1 ) 318 D( N ) = Z( 2 ) 319 ELSE 320 E( IZERO-1 ) = Z( 1 ) 321 D( IZERO ) = Z( 2 ) 322 E( IZERO ) = Z( 3 ) 323 END IF 324 END IF 325* 326* For types 8-10, set one row and column of the matrix to 327* zero. 328* 329 IZERO = 0 330 IF( IMAT.EQ.8 ) THEN 331 IZERO = 1 332 Z( 2 ) = D( 1 ) 333 D( 1 ) = ZERO 334 IF( N.GT.1 ) THEN 335 Z( 3 ) = E( 1 ) 336 E( 1 ) = ZERO 337 END IF 338 ELSE IF( IMAT.EQ.9 ) THEN 339 IZERO = N 340 IF( N.GT.1 ) THEN 341 Z( 1 ) = E( N-1 ) 342 E( N-1 ) = ZERO 343 END IF 344 Z( 2 ) = D( N ) 345 D( N ) = ZERO 346 ELSE IF( IMAT.EQ.10 ) THEN 347 IZERO = ( N+1 ) / 2 348 IF( IZERO.GT.1 ) THEN 349 Z( 1 ) = E( IZERO-1 ) 350 Z( 3 ) = E( IZERO ) 351 E( IZERO-1 ) = ZERO 352 E( IZERO ) = ZERO 353 END IF 354 Z( 2 ) = D( IZERO ) 355 D( IZERO ) = ZERO 356 END IF 357 END IF 358* 359* Generate NRHS random solution vectors. 360* 361 IX = 1 362 DO 40 J = 1, NRHS 363 CALL DLARNV( 2, ISEED, N, XACT( IX ) ) 364 IX = IX + LDA 365 40 CONTINUE 366* 367* Set the right hand side. 368* 369 CALL DLAPTM( N, NRHS, ONE, D, E, XACT, LDA, ZERO, B, LDA ) 370* 371 DO 100 IFACT = 1, 2 372 IF( IFACT.EQ.1 ) THEN 373 FACT = 'F' 374 ELSE 375 FACT = 'N' 376 END IF 377* 378* Compute the condition number for comparison with 379* the value returned by DPTSVX. 380* 381 IF( ZEROT ) THEN 382 IF( IFACT.EQ.1 ) 383 $ GO TO 100 384 RCONDC = ZERO 385* 386 ELSE IF( IFACT.EQ.1 ) THEN 387* 388* Compute the 1-norm of A. 389* 390 ANORM = DLANST( '1', N, D, E ) 391* 392 CALL DCOPY( N, D, 1, D( N+1 ), 1 ) 393 IF( N.GT.1 ) 394 $ CALL DCOPY( N-1, E, 1, E( N+1 ), 1 ) 395* 396* Factor the matrix A. 397* 398 CALL DPTTRF( N, D( N+1 ), E( N+1 ), INFO ) 399* 400* Use DPTTRS to solve for one column at a time of 401* inv(A), computing the maximum column sum as we go. 402* 403 AINVNM = ZERO 404 DO 60 I = 1, N 405 DO 50 J = 1, N 406 X( J ) = ZERO 407 50 CONTINUE 408 X( I ) = ONE 409 CALL DPTTRS( N, 1, D( N+1 ), E( N+1 ), X, LDA, 410 $ INFO ) 411 AINVNM = MAX( AINVNM, DASUM( N, X, 1 ) ) 412 60 CONTINUE 413* 414* Compute the 1-norm condition number of A. 415* 416 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN 417 RCONDC = ONE 418 ELSE 419 RCONDC = ( ONE / ANORM ) / AINVNM 420 END IF 421 END IF 422* 423 IF( IFACT.EQ.2 ) THEN 424* 425* --- Test DPTSV -- 426* 427 CALL DCOPY( N, D, 1, D( N+1 ), 1 ) 428 IF( N.GT.1 ) 429 $ CALL DCOPY( N-1, E, 1, E( N+1 ), 1 ) 430 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) 431* 432* Factor A as L*D*L' and solve the system A*X = B. 433* 434 SRNAMT = 'DPTSV ' 435 CALL DPTSV( N, NRHS, D( N+1 ), E( N+1 ), X, LDA, 436 $ INFO ) 437* 438* Check error code from DPTSV . 439* 440 IF( INFO.NE.IZERO ) 441 $ CALL ALAERH( PATH, 'DPTSV ', INFO, IZERO, ' ', N, 442 $ N, 1, 1, NRHS, IMAT, NFAIL, NERRS, 443 $ NOUT ) 444 NT = 0 445 IF( IZERO.EQ.0 ) THEN 446* 447* Check the factorization by computing the ratio 448* norm(L*D*L' - A) / (n * norm(A) * EPS ) 449* 450 CALL DPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK, 451 $ RESULT( 1 ) ) 452* 453* Compute the residual in the solution. 454* 455 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) 456 CALL DPTT02( N, NRHS, D, E, X, LDA, WORK, LDA, 457 $ RESULT( 2 ) ) 458* 459* Check solution from generated exact solution. 460* 461 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 462 $ RESULT( 3 ) ) 463 NT = 3 464 END IF 465* 466* Print information about the tests that did not pass 467* the threshold. 468* 469 DO 70 K = 1, NT 470 IF( RESULT( K ).GE.THRESH ) THEN 471 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 472 $ CALL ALADHD( NOUT, PATH ) 473 WRITE( NOUT, FMT = 9999 )'DPTSV ', N, IMAT, K, 474 $ RESULT( K ) 475 NFAIL = NFAIL + 1 476 END IF 477 70 CONTINUE 478 NRUN = NRUN + NT 479 END IF 480* 481* --- Test DPTSVX --- 482* 483 IF( IFACT.GT.1 ) THEN 484* 485* Initialize D( N+1:2*N ) and E( N+1:2*N ) to zero. 486* 487 DO 80 I = 1, N - 1 488 D( N+I ) = ZERO 489 E( N+I ) = ZERO 490 80 CONTINUE 491 IF( N.GT.0 ) 492 $ D( N+N ) = ZERO 493 END IF 494* 495 CALL DLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA ) 496* 497* Solve the system and compute the condition number and 498* error bounds using DPTSVX. 499* 500 SRNAMT = 'DPTSVX' 501 CALL DPTSVX( FACT, N, NRHS, D, E, D( N+1 ), E( N+1 ), B, 502 $ LDA, X, LDA, RCOND, RWORK, RWORK( NRHS+1 ), 503 $ WORK, INFO ) 504* 505* Check the error code from DPTSVX. 506* 507 IF( INFO.NE.IZERO ) 508 $ CALL ALAERH( PATH, 'DPTSVX', INFO, IZERO, FACT, N, N, 509 $ 1, 1, NRHS, IMAT, NFAIL, NERRS, NOUT ) 510 IF( IZERO.EQ.0 ) THEN 511 IF( IFACT.EQ.2 ) THEN 512* 513* Check the factorization by computing the ratio 514* norm(L*D*L' - A) / (n * norm(A) * EPS ) 515* 516 K1 = 1 517 CALL DPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK, 518 $ RESULT( 1 ) ) 519 ELSE 520 K1 = 2 521 END IF 522* 523* Compute the residual in the solution. 524* 525 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) 526 CALL DPTT02( N, NRHS, D, E, X, LDA, WORK, LDA, 527 $ RESULT( 2 ) ) 528* 529* Check solution from generated exact solution. 530* 531 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 532 $ RESULT( 3 ) ) 533* 534* Check error bounds from iterative refinement. 535* 536 CALL DPTT05( N, NRHS, D, E, B, LDA, X, LDA, XACT, LDA, 537 $ RWORK, RWORK( NRHS+1 ), RESULT( 4 ) ) 538 ELSE 539 K1 = 6 540 END IF 541* 542* Check the reciprocal of the condition number. 543* 544 RESULT( 6 ) = DGET06( RCOND, RCONDC ) 545* 546* Print information about the tests that did not pass 547* the threshold. 548* 549 DO 90 K = K1, 6 550 IF( RESULT( K ).GE.THRESH ) THEN 551 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 552 $ CALL ALADHD( NOUT, PATH ) 553 WRITE( NOUT, FMT = 9998 )'DPTSVX', FACT, N, IMAT, 554 $ K, RESULT( K ) 555 NFAIL = NFAIL + 1 556 END IF 557 90 CONTINUE 558 NRUN = NRUN + 7 - K1 559 100 CONTINUE 560 110 CONTINUE 561 120 CONTINUE 562* 563* Print a summary of the results. 564* 565 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) 566* 567 9999 FORMAT( 1X, A, ', N =', I5, ', type ', I2, ', test ', I2, 568 $ ', ratio = ', G12.5 ) 569 9998 FORMAT( 1X, A, ', FACT=''', A1, ''', N =', I5, ', type ', I2, 570 $ ', test ', I2, ', ratio = ', G12.5 ) 571 RETURN 572* 573* End of DDRVPT 574* 575 END 576