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