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