1*> \brief \b CDRVRFP 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 CDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL, 12* + THRESH, A, ASAV, AFAC, AINV, B, 13* + BSAV, XACT, X, ARF, ARFINV, 14* + C_WORK_CLATMS, C_WORK_CPOT02, 15* + C_WORK_CPOT03, S_WORK_CLATMS, S_WORK_CLANHE, 16* + S_WORK_CPOT01, S_WORK_CPOT02, S_WORK_CPOT03 ) 17* 18* .. Scalar Arguments .. 19* INTEGER NN, NNS, NNT, NOUT 20* REAL THRESH 21* .. 22* .. Array Arguments .. 23* INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT ) 24* COMPLEX A( * ) 25* COMPLEX AINV( * ) 26* COMPLEX ASAV( * ) 27* COMPLEX B( * ) 28* COMPLEX BSAV( * ) 29* COMPLEX AFAC( * ) 30* COMPLEX ARF( * ) 31* COMPLEX ARFINV( * ) 32* COMPLEX XACT( * ) 33* COMPLEX X( * ) 34* COMPLEX C_WORK_CLATMS( * ) 35* COMPLEX C_WORK_CPOT02( * ) 36* COMPLEX C_WORK_CPOT03( * ) 37* REAL S_WORK_CLATMS( * ) 38* REAL S_WORK_CLANHE( * ) 39* REAL S_WORK_CPOT01( * ) 40* REAL S_WORK_CPOT02( * ) 41* REAL S_WORK_CPOT03( * ) 42* .. 43* 44* 45*> \par Purpose: 46* ============= 47*> 48*> \verbatim 49*> 50*> CDRVRFP tests the LAPACK RFP routines: 51*> CPFTRF, CPFTRS, and CPFTRI. 52*> 53*> This testing routine follow the same tests as CDRVPO (test for the full 54*> format Symmetric Positive Definite solver). 55*> 56*> The tests are performed in Full Format, conversion back and forth from 57*> full format to RFP format are performed using the routines CTRTTF and 58*> CTFTTR. 59*> 60*> First, a specific matrix A of size N is created. There is nine types of 61*> different matrixes possible. 62*> 1. Diagonal 6. Random, CNDNUM = sqrt(0.1/EPS) 63*> 2. Random, CNDNUM = 2 7. Random, CNDNUM = 0.1/EPS 64*> *3. First row and column zero 8. Scaled near underflow 65*> *4. Last row and column zero 9. Scaled near overflow 66*> *5. Middle row and column zero 67*> (* - tests error exits from CPFTRF, no test ratios are computed) 68*> A solution XACT of size N-by-NRHS is created and the associated right 69*> hand side B as well. Then CPFTRF is called to compute L (or U), the 70*> Cholesky factor of A. Then L (or U) is used to solve the linear system 71*> of equations AX = B. This gives X. Then L (or U) is used to compute the 72*> inverse of A, AINV. The following four tests are then performed: 73*> (1) norm( L*L' - A ) / ( N * norm(A) * EPS ) or 74*> norm( U'*U - A ) / ( N * norm(A) * EPS ), 75*> (2) norm(B - A*X) / ( norm(A) * norm(X) * EPS ), 76*> (3) norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ), 77*> (4) ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ), 78*> where EPS is the machine precision, RCOND the condition number of A, and 79*> norm( . ) the 1-norm for (1,2,3) and the inf-norm for (4). 80*> Errors occur when INFO parameter is not as expected. Failures occur when 81*> a test ratios is greater than THRES. 82*> \endverbatim 83* 84* Arguments: 85* ========== 86* 87*> \param[in] NOUT 88*> \verbatim 89*> NOUT is INTEGER 90*> The unit number for output. 91*> \endverbatim 92*> 93*> \param[in] NN 94*> \verbatim 95*> NN is INTEGER 96*> The number of values of N contained in the vector NVAL. 97*> \endverbatim 98*> 99*> \param[in] NVAL 100*> \verbatim 101*> NVAL is INTEGER array, dimension (NN) 102*> The values of the matrix dimension N. 103*> \endverbatim 104*> 105*> \param[in] NNS 106*> \verbatim 107*> NNS is INTEGER 108*> The number of values of NRHS contained in the vector NSVAL. 109*> \endverbatim 110*> 111*> \param[in] NSVAL 112*> \verbatim 113*> NSVAL is INTEGER array, dimension (NNS) 114*> The values of the number of right-hand sides NRHS. 115*> \endverbatim 116*> 117*> \param[in] NNT 118*> \verbatim 119*> NNT is INTEGER 120*> The number of values of MATRIX TYPE contained in the vector NTVAL. 121*> \endverbatim 122*> 123*> \param[in] NTVAL 124*> \verbatim 125*> NTVAL is INTEGER array, dimension (NNT) 126*> The values of matrix type (between 0 and 9 for PO/PP/PF matrices). 127*> \endverbatim 128*> 129*> \param[in] THRESH 130*> \verbatim 131*> THRESH is REAL 132*> The threshold value for the test ratios. A result is 133*> included in the output file if RESULT >= THRESH. To have 134*> every test ratio printed, use THRESH = 0. 135*> \endverbatim 136*> 137*> \param[out] A 138*> \verbatim 139*> A is COMPLEX array, dimension (NMAX*NMAX) 140*> \endverbatim 141*> 142*> \param[out] ASAV 143*> \verbatim 144*> ASAV is COMPLEX array, dimension (NMAX*NMAX) 145*> \endverbatim 146*> 147*> \param[out] AFAC 148*> \verbatim 149*> AFAC is COMPLEX array, dimension (NMAX*NMAX) 150*> \endverbatim 151*> 152*> \param[out] AINV 153*> \verbatim 154*> AINV is COMPLEX array, dimension (NMAX*NMAX) 155*> \endverbatim 156*> 157*> \param[out] B 158*> \verbatim 159*> B is COMPLEX array, dimension (NMAX*MAXRHS) 160*> \endverbatim 161*> 162*> \param[out] BSAV 163*> \verbatim 164*> BSAV is COMPLEX array, dimension (NMAX*MAXRHS) 165*> \endverbatim 166*> 167*> \param[out] XACT 168*> \verbatim 169*> XACT is COMPLEX array, dimension (NMAX*MAXRHS) 170*> \endverbatim 171*> 172*> \param[out] X 173*> \verbatim 174*> X is COMPLEX array, dimension (NMAX*MAXRHS) 175*> \endverbatim 176*> 177*> \param[out] ARF 178*> \verbatim 179*> ARF is COMPLEX array, dimension ((NMAX*(NMAX+1))/2) 180*> \endverbatim 181*> 182*> \param[out] ARFINV 183*> \verbatim 184*> ARFINV is COMPLEX array, dimension ((NMAX*(NMAX+1))/2) 185*> \endverbatim 186*> 187*> \param[out] C_WORK_CLATMS 188*> \verbatim 189*> C_WORK_CLATMS is COMPLEX array, dimension ( 3*NMAX ) 190*> \endverbatim 191*> 192*> \param[out] C_WORK_CPOT02 193*> \verbatim 194*> C_WORK_CPOT02 is COMPLEX array, dimension ( NMAX*MAXRHS ) 195*> \endverbatim 196*> 197*> \param[out] C_WORK_CPOT03 198*> \verbatim 199*> C_WORK_CPOT03 is COMPLEX array, dimension ( NMAX*NMAX ) 200*> \endverbatim 201*> 202*> \param[out] S_WORK_CLATMS 203*> \verbatim 204*> S_WORK_CLATMS is REAL array, dimension ( NMAX ) 205*> \endverbatim 206*> 207*> \param[out] S_WORK_CLANHE 208*> \verbatim 209*> S_WORK_CLANHE is REAL array, dimension ( NMAX ) 210*> \endverbatim 211*> 212*> \param[out] S_WORK_CPOT01 213*> \verbatim 214*> S_WORK_CPOT01 is REAL array, dimension ( NMAX ) 215*> \endverbatim 216*> 217*> \param[out] S_WORK_CPOT02 218*> \verbatim 219*> S_WORK_CPOT02 is REAL array, dimension ( NMAX ) 220*> \endverbatim 221*> 222*> \param[out] S_WORK_CPOT03 223*> \verbatim 224*> S_WORK_CPOT03 is REAL array, dimension ( NMAX ) 225*> \endverbatim 226* 227* Authors: 228* ======== 229* 230*> \author Univ. of Tennessee 231*> \author Univ. of California Berkeley 232*> \author Univ. of Colorado Denver 233*> \author NAG Ltd. 234* 235*> \ingroup complex_lin 236* 237* ===================================================================== 238 SUBROUTINE CDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL, 239 + THRESH, A, ASAV, AFAC, AINV, B, 240 + BSAV, XACT, X, ARF, ARFINV, 241 + C_WORK_CLATMS, C_WORK_CPOT02, 242 + C_WORK_CPOT03, S_WORK_CLATMS, S_WORK_CLANHE, 243 + S_WORK_CPOT01, S_WORK_CPOT02, S_WORK_CPOT03 ) 244* 245* -- LAPACK test routine -- 246* -- LAPACK is a software package provided by Univ. of Tennessee, -- 247* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 248* 249* .. Scalar Arguments .. 250 INTEGER NN, NNS, NNT, NOUT 251 REAL THRESH 252* .. 253* .. Array Arguments .. 254 INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT ) 255 COMPLEX A( * ) 256 COMPLEX AINV( * ) 257 COMPLEX ASAV( * ) 258 COMPLEX B( * ) 259 COMPLEX BSAV( * ) 260 COMPLEX AFAC( * ) 261 COMPLEX ARF( * ) 262 COMPLEX ARFINV( * ) 263 COMPLEX XACT( * ) 264 COMPLEX X( * ) 265 COMPLEX C_WORK_CLATMS( * ) 266 COMPLEX C_WORK_CPOT02( * ) 267 COMPLEX C_WORK_CPOT03( * ) 268 REAL S_WORK_CLATMS( * ) 269 REAL S_WORK_CLANHE( * ) 270 REAL S_WORK_CPOT01( * ) 271 REAL S_WORK_CPOT02( * ) 272 REAL S_WORK_CPOT03( * ) 273* .. 274* 275* ===================================================================== 276* 277* .. Parameters .. 278 REAL ONE, ZERO 279 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 280 INTEGER NTESTS 281 PARAMETER ( NTESTS = 4 ) 282* .. 283* .. Local Scalars .. 284 LOGICAL ZEROT 285 INTEGER I, INFO, IUPLO, LDA, LDB, IMAT, NERRS, NFAIL, 286 + NRHS, NRUN, IZERO, IOFF, K, NT, N, IFORM, IIN, 287 + IIT, IIS 288 CHARACTER DIST, CTYPE, UPLO, CFORM 289 INTEGER KL, KU, MODE 290 REAL ANORM, AINVNM, CNDNUM, RCONDC 291* .. 292* .. Local Arrays .. 293 CHARACTER UPLOS( 2 ), FORMS( 2 ) 294 INTEGER ISEED( 4 ), ISEEDY( 4 ) 295 REAL RESULT( NTESTS ) 296* .. 297* .. External Functions .. 298 REAL CLANHE 299 EXTERNAL CLANHE 300* .. 301* .. External Subroutines .. 302 EXTERNAL ALADHD, ALAERH, ALASVM, CGET04, CTFTTR, CLACPY, 303 + CLAIPD, CLARHS, CLATB4, CLATMS, CPFTRI, CPFTRF, 304 + CPFTRS, CPOT01, CPOT02, CPOT03, CPOTRI, CPOTRF, 305 + CTRTTF 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', 'C' / 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.GE.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 CLATB4 and generate a test 363* matrix with CLATMS. 364* 365 CALL CLATB4( 'CPO', IMAT, N, N, CTYPE, KL, KU, 366 + ANORM, MODE, CNDNUM, DIST ) 367* 368 SRNAMT = 'CLATMS' 369 CALL CLATMS( N, N, DIST, ISEED, CTYPE, 370 + S_WORK_CLATMS, 371 + MODE, CNDNUM, ANORM, KL, KU, UPLO, A, 372 + LDA, C_WORK_CLATMS, INFO ) 373* 374* Check error code from CLATMS. 375* 376 IF( INFO.NE.0 ) THEN 377 CALL ALAERH( 'CPF', 'CLATMS', 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* Set the imaginary part of the diagonals. 424* 425 CALL CLAIPD( N, A, LDA+1, 0 ) 426* 427* Save a copy of the matrix A in ASAV. 428* 429 CALL CLACPY( UPLO, N, N, A, LDA, ASAV, LDA ) 430* 431* Compute the condition number of A (RCONDC). 432* 433 IF( ZEROT ) THEN 434 RCONDC = ZERO 435 ELSE 436* 437* Compute the 1-norm of A. 438* 439 ANORM = CLANHE( '1', UPLO, N, A, LDA, 440 + S_WORK_CLANHE ) 441* 442* Factor the matrix A. 443* 444 CALL CPOTRF( UPLO, N, A, LDA, INFO ) 445* 446* Form the inverse of A. 447* 448 CALL CPOTRI( UPLO, N, A, LDA, INFO ) 449 450 IF ( N .NE. 0 ) THEN 451* 452* Compute the 1-norm condition number of A. 453* 454 AINVNM = CLANHE( '1', UPLO, N, A, LDA, 455 + S_WORK_CLANHE ) 456 RCONDC = ( ONE / ANORM ) / AINVNM 457* 458* Restore the matrix A. 459* 460 CALL CLACPY( UPLO, N, N, ASAV, LDA, A, LDA ) 461 END IF 462* 463 END IF 464* 465* Form an exact solution and set the right hand side. 466* 467 SRNAMT = 'CLARHS' 468 CALL CLARHS( 'CPO', 'N', UPLO, ' ', N, N, KL, KU, 469 + NRHS, A, LDA, XACT, LDA, B, LDA, 470 + ISEED, INFO ) 471 CALL CLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA ) 472* 473* Compute the L*L' or U'*U factorization of the 474* matrix and solve the system. 475* 476 CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) 477 CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDB ) 478* 479 SRNAMT = 'CTRTTF' 480 CALL CTRTTF( CFORM, UPLO, N, AFAC, LDA, ARF, INFO ) 481 SRNAMT = 'CPFTRF' 482 CALL CPFTRF( CFORM, UPLO, N, ARF, INFO ) 483* 484* Check error code from CPFTRF. 485* 486 IF( INFO.NE.IZERO ) THEN 487* 488* LANGOU: there is a small hick here: IZERO should 489* always be INFO however if INFO is ZERO, ALAERH does not 490* complain. 491* 492 CALL ALAERH( 'CPF', 'CPFSV ', INFO, IZERO, 493 + UPLO, N, N, -1, -1, NRHS, IIT, 494 + NFAIL, NERRS, NOUT ) 495 GO TO 100 496 END IF 497* 498* Skip the tests if INFO is not 0. 499* 500 IF( INFO.NE.0 ) THEN 501 GO TO 100 502 END IF 503* 504 SRNAMT = 'CPFTRS' 505 CALL CPFTRS( CFORM, UPLO, N, NRHS, ARF, X, LDB, 506 + INFO ) 507* 508 SRNAMT = 'CTFTTR' 509 CALL CTFTTR( CFORM, UPLO, N, ARF, AFAC, LDA, INFO ) 510* 511* Reconstruct matrix from factors and compute 512* residual. 513* 514 CALL CLACPY( UPLO, N, N, AFAC, LDA, ASAV, LDA ) 515 CALL CPOT01( UPLO, N, A, LDA, AFAC, LDA, 516 + S_WORK_CPOT01, RESULT( 1 ) ) 517 CALL CLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA ) 518* 519* Form the inverse and compute the residual. 520* 521 IF(MOD(N,2).EQ.0)THEN 522 CALL CLACPY( 'A', N+1, N/2, ARF, N+1, ARFINV, 523 + N+1 ) 524 ELSE 525 CALL CLACPY( 'A', N, (N+1)/2, ARF, N, ARFINV, 526 + N ) 527 END IF 528* 529 SRNAMT = 'CPFTRI' 530 CALL CPFTRI( CFORM, UPLO, N, ARFINV , INFO ) 531* 532 SRNAMT = 'CTFTTR' 533 CALL CTFTTR( CFORM, UPLO, N, ARFINV, AINV, LDA, 534 + INFO ) 535* 536* Check error code from CPFTRI. 537* 538 IF( INFO.NE.0 ) 539 + CALL ALAERH( 'CPO', 'CPFTRI', INFO, 0, UPLO, N, 540 + N, -1, -1, -1, IMAT, NFAIL, NERRS, 541 + NOUT ) 542* 543 CALL CPOT03( UPLO, N, A, LDA, AINV, LDA, 544 + C_WORK_CPOT03, LDA, S_WORK_CPOT03, 545 + RCONDC, RESULT( 2 ) ) 546* 547* Compute residual of the computed solution. 548* 549 CALL CLACPY( 'Full', N, NRHS, B, LDA, 550 + C_WORK_CPOT02, LDA ) 551 CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA, 552 + C_WORK_CPOT02, LDA, S_WORK_CPOT02, 553 + RESULT( 3 ) ) 554* 555* Check solution from generated exact solution. 556* 557 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 558 + RESULT( 4 ) ) 559 NT = 4 560* 561* Print information about the tests that did not 562* pass the threshold. 563* 564 DO 60 K = 1, NT 565 IF( RESULT( K ).GE.THRESH ) THEN 566 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 567 + CALL ALADHD( NOUT, 'CPF' ) 568 WRITE( NOUT, FMT = 9999 )'CPFSV ', UPLO, 569 + N, IIT, K, RESULT( K ) 570 NFAIL = NFAIL + 1 571 END IF 572 60 CONTINUE 573 NRUN = NRUN + NT 574 100 CONTINUE 575 110 CONTINUE 576 120 CONTINUE 577 980 CONTINUE 578 130 CONTINUE 579* 580* Print a summary of the results. 581* 582 CALL ALASVM( 'CPF', NOUT, NFAIL, NRUN, NERRS ) 583* 584 9999 FORMAT( 1X, A6, ', UPLO=''', A1, ''', N =', I5, ', type ', I1, 585 + ', test(', I1, ')=', G12.5 ) 586* 587 RETURN 588* 589* End of CDRVRFP 590* 591 END 592