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*> \date December 2016 236* 237*> \ingroup complex_lin 238* 239* ===================================================================== 240 SUBROUTINE CDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL, 241 + THRESH, A, ASAV, AFAC, AINV, B, 242 + BSAV, XACT, X, ARF, ARFINV, 243 + C_WORK_CLATMS, C_WORK_CPOT02, 244 + C_WORK_CPOT03, S_WORK_CLATMS, S_WORK_CLANHE, 245 + S_WORK_CPOT01, S_WORK_CPOT02, S_WORK_CPOT03 ) 246* 247* -- LAPACK test routine (version 3.7.0) -- 248* -- LAPACK is a software package provided by Univ. of Tennessee, -- 249* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 250* December 2016 251* 252* .. Scalar Arguments .. 253 INTEGER NN, NNS, NNT, NOUT 254 REAL THRESH 255* .. 256* .. Array Arguments .. 257 INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT ) 258 COMPLEX A( * ) 259 COMPLEX AINV( * ) 260 COMPLEX ASAV( * ) 261 COMPLEX B( * ) 262 COMPLEX BSAV( * ) 263 COMPLEX AFAC( * ) 264 COMPLEX ARF( * ) 265 COMPLEX ARFINV( * ) 266 COMPLEX XACT( * ) 267 COMPLEX X( * ) 268 COMPLEX C_WORK_CLATMS( * ) 269 COMPLEX C_WORK_CPOT02( * ) 270 COMPLEX C_WORK_CPOT03( * ) 271 REAL S_WORK_CLATMS( * ) 272 REAL S_WORK_CLANHE( * ) 273 REAL S_WORK_CPOT01( * ) 274 REAL S_WORK_CPOT02( * ) 275 REAL S_WORK_CPOT03( * ) 276* .. 277* 278* ===================================================================== 279* 280* .. Parameters .. 281 REAL ONE, ZERO 282 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 283 INTEGER NTESTS 284 PARAMETER ( NTESTS = 4 ) 285* .. 286* .. Local Scalars .. 287 LOGICAL ZEROT 288 INTEGER I, INFO, IUPLO, LDA, LDB, IMAT, NERRS, NFAIL, 289 + NRHS, NRUN, IZERO, IOFF, K, NT, N, IFORM, IIN, 290 + IIT, IIS 291 CHARACTER DIST, CTYPE, UPLO, CFORM 292 INTEGER KL, KU, MODE 293 REAL ANORM, AINVNM, CNDNUM, RCONDC 294* .. 295* .. Local Arrays .. 296 CHARACTER UPLOS( 2 ), FORMS( 2 ) 297 INTEGER ISEED( 4 ), ISEEDY( 4 ) 298 REAL RESULT( NTESTS ) 299* .. 300* .. External Functions .. 301 REAL CLANHE 302 EXTERNAL CLANHE 303* .. 304* .. External Subroutines .. 305 EXTERNAL ALADHD, ALAERH, ALASVM, CGET04, CTFTTR, CLACPY, 306 + CLAIPD, CLARHS, CLATB4, CLATMS, CPFTRI, CPFTRF, 307 + CPFTRS, CPOT01, CPOT02, CPOT03, CPOTRI, CPOTRF, 308 + CTRTTF 309* .. 310* .. Scalars in Common .. 311 CHARACTER*32 SRNAMT 312* .. 313* .. Common blocks .. 314 COMMON / SRNAMC / SRNAMT 315* .. 316* .. Data statements .. 317 DATA ISEEDY / 1988, 1989, 1990, 1991 / 318 DATA UPLOS / 'U', 'L' / 319 DATA FORMS / 'N', 'C' / 320* .. 321* .. Executable Statements .. 322* 323* Initialize constants and the random number seed. 324* 325 NRUN = 0 326 NFAIL = 0 327 NERRS = 0 328 DO 10 I = 1, 4 329 ISEED( I ) = ISEEDY( I ) 330 10 CONTINUE 331* 332 DO 130 IIN = 1, NN 333* 334 N = NVAL( IIN ) 335 LDA = MAX( N, 1 ) 336 LDB = MAX( N, 1 ) 337* 338 DO 980 IIS = 1, NNS 339* 340 NRHS = NSVAL( IIS ) 341* 342 DO 120 IIT = 1, NNT 343* 344 IMAT = NTVAL( IIT ) 345* 346* If N.EQ.0, only consider the first type 347* 348 IF( N.EQ.0 .AND. IIT.GE.1 ) GO TO 120 349* 350* Skip types 3, 4, or 5 if the matrix size is too small. 351* 352 IF( IMAT.EQ.4 .AND. N.LE.1 ) GO TO 120 353 IF( IMAT.EQ.5 .AND. N.LE.2 ) GO TO 120 354* 355* Do first for UPLO = 'U', then for UPLO = 'L' 356* 357 DO 110 IUPLO = 1, 2 358 UPLO = UPLOS( IUPLO ) 359* 360* Do first for CFORM = 'N', then for CFORM = 'C' 361* 362 DO 100 IFORM = 1, 2 363 CFORM = FORMS( IFORM ) 364* 365* Set up parameters with CLATB4 and generate a test 366* matrix with CLATMS. 367* 368 CALL CLATB4( 'CPO', IMAT, N, N, CTYPE, KL, KU, 369 + ANORM, MODE, CNDNUM, DIST ) 370* 371 SRNAMT = 'CLATMS' 372 CALL CLATMS( N, N, DIST, ISEED, CTYPE, 373 + S_WORK_CLATMS, 374 + MODE, CNDNUM, ANORM, KL, KU, UPLO, A, 375 + LDA, C_WORK_CLATMS, INFO ) 376* 377* Check error code from CLATMS. 378* 379 IF( INFO.NE.0 ) THEN 380 CALL ALAERH( 'CPF', 'CLATMS', INFO, 0, UPLO, N, 381 + N, -1, -1, -1, IIT, NFAIL, NERRS, 382 + NOUT ) 383 GO TO 100 384 END IF 385* 386* For types 3-5, zero one row and column of the matrix to 387* test that INFO is returned correctly. 388* 389 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5 390 IF( ZEROT ) THEN 391 IF( IIT.EQ.3 ) THEN 392 IZERO = 1 393 ELSE IF( IIT.EQ.4 ) THEN 394 IZERO = N 395 ELSE 396 IZERO = N / 2 + 1 397 END IF 398 IOFF = ( IZERO-1 )*LDA 399* 400* Set row and column IZERO of A to 0. 401* 402 IF( IUPLO.EQ.1 ) THEN 403 DO 20 I = 1, IZERO - 1 404 A( IOFF+I ) = ZERO 405 20 CONTINUE 406 IOFF = IOFF + IZERO 407 DO 30 I = IZERO, N 408 A( IOFF ) = ZERO 409 IOFF = IOFF + LDA 410 30 CONTINUE 411 ELSE 412 IOFF = IZERO 413 DO 40 I = 1, IZERO - 1 414 A( IOFF ) = ZERO 415 IOFF = IOFF + LDA 416 40 CONTINUE 417 IOFF = IOFF - IZERO 418 DO 50 I = IZERO, N 419 A( IOFF+I ) = ZERO 420 50 CONTINUE 421 END IF 422 ELSE 423 IZERO = 0 424 END IF 425* 426* Set the imaginary part of the diagonals. 427* 428 CALL CLAIPD( N, A, LDA+1, 0 ) 429* 430* Save a copy of the matrix A in ASAV. 431* 432 CALL CLACPY( UPLO, N, N, A, LDA, ASAV, LDA ) 433* 434* Compute the condition number of A (RCONDC). 435* 436 IF( ZEROT ) THEN 437 RCONDC = ZERO 438 ELSE 439* 440* Compute the 1-norm of A. 441* 442 ANORM = CLANHE( '1', UPLO, N, A, LDA, 443 + S_WORK_CLANHE ) 444* 445* Factor the matrix A. 446* 447 CALL CPOTRF( UPLO, N, A, LDA, INFO ) 448* 449* Form the inverse of A. 450* 451 CALL CPOTRI( UPLO, N, A, LDA, INFO ) 452 453 IF ( N .NE. 0 ) THEN 454* 455* Compute the 1-norm condition number of A. 456* 457 AINVNM = CLANHE( '1', UPLO, N, A, LDA, 458 + S_WORK_CLANHE ) 459 RCONDC = ( ONE / ANORM ) / AINVNM 460* 461* Restore the matrix A. 462* 463 CALL CLACPY( UPLO, N, N, ASAV, LDA, A, LDA ) 464 END IF 465* 466 END IF 467* 468* Form an exact solution and set the right hand side. 469* 470 SRNAMT = 'CLARHS' 471 CALL CLARHS( 'CPO', 'N', UPLO, ' ', N, N, KL, KU, 472 + NRHS, A, LDA, XACT, LDA, B, LDA, 473 + ISEED, INFO ) 474 CALL CLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA ) 475* 476* Compute the L*L' or U'*U factorization of the 477* matrix and solve the system. 478* 479 CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) 480 CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDB ) 481* 482 SRNAMT = 'CTRTTF' 483 CALL CTRTTF( CFORM, UPLO, N, AFAC, LDA, ARF, INFO ) 484 SRNAMT = 'CPFTRF' 485 CALL CPFTRF( CFORM, UPLO, N, ARF, INFO ) 486* 487* Check error code from CPFTRF. 488* 489 IF( INFO.NE.IZERO ) THEN 490* 491* LANGOU: there is a small hick here: IZERO should 492* always be INFO however if INFO is ZERO, ALAERH does not 493* complain. 494* 495 CALL ALAERH( 'CPF', 'CPFSV ', INFO, IZERO, 496 + UPLO, N, N, -1, -1, NRHS, IIT, 497 + NFAIL, NERRS, NOUT ) 498 GO TO 100 499 END IF 500* 501* Skip the tests if INFO is not 0. 502* 503 IF( INFO.NE.0 ) THEN 504 GO TO 100 505 END IF 506* 507 SRNAMT = 'CPFTRS' 508 CALL CPFTRS( CFORM, UPLO, N, NRHS, ARF, X, LDB, 509 + INFO ) 510* 511 SRNAMT = 'CTFTTR' 512 CALL CTFTTR( CFORM, UPLO, N, ARF, AFAC, LDA, INFO ) 513* 514* Reconstruct matrix from factors and compute 515* residual. 516* 517 CALL CLACPY( UPLO, N, N, AFAC, LDA, ASAV, LDA ) 518 CALL CPOT01( UPLO, N, A, LDA, AFAC, LDA, 519 + S_WORK_CPOT01, RESULT( 1 ) ) 520 CALL CLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA ) 521* 522* Form the inverse and compute the residual. 523* 524 IF(MOD(N,2).EQ.0)THEN 525 CALL CLACPY( 'A', N+1, N/2, ARF, N+1, ARFINV, 526 + N+1 ) 527 ELSE 528 CALL CLACPY( 'A', N, (N+1)/2, ARF, N, ARFINV, 529 + N ) 530 END IF 531* 532 SRNAMT = 'CPFTRI' 533 CALL CPFTRI( CFORM, UPLO, N, ARFINV , INFO ) 534* 535 SRNAMT = 'CTFTTR' 536 CALL CTFTTR( CFORM, UPLO, N, ARFINV, AINV, LDA, 537 + INFO ) 538* 539* Check error code from CPFTRI. 540* 541 IF( INFO.NE.0 ) 542 + CALL ALAERH( 'CPO', 'CPFTRI', INFO, 0, UPLO, N, 543 + N, -1, -1, -1, IMAT, NFAIL, NERRS, 544 + NOUT ) 545* 546 CALL CPOT03( UPLO, N, A, LDA, AINV, LDA, 547 + C_WORK_CPOT03, LDA, S_WORK_CPOT03, 548 + RCONDC, RESULT( 2 ) ) 549* 550* Compute residual of the computed solution. 551* 552 CALL CLACPY( 'Full', N, NRHS, B, LDA, 553 + C_WORK_CPOT02, LDA ) 554 CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA, 555 + C_WORK_CPOT02, LDA, S_WORK_CPOT02, 556 + RESULT( 3 ) ) 557* 558* Check solution from generated exact solution. 559* 560 CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, 561 + RESULT( 4 ) ) 562 NT = 4 563* 564* Print information about the tests that did not 565* pass the threshold. 566* 567 DO 60 K = 1, NT 568 IF( RESULT( K ).GE.THRESH ) THEN 569 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) 570 + CALL ALADHD( NOUT, 'CPF' ) 571 WRITE( NOUT, FMT = 9999 )'CPFSV ', UPLO, 572 + N, IIT, K, RESULT( K ) 573 NFAIL = NFAIL + 1 574 END IF 575 60 CONTINUE 576 NRUN = NRUN + NT 577 100 CONTINUE 578 110 CONTINUE 579 120 CONTINUE 580 980 CONTINUE 581 130 CONTINUE 582* 583* Print a summary of the results. 584* 585 CALL ALASVM( 'CPF', NOUT, NFAIL, NRUN, NERRS ) 586* 587 9999 FORMAT( 1X, A6, ', UPLO=''', A1, ''', N =', I5, ', type ', I1, 588 + ', test(', I1, ')=', G12.5 ) 589* 590 RETURN 591* 592* End of CDRVRFP 593* 594 END 595