1*> \brief \b DDRVRFP 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 DDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL, 12* + THRESH, A, ASAV, AFAC, AINV, B, 13* + BSAV, XACT, X, ARF, ARFINV, 14* + D_WORK_DLATMS, D_WORK_DPOT01, D_TEMP_DPOT02, 15* + D_TEMP_DPOT03, D_WORK_DLANSY, 16* + D_WORK_DPOT02, D_WORK_DPOT03 ) 17* 18* .. Scalar Arguments .. 19* INTEGER NN, NNS, NNT, NOUT 20* DOUBLE PRECISION THRESH 21* .. 22* .. Array Arguments .. 23* INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT ) 24* DOUBLE PRECISION A( * ) 25* DOUBLE PRECISION AINV( * ) 26* DOUBLE PRECISION ASAV( * ) 27* DOUBLE PRECISION B( * ) 28* DOUBLE PRECISION BSAV( * ) 29* DOUBLE PRECISION AFAC( * ) 30* DOUBLE PRECISION ARF( * ) 31* DOUBLE PRECISION ARFINV( * ) 32* DOUBLE PRECISION XACT( * ) 33* DOUBLE PRECISION X( * ) 34* DOUBLE PRECISION D_WORK_DLATMS( * ) 35* DOUBLE PRECISION D_WORK_DPOT01( * ) 36* DOUBLE PRECISION D_TEMP_DPOT02( * ) 37* DOUBLE PRECISION D_TEMP_DPOT03( * ) 38* DOUBLE PRECISION D_WORK_DLANSY( * ) 39* DOUBLE PRECISION D_WORK_DPOT02( * ) 40* DOUBLE PRECISION D_WORK_DPOT03( * ) 41* .. 42* 43* 44*> \par Purpose: 45* ============= 46*> 47*> \verbatim 48*> 49*> DDRVRFP tests the LAPACK RFP routines: 50*> DPFTRF, DPFTRS, and DPFTRI. 51*> 52*> This testing routine follow the same tests as DDRVPO (test for the full 53*> format Symmetric Positive Definite solver). 54*> 55*> The tests are performed in Full Format, convertion back and forth from 56*> full format to RFP format are performed using the routines DTRTTF and 57*> DTFTTR. 58*> 59*> First, a specific matrix A of size N is created. There is nine types of 60*> different matrixes possible. 61*> 1. Diagonal 6. Random, CNDNUM = sqrt(0.1/EPS) 62*> 2. Random, CNDNUM = 2 7. Random, CNDNUM = 0.1/EPS 63*> *3. First row and column zero 8. Scaled near underflow 64*> *4. Last row and column zero 9. Scaled near overflow 65*> *5. Middle row and column zero 66*> (* - tests error exits from DPFTRF, no test ratios are computed) 67*> A solution XACT of size N-by-NRHS is created and the associated right 68*> hand side B as well. Then DPFTRF is called to compute L (or U), the 69*> Cholesky factor of A. Then L (or U) is used to solve the linear system 70*> of equations AX = B. This gives X. Then L (or U) is used to compute the 71*> inverse of A, AINV. The following four tests are then performed: 72*> (1) norm( L*L' - A ) / ( N * norm(A) * EPS ) or 73*> norm( U'*U - A ) / ( N * norm(A) * EPS ), 74*> (2) norm(B - A*X) / ( norm(A) * norm(X) * EPS ), 75*> (3) norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ), 76*> (4) ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ), 77*> where EPS is the machine precision, RCOND the condition number of A, and 78*> norm( . ) the 1-norm for (1,2,3) and the inf-norm for (4). 79*> Errors occur when INFO parameter is not as expected. Failures occur when 80*> a test ratios is greater than THRES. 81*> \endverbatim 82* 83* Arguments: 84* ========== 85* 86*> \param[in] NOUT 87*> \verbatim 88*> NOUT is INTEGER 89*> The unit number for output. 90*> \endverbatim 91*> 92*> \param[in] NN 93*> \verbatim 94*> NN is INTEGER 95*> The number of values of N contained in the vector NVAL. 96*> \endverbatim 97*> 98*> \param[in] NVAL 99*> \verbatim 100*> NVAL is INTEGER array, dimension (NN) 101*> The values of the matrix dimension N. 102*> \endverbatim 103*> 104*> \param[in] NNS 105*> \verbatim 106*> NNS is INTEGER 107*> The number of values of NRHS contained in the vector NSVAL. 108*> \endverbatim 109*> 110*> \param[in] NSVAL 111*> \verbatim 112*> NSVAL is INTEGER array, dimension (NNS) 113*> The values of the number of right-hand sides NRHS. 114*> \endverbatim 115*> 116*> \param[in] NNT 117*> \verbatim 118*> NNT is INTEGER 119*> The number of values of MATRIX TYPE contained in the vector NTVAL. 120*> \endverbatim 121*> 122*> \param[in] NTVAL 123*> \verbatim 124*> NTVAL is INTEGER array, dimension (NNT) 125*> The values of matrix type (between 0 and 9 for PO/PP/PF matrices). 126*> \endverbatim 127*> 128*> \param[in] THRESH 129*> \verbatim 130*> THRESH is DOUBLE PRECISION 131*> The threshold value for the test ratios. A result is 132*> included in the output file if RESULT >= THRESH. To have 133*> every test ratio printed, use THRESH = 0. 134*> \endverbatim 135*> 136*> \param[out] A 137*> \verbatim 138*> A is DOUBLE PRECISION array, dimension (NMAX*NMAX) 139*> \endverbatim 140*> 141*> \param[out] ASAV 142*> \verbatim 143*> ASAV is DOUBLE PRECISION array, dimension (NMAX*NMAX) 144*> \endverbatim 145*> 146*> \param[out] AFAC 147*> \verbatim 148*> AFAC is DOUBLE PRECISION array, dimension (NMAX*NMAX) 149*> \endverbatim 150*> 151*> \param[out] AINV 152*> \verbatim 153*> AINV is DOUBLE PRECISION array, dimension (NMAX*NMAX) 154*> \endverbatim 155*> 156*> \param[out] B 157*> \verbatim 158*> B is DOUBLE PRECISION array, dimension (NMAX*MAXRHS) 159*> \endverbatim 160*> 161*> \param[out] BSAV 162*> \verbatim 163*> BSAV is DOUBLE PRECISION array, dimension (NMAX*MAXRHS) 164*> \endverbatim 165*> 166*> \param[out] XACT 167*> \verbatim 168*> XACT is DOUBLE PRECISION array, dimension (NMAX*MAXRHS) 169*> \endverbatim 170*> 171*> \param[out] X 172*> \verbatim 173*> X is DOUBLE PRECISION array, dimension (NMAX*MAXRHS) 174*> \endverbatim 175*> 176*> \param[out] ARF 177*> \verbatim 178*> ARF is DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2) 179*> \endverbatim 180*> 181*> \param[out] ARFINV 182*> \verbatim 183*> ARFINV is DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2) 184*> \endverbatim 185*> 186*> \param[out] D_WORK_DLATMS 187*> \verbatim 188*> D_WORK_DLATMS is DOUBLE PRECISION array, dimension ( 3*NMAX ) 189*> \endverbatim 190*> 191*> \param[out] D_WORK_DPOT01 192*> \verbatim 193*> D_WORK_DPOT01 is DOUBLE PRECISION array, dimension ( NMAX ) 194*> \endverbatim 195*> 196*> \param[out] D_TEMP_DPOT02 197*> \verbatim 198*> D_TEMP_DPOT02 is DOUBLE PRECISION array, dimension ( NMAX*MAXRHS ) 199*> \endverbatim 200*> 201*> \param[out] D_TEMP_DPOT03 202*> \verbatim 203*> D_TEMP_DPOT03 is DOUBLE PRECISION array, dimension ( NMAX*NMAX ) 204*> \endverbatim 205*> 206*> \param[out] D_WORK_DLATMS 207*> \verbatim 208*> D_WORK_DLATMS is DOUBLE PRECISION array, dimension ( NMAX ) 209*> \endverbatim 210*> 211*> \param[out] D_WORK_DLANSY 212*> \verbatim 213*> D_WORK_DLANSY is DOUBLE PRECISION array, dimension ( NMAX ) 214*> \endverbatim 215*> 216*> \param[out] D_WORK_DPOT02 217*> \verbatim 218*> D_WORK_DPOT02 is DOUBLE PRECISION array, dimension ( NMAX ) 219*> \endverbatim 220*> 221*> \param[out] D_WORK_DPOT03 222*> \verbatim 223*> D_WORK_DPOT03 is DOUBLE PRECISION array, dimension ( NMAX ) 224*> \endverbatim 225* 226* Authors: 227* ======== 228* 229*> \author Univ. of Tennessee 230*> \author Univ. of California Berkeley 231*> \author Univ. of Colorado Denver 232*> \author NAG Ltd. 233* 234*> \date November 2011 235* 236*> \ingroup double_lin 237* 238* ===================================================================== 239 SUBROUTINE DDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL, 240 + THRESH, A, ASAV, AFAC, AINV, B, 241 + BSAV, XACT, X, ARF, ARFINV, 242 + D_WORK_DLATMS, D_WORK_DPOT01, D_TEMP_DPOT02, 243 + D_TEMP_DPOT03, D_WORK_DLANSY, 244 + D_WORK_DPOT02, D_WORK_DPOT03 ) 245* 246* -- LAPACK test routine (version 3.4.0) -- 247* -- LAPACK is a software package provided by Univ. of Tennessee, -- 248* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 249* November 2011 250* 251* .. Scalar Arguments .. 252 INTEGER NN, NNS, NNT, NOUT 253 DOUBLE PRECISION THRESH 254* .. 255* .. Array Arguments .. 256 INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT ) 257 DOUBLE PRECISION A( * ) 258 DOUBLE PRECISION AINV( * ) 259 DOUBLE PRECISION ASAV( * ) 260 DOUBLE PRECISION B( * ) 261 DOUBLE PRECISION BSAV( * ) 262 DOUBLE PRECISION AFAC( * ) 263 DOUBLE PRECISION ARF( * ) 264 DOUBLE PRECISION ARFINV( * ) 265 DOUBLE PRECISION XACT( * ) 266 DOUBLE PRECISION X( * ) 267 DOUBLE PRECISION D_WORK_DLATMS( * ) 268 DOUBLE PRECISION D_WORK_DPOT01( * ) 269 DOUBLE PRECISION D_TEMP_DPOT02( * ) 270 DOUBLE PRECISION D_TEMP_DPOT03( * ) 271 DOUBLE PRECISION D_WORK_DLANSY( * ) 272 DOUBLE PRECISION D_WORK_DPOT02( * ) 273 DOUBLE PRECISION D_WORK_DPOT03( * ) 274* .. 275* 276* ===================================================================== 277* 278* .. Parameters .. 279 DOUBLE PRECISION ONE, ZERO 280 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 281 INTEGER NTESTS 282 PARAMETER ( NTESTS = 4 ) 283* .. 284* .. Local Scalars .. 285 LOGICAL ZEROT 286 INTEGER I, INFO, IUPLO, LDA, LDB, IMAT, NERRS, NFAIL, 287 + NRHS, NRUN, IZERO, IOFF, K, NT, N, IFORM, IIN, 288 + IIT, IIS 289 CHARACTER DIST, CTYPE, UPLO, CFORM 290 INTEGER KL, KU, MODE 291 DOUBLE PRECISION ANORM, AINVNM, CNDNUM, RCONDC 292* .. 293* .. Local Arrays .. 294 CHARACTER UPLOS( 2 ), FORMS( 2 ) 295 INTEGER ISEED( 4 ), ISEEDY( 4 ) 296 DOUBLE PRECISION RESULT( NTESTS ) 297* .. 298* .. External Functions .. 299 DOUBLE PRECISION DLANSY 300 EXTERNAL DLANSY 301* .. 302* .. External Subroutines .. 303 EXTERNAL ALADHD, ALAERH, ALASVM, DGET04, DTFTTR, DLACPY, 304 + DLARHS, DLATB4, DLATMS, DPFTRI, DPFTRF, DPFTRS, 305 + DPOT01, DPOT02, DPOT03, DPOTRI, DPOTRF, DTRTTF 306* .. 307* .. Scalars in Common .. 308 CHARACTER*32 SRNAMT 309* .. 310* .. Common blocks .. 311 COMMON / SRNAMC / SRNAMT 312* .. 313* .. Data statements .. 314 DATA ISEEDY / 1988, 1989, 1990, 1991 / 315 DATA UPLOS / 'U', 'L' / 316 DATA FORMS / 'N', 'T' / 317* .. 318* .. Executable Statements .. 319* 320* Initialize constants and the random number seed. 321* 322 NRUN = 0 323 NFAIL = 0 324 NERRS = 0 325 DO 10 I = 1, 4 326 ISEED( I ) = ISEEDY( I ) 327 10 CONTINUE 328* 329 DO 130 IIN = 1, NN 330* 331 N = NVAL( IIN ) 332 LDA = MAX( N, 1 ) 333 LDB = MAX( N, 1 ) 334* 335 DO 980 IIS = 1, NNS 336* 337 NRHS = NSVAL( IIS ) 338* 339 DO 120 IIT = 1, NNT 340* 341 IMAT = NTVAL( IIT ) 342* 343* If N.EQ.0, only consider the first type 344* 345 IF( N.EQ.0 .AND. IIT.GT.1 ) GO TO 120 346* 347* Skip types 3, 4, or 5 if the matrix size is too small. 348* 349 IF( IMAT.EQ.4 .AND. N.LE.1 ) GO TO 120 350 IF( IMAT.EQ.5 .AND. N.LE.2 ) GO TO 120 351* 352* Do first for UPLO = 'U', then for UPLO = 'L' 353* 354 DO 110 IUPLO = 1, 2 355 UPLO = UPLOS( IUPLO ) 356* 357* Do first for CFORM = 'N', then for CFORM = 'C' 358* 359 DO 100 IFORM = 1, 2 360 CFORM = FORMS( IFORM ) 361* 362* Set up parameters with DLATB4 and generate a test 363* matrix with DLATMS. 364* 365 CALL DLATB4( 'DPO', IMAT, N, N, CTYPE, KL, KU, 366 + ANORM, MODE, CNDNUM, DIST ) 367* 368 SRNAMT = 'DLATMS' 369 CALL DLATMS( N, N, DIST, ISEED, CTYPE, 370 + D_WORK_DLATMS, 371 + MODE, CNDNUM, ANORM, KL, KU, UPLO, A, 372 + LDA, D_WORK_DLATMS, INFO ) 373* 374* Check error code from DLATMS. 375* 376 IF( INFO.NE.0 ) THEN 377 CALL ALAERH( 'DPF', 'DLATMS', INFO, 0, UPLO, N, 378 + N, -1, -1, -1, IIT, NFAIL, NERRS, 379 + NOUT ) 380 GO TO 100 381 END IF 382* 383* For types 3-5, zero one row and column of the matrix to 384* test that INFO is returned correctly. 385* 386 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5 387 IF( ZEROT ) THEN 388 IF( IIT.EQ.3 ) THEN 389 IZERO = 1 390 ELSE IF( IIT.EQ.4 ) THEN 391 IZERO = N 392 ELSE 393 IZERO = N / 2 + 1 394 END IF 395 IOFF = ( IZERO-1 )*LDA 396* 397* Set row and column IZERO of A to 0. 398* 399 IF( IUPLO.EQ.1 ) THEN 400 DO 20 I = 1, IZERO - 1 401 A( IOFF+I ) = ZERO 402 20 CONTINUE 403 IOFF = IOFF + IZERO 404 DO 30 I = IZERO, N 405 A( IOFF ) = ZERO 406 IOFF = IOFF + LDA 407 30 CONTINUE 408 ELSE 409 IOFF = IZERO 410 DO 40 I = 1, IZERO - 1 411 A( IOFF ) = ZERO 412 IOFF = IOFF + LDA 413 40 CONTINUE 414 IOFF = IOFF - IZERO 415 DO 50 I = IZERO, N 416 A( IOFF+I ) = ZERO 417 50 CONTINUE 418 END IF 419 ELSE 420 IZERO = 0 421 END IF 422* 423* Save a copy of the matrix A in ASAV. 424* 425 CALL DLACPY( UPLO, N, N, A, LDA, ASAV, LDA ) 426* 427* Compute the condition number of A (RCONDC). 428* 429 IF( ZEROT ) THEN 430 RCONDC = ZERO 431 ELSE 432* 433* Compute the 1-norm of A. 434* 435 ANORM = DLANSY( '1', UPLO, N, A, LDA, 436 + D_WORK_DLANSY ) 437* 438* Factor the matrix A. 439* 440 CALL DPOTRF( UPLO, N, A, LDA, INFO ) 441* 442* Form the inverse of A. 443* 444 CALL DPOTRI( UPLO, N, A, LDA, INFO ) 445* 446* Compute the 1-norm condition number of A. 447* 448 AINVNM = DLANSY( '1', UPLO, N, A, LDA, 449 + D_WORK_DLANSY ) 450 RCONDC = ( ONE / ANORM ) / AINVNM 451* 452* Restore the matrix A. 453* 454 CALL DLACPY( UPLO, N, N, ASAV, LDA, A, LDA ) 455* 456 END IF 457* 458* Form an exact solution and set the right hand side. 459* 460 SRNAMT = 'DLARHS' 461 CALL DLARHS( 'DPO', 'N', UPLO, ' ', N, N, KL, KU, 462 + NRHS, A, LDA, XACT, LDA, B, LDA, 463 + ISEED, INFO ) 464 CALL DLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA ) 465* 466* Compute the L*L' or U'*U factorization of the 467* matrix and solve the system. 468* 469 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) 470 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDB ) 471* 472 SRNAMT = 'DTRTTF' 473 CALL DTRTTF( CFORM, UPLO, N, AFAC, LDA, ARF, INFO ) 474 SRNAMT = 'DPFTRF' 475 CALL DPFTRF( CFORM, UPLO, N, ARF, INFO ) 476* 477* Check error code from DPFTRF. 478* 479 IF( INFO.NE.IZERO ) THEN 480* 481* LANGOU: there is a small hick here: IZERO should 482* always be INFO however if INFO is ZERO, ALAERH does not 483* complain. 484* 485 CALL ALAERH( 'DPF', 'DPFSV ', INFO, IZERO, 486 + UPLO, N, N, -1, -1, NRHS, IIT, 487 + NFAIL, NERRS, NOUT ) 488 GO TO 100 489 END IF 490* 491* Skip the tests if INFO is not 0. 492* 493 IF( INFO.NE.0 ) THEN 494 GO TO 100 495 END IF 496* 497 SRNAMT = 'DPFTRS' 498 CALL DPFTRS( CFORM, UPLO, N, NRHS, ARF, X, LDB, 499 + INFO ) 500* 501 SRNAMT = 'DTFTTR' 502 CALL DTFTTR( CFORM, UPLO, N, ARF, AFAC, LDA, INFO ) 503* 504* Reconstruct matrix from factors and compute 505* residual. 506* 507 CALL DLACPY( UPLO, N, N, AFAC, LDA, ASAV, LDA ) 508 CALL DPOT01( UPLO, N, A, LDA, AFAC, LDA, 509 + D_WORK_DPOT01, RESULT( 1 ) ) 510 CALL DLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA ) 511* 512* Form the inverse and compute the residual. 513* 514 IF(MOD(N,2).EQ.0)THEN 515 CALL DLACPY( 'A', N+1, N/2, ARF, N+1, ARFINV, 516 + N+1 ) 517 ELSE 518 CALL DLACPY( 'A', N, (N+1)/2, ARF, N, ARFINV, 519 + N ) 520 END IF 521* 522 SRNAMT = 'DPFTRI' 523 CALL DPFTRI( CFORM, UPLO, N, ARFINV , INFO ) 524* 525 SRNAMT = 'DTFTTR' 526 CALL DTFTTR( CFORM, UPLO, N, ARFINV, AINV, LDA, 527 + INFO ) 528* 529* Check error code from DPFTRI. 530* 531 IF( INFO.NE.0 ) 532 + CALL ALAERH( 'DPO', 'DPFTRI', INFO, 0, UPLO, N, 533 + N, -1, -1, -1, IMAT, NFAIL, NERRS, 534 + NOUT ) 535* 536 CALL DPOT03( UPLO, N, A, LDA, AINV, LDA, 537 + D_TEMP_DPOT03, LDA, D_WORK_DPOT03, 538 + RCONDC, RESULT( 2 ) ) 539* 540* Compute residual of the computed solution. 541* 542 CALL DLACPY( 'Full', N, NRHS, B, LDA, 543 + D_TEMP_DPOT02, LDA ) 544 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, 545 + D_TEMP_DPOT02, LDA, D_WORK_DPOT02, 546 + RESULT( 3 ) ) 547* 548* Check solution from generated exact solution. 549 550 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 551 + RESULT( 4 ) ) 552 NT = 4 553* 554* Print information about the tests that did not 555* pass the threshold. 556* 557 DO 60 K = 1, NT 558 IF( RESULT( K ).GE.THRESH ) THEN 559 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 560 + CALL ALADHD( NOUT, 'DPF' ) 561 WRITE( NOUT, FMT = 9999 )'DPFSV ', UPLO, 562 + N, IIT, K, RESULT( K ) 563 NFAIL = NFAIL + 1 564 END IF 565 60 CONTINUE 566 NRUN = NRUN + NT 567 100 CONTINUE 568 110 CONTINUE 569 120 CONTINUE 570 980 CONTINUE 571 130 CONTINUE 572* 573* Print a summary of the results. 574* 575 CALL ALASVM( 'DPF', NOUT, NFAIL, NRUN, NERRS ) 576* 577 9999 FORMAT( 1X, A6, ', UPLO=''', A1, ''', N =', I5, ', type ', I1, 578 + ', test(', I1, ')=', G12.5 ) 579* 580 RETURN 581* 582* End of DDRVRFP 583* 584 END 585