1*> \brief \b CLATTP 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 CLATTP( IMAT, UPLO, TRANS, DIAG, ISEED, N, AP, B, WORK, 12* RWORK, INFO ) 13* 14* .. Scalar Arguments .. 15* CHARACTER DIAG, TRANS, UPLO 16* INTEGER IMAT, INFO, N 17* .. 18* .. Array Arguments .. 19* INTEGER ISEED( 4 ) 20* REAL RWORK( * ) 21* COMPLEX AP( * ), B( * ), WORK( * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> CLATTP generates a triangular test matrix in packed storage. 31*> IMAT and UPLO uniquely specify the properties of the test matrix, 32*> which is returned in the array AP. 33*> \endverbatim 34* 35* Arguments: 36* ========== 37* 38*> \param[in] IMAT 39*> \verbatim 40*> IMAT is INTEGER 41*> An integer key describing which matrix to generate for this 42*> path. 43*> \endverbatim 44*> 45*> \param[in] UPLO 46*> \verbatim 47*> UPLO is CHARACTER*1 48*> Specifies whether the matrix A will be upper or lower 49*> triangular. 50*> = 'U': Upper triangular 51*> = 'L': Lower triangular 52*> \endverbatim 53*> 54*> \param[in] TRANS 55*> \verbatim 56*> TRANS is CHARACTER*1 57*> Specifies whether the matrix or its transpose will be used. 58*> = 'N': No transpose 59*> = 'T': Transpose 60*> = 'C': Conjugate transpose 61*> \endverbatim 62*> 63*> \param[out] DIAG 64*> \verbatim 65*> DIAG is CHARACTER*1 66*> Specifies whether or not the matrix A is unit triangular. 67*> = 'N': Non-unit triangular 68*> = 'U': Unit triangular 69*> \endverbatim 70*> 71*> \param[in,out] ISEED 72*> \verbatim 73*> ISEED is INTEGER array, dimension (4) 74*> The seed vector for the random number generator (used in 75*> CLATMS). Modified on exit. 76*> \endverbatim 77*> 78*> \param[in] N 79*> \verbatim 80*> N is INTEGER 81*> The order of the matrix to be generated. 82*> \endverbatim 83*> 84*> \param[out] AP 85*> \verbatim 86*> AP is COMPLEX array, dimension (N*(N+1)/2) 87*> The upper or lower triangular matrix A, packed columnwise in 88*> a linear array. The j-th column of A is stored in the array 89*> AP as follows: 90*> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j; 91*> if UPLO = 'L', 92*> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n. 93*> \endverbatim 94*> 95*> \param[out] B 96*> \verbatim 97*> B is COMPLEX array, dimension (N) 98*> The right hand side vector, if IMAT > 10. 99*> \endverbatim 100*> 101*> \param[out] WORK 102*> \verbatim 103*> WORK is COMPLEX array, dimension (2*N) 104*> \endverbatim 105*> 106*> \param[out] RWORK 107*> \verbatim 108*> RWORK is REAL array, dimension (N) 109*> \endverbatim 110*> 111*> \param[out] INFO 112*> \verbatim 113*> INFO is INTEGER 114*> = 0: successful exit 115*> < 0: if INFO = -i, the i-th argument had an illegal value 116*> \endverbatim 117* 118* Authors: 119* ======== 120* 121*> \author Univ. of Tennessee 122*> \author Univ. of California Berkeley 123*> \author Univ. of Colorado Denver 124*> \author NAG Ltd. 125* 126*> \ingroup complex_lin 127* 128* ===================================================================== 129 SUBROUTINE CLATTP( IMAT, UPLO, TRANS, DIAG, ISEED, N, AP, B, WORK, 130 $ RWORK, INFO ) 131* 132* -- LAPACK test routine -- 133* -- LAPACK is a software package provided by Univ. of Tennessee, -- 134* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 135* 136* .. Scalar Arguments .. 137 CHARACTER DIAG, TRANS, UPLO 138 INTEGER IMAT, INFO, N 139* .. 140* .. Array Arguments .. 141 INTEGER ISEED( 4 ) 142 REAL RWORK( * ) 143 COMPLEX AP( * ), B( * ), WORK( * ) 144* .. 145* 146* ===================================================================== 147* 148* .. Parameters .. 149 REAL ONE, TWO, ZERO 150 PARAMETER ( ONE = 1.0E+0, TWO = 2.0E+0, ZERO = 0.0E+0 ) 151* .. 152* .. Local Scalars .. 153 LOGICAL UPPER 154 CHARACTER DIST, PACKIT, TYPE 155 CHARACTER*3 PATH 156 INTEGER I, IY, J, JC, JCNEXT, JCOUNT, JJ, JL, JR, JX, 157 $ KL, KU, MODE 158 REAL ANORM, BIGNUM, BNORM, BSCAL, C, CNDNUM, REXP, 159 $ SFAC, SMLNUM, T, TEXP, TLEFT, TSCAL, ULP, UNFL, 160 $ X, Y, Z 161 COMPLEX CTEMP, PLUS1, PLUS2, RA, RB, S, STAR1 162* .. 163* .. External Functions .. 164 LOGICAL LSAME 165 INTEGER ICAMAX 166 REAL SLAMCH 167 COMPLEX CLARND 168 EXTERNAL LSAME, ICAMAX, SLAMCH, CLARND 169* .. 170* .. External Subroutines .. 171 EXTERNAL CLARNV, CLATB4, CLATMS, CROT, CROTG, CSSCAL, 172 $ SLABAD, SLARNV 173* .. 174* .. Intrinsic Functions .. 175 INTRINSIC ABS, CMPLX, CONJG, MAX, REAL, SQRT 176* .. 177* .. Executable Statements .. 178* 179 PATH( 1: 1 ) = 'Complex precision' 180 PATH( 2: 3 ) = 'TP' 181 UNFL = SLAMCH( 'Safe minimum' ) 182 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 183 SMLNUM = UNFL 184 BIGNUM = ( ONE-ULP ) / SMLNUM 185 CALL SLABAD( SMLNUM, BIGNUM ) 186 IF( ( IMAT.GE.7 .AND. IMAT.LE.10 ) .OR. IMAT.EQ.18 ) THEN 187 DIAG = 'U' 188 ELSE 189 DIAG = 'N' 190 END IF 191 INFO = 0 192* 193* Quick return if N.LE.0. 194* 195 IF( N.LE.0 ) 196 $ RETURN 197* 198* Call CLATB4 to set parameters for CLATMS. 199* 200 UPPER = LSAME( UPLO, 'U' ) 201 IF( UPPER ) THEN 202 CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 203 $ CNDNUM, DIST ) 204 PACKIT = 'C' 205 ELSE 206 CALL CLATB4( PATH, -IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 207 $ CNDNUM, DIST ) 208 PACKIT = 'R' 209 END IF 210* 211* IMAT <= 6: Non-unit triangular matrix 212* 213 IF( IMAT.LE.6 ) THEN 214 CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, CNDNUM, 215 $ ANORM, KL, KU, PACKIT, AP, N, WORK, INFO ) 216* 217* IMAT > 6: Unit triangular matrix 218* The diagonal is deliberately set to something other than 1. 219* 220* IMAT = 7: Matrix is the identity 221* 222 ELSE IF( IMAT.EQ.7 ) THEN 223 IF( UPPER ) THEN 224 JC = 1 225 DO 20 J = 1, N 226 DO 10 I = 1, J - 1 227 AP( JC+I-1 ) = ZERO 228 10 CONTINUE 229 AP( JC+J-1 ) = J 230 JC = JC + J 231 20 CONTINUE 232 ELSE 233 JC = 1 234 DO 40 J = 1, N 235 AP( JC ) = J 236 DO 30 I = J + 1, N 237 AP( JC+I-J ) = ZERO 238 30 CONTINUE 239 JC = JC + N - J + 1 240 40 CONTINUE 241 END IF 242* 243* IMAT > 7: Non-trivial unit triangular matrix 244* 245* Generate a unit triangular matrix T with condition CNDNUM by 246* forming a triangular matrix with known singular values and 247* filling in the zero entries with Givens rotations. 248* 249 ELSE IF( IMAT.LE.10 ) THEN 250 IF( UPPER ) THEN 251 JC = 0 252 DO 60 J = 1, N 253 DO 50 I = 1, J - 1 254 AP( JC+I ) = ZERO 255 50 CONTINUE 256 AP( JC+J ) = J 257 JC = JC + J 258 60 CONTINUE 259 ELSE 260 JC = 1 261 DO 80 J = 1, N 262 AP( JC ) = J 263 DO 70 I = J + 1, N 264 AP( JC+I-J ) = ZERO 265 70 CONTINUE 266 JC = JC + N - J + 1 267 80 CONTINUE 268 END IF 269* 270* Since the trace of a unit triangular matrix is 1, the product 271* of its singular values must be 1. Let s = sqrt(CNDNUM), 272* x = sqrt(s) - 1/sqrt(s), y = sqrt(2/(n-2))*x, and z = x**2. 273* The following triangular matrix has singular values s, 1, 1, 274* ..., 1, 1/s: 275* 276* 1 y y y ... y y z 277* 1 0 0 ... 0 0 y 278* 1 0 ... 0 0 y 279* . ... . . . 280* . . . . 281* 1 0 y 282* 1 y 283* 1 284* 285* To fill in the zeros, we first multiply by a matrix with small 286* condition number of the form 287* 288* 1 0 0 0 0 ... 289* 1 + * 0 0 ... 290* 1 + 0 0 0 291* 1 + * 0 0 292* 1 + 0 0 293* ... 294* 1 + 0 295* 1 0 296* 1 297* 298* Each element marked with a '*' is formed by taking the product 299* of the adjacent elements marked with '+'. The '*'s can be 300* chosen freely, and the '+'s are chosen so that the inverse of 301* T will have elements of the same magnitude as T. If the *'s in 302* both T and inv(T) have small magnitude, T is well conditioned. 303* The two offdiagonals of T are stored in WORK. 304* 305* The product of these two matrices has the form 306* 307* 1 y y y y y . y y z 308* 1 + * 0 0 . 0 0 y 309* 1 + 0 0 . 0 0 y 310* 1 + * . . . . 311* 1 + . . . . 312* . . . . . 313* . . . . 314* 1 + y 315* 1 y 316* 1 317* 318* Now we multiply by Givens rotations, using the fact that 319* 320* [ c s ] [ 1 w ] [ -c -s ] = [ 1 -w ] 321* [ -s c ] [ 0 1 ] [ s -c ] [ 0 1 ] 322* and 323* [ -c -s ] [ 1 0 ] [ c s ] = [ 1 0 ] 324* [ s -c ] [ w 1 ] [ -s c ] [ -w 1 ] 325* 326* where c = w / sqrt(w**2+4) and s = 2 / sqrt(w**2+4). 327* 328 STAR1 = 0.25*CLARND( 5, ISEED ) 329 SFAC = 0.5 330 PLUS1 = SFAC*CLARND( 5, ISEED ) 331 DO 90 J = 1, N, 2 332 PLUS2 = STAR1 / PLUS1 333 WORK( J ) = PLUS1 334 WORK( N+J ) = STAR1 335 IF( J+1.LE.N ) THEN 336 WORK( J+1 ) = PLUS2 337 WORK( N+J+1 ) = ZERO 338 PLUS1 = STAR1 / PLUS2 339 REXP = CLARND( 2, ISEED ) 340 IF( REXP.LT.ZERO ) THEN 341 STAR1 = -SFAC**( ONE-REXP )*CLARND( 5, ISEED ) 342 ELSE 343 STAR1 = SFAC**( ONE+REXP )*CLARND( 5, ISEED ) 344 END IF 345 END IF 346 90 CONTINUE 347* 348 X = SQRT( CNDNUM ) - ONE / SQRT( CNDNUM ) 349 IF( N.GT.2 ) THEN 350 Y = SQRT( TWO / REAL( N-2 ) )*X 351 ELSE 352 Y = ZERO 353 END IF 354 Z = X*X 355* 356 IF( UPPER ) THEN 357* 358* Set the upper triangle of A with a unit triangular matrix 359* of known condition number. 360* 361 JC = 1 362 DO 100 J = 2, N 363 AP( JC+1 ) = Y 364 IF( J.GT.2 ) 365 $ AP( JC+J-1 ) = WORK( J-2 ) 366 IF( J.GT.3 ) 367 $ AP( JC+J-2 ) = WORK( N+J-3 ) 368 JC = JC + J 369 100 CONTINUE 370 JC = JC - N 371 AP( JC+1 ) = Z 372 DO 110 J = 2, N - 1 373 AP( JC+J ) = Y 374 110 CONTINUE 375 ELSE 376* 377* Set the lower triangle of A with a unit triangular matrix 378* of known condition number. 379* 380 DO 120 I = 2, N - 1 381 AP( I ) = Y 382 120 CONTINUE 383 AP( N ) = Z 384 JC = N + 1 385 DO 130 J = 2, N - 1 386 AP( JC+1 ) = WORK( J-1 ) 387 IF( J.LT.N-1 ) 388 $ AP( JC+2 ) = WORK( N+J-1 ) 389 AP( JC+N-J ) = Y 390 JC = JC + N - J + 1 391 130 CONTINUE 392 END IF 393* 394* Fill in the zeros using Givens rotations 395* 396 IF( UPPER ) THEN 397 JC = 1 398 DO 150 J = 1, N - 1 399 JCNEXT = JC + J 400 RA = AP( JCNEXT+J-1 ) 401 RB = TWO 402 CALL CROTG( RA, RB, C, S ) 403* 404* Multiply by [ c s; -conjg(s) c] on the left. 405* 406 IF( N.GT.J+1 ) THEN 407 JX = JCNEXT + J 408 DO 140 I = J + 2, N 409 CTEMP = C*AP( JX+J ) + S*AP( JX+J+1 ) 410 AP( JX+J+1 ) = -CONJG( S )*AP( JX+J ) + 411 $ C*AP( JX+J+1 ) 412 AP( JX+J ) = CTEMP 413 JX = JX + I 414 140 CONTINUE 415 END IF 416* 417* Multiply by [-c -s; conjg(s) -c] on the right. 418* 419 IF( J.GT.1 ) 420 $ CALL CROT( J-1, AP( JCNEXT ), 1, AP( JC ), 1, -C, -S ) 421* 422* Negate A(J,J+1). 423* 424 AP( JCNEXT+J-1 ) = -AP( JCNEXT+J-1 ) 425 JC = JCNEXT 426 150 CONTINUE 427 ELSE 428 JC = 1 429 DO 170 J = 1, N - 1 430 JCNEXT = JC + N - J + 1 431 RA = AP( JC+1 ) 432 RB = TWO 433 CALL CROTG( RA, RB, C, S ) 434 S = CONJG( S ) 435* 436* Multiply by [ c -s; conjg(s) c] on the right. 437* 438 IF( N.GT.J+1 ) 439 $ CALL CROT( N-J-1, AP( JCNEXT+1 ), 1, AP( JC+2 ), 1, C, 440 $ -S ) 441* 442* Multiply by [-c s; -conjg(s) -c] on the left. 443* 444 IF( J.GT.1 ) THEN 445 JX = 1 446 DO 160 I = 1, J - 1 447 CTEMP = -C*AP( JX+J-I ) + S*AP( JX+J-I+1 ) 448 AP( JX+J-I+1 ) = -CONJG( S )*AP( JX+J-I ) - 449 $ C*AP( JX+J-I+1 ) 450 AP( JX+J-I ) = CTEMP 451 JX = JX + N - I + 1 452 160 CONTINUE 453 END IF 454* 455* Negate A(J+1,J). 456* 457 AP( JC+1 ) = -AP( JC+1 ) 458 JC = JCNEXT 459 170 CONTINUE 460 END IF 461* 462* IMAT > 10: Pathological test cases. These triangular matrices 463* are badly scaled or badly conditioned, so when used in solving a 464* triangular system they may cause overflow in the solution vector. 465* 466 ELSE IF( IMAT.EQ.11 ) THEN 467* 468* Type 11: Generate a triangular matrix with elements between 469* -1 and 1. Give the diagonal norm 2 to make it well-conditioned. 470* Make the right hand side large so that it requires scaling. 471* 472 IF( UPPER ) THEN 473 JC = 1 474 DO 180 J = 1, N 475 CALL CLARNV( 4, ISEED, J-1, AP( JC ) ) 476 AP( JC+J-1 ) = CLARND( 5, ISEED )*TWO 477 JC = JC + J 478 180 CONTINUE 479 ELSE 480 JC = 1 481 DO 190 J = 1, N 482 IF( J.LT.N ) 483 $ CALL CLARNV( 4, ISEED, N-J, AP( JC+1 ) ) 484 AP( JC ) = CLARND( 5, ISEED )*TWO 485 JC = JC + N - J + 1 486 190 CONTINUE 487 END IF 488* 489* Set the right hand side so that the largest value is BIGNUM. 490* 491 CALL CLARNV( 2, ISEED, N, B ) 492 IY = ICAMAX( N, B, 1 ) 493 BNORM = ABS( B( IY ) ) 494 BSCAL = BIGNUM / MAX( ONE, BNORM ) 495 CALL CSSCAL( N, BSCAL, B, 1 ) 496* 497 ELSE IF( IMAT.EQ.12 ) THEN 498* 499* Type 12: Make the first diagonal element in the solve small to 500* cause immediate overflow when dividing by T(j,j). 501* In type 12, the offdiagonal elements are small (CNORM(j) < 1). 502* 503 CALL CLARNV( 2, ISEED, N, B ) 504 TSCAL = ONE / MAX( ONE, REAL( N-1 ) ) 505 IF( UPPER ) THEN 506 JC = 1 507 DO 200 J = 1, N 508 CALL CLARNV( 4, ISEED, J-1, AP( JC ) ) 509 CALL CSSCAL( J-1, TSCAL, AP( JC ), 1 ) 510 AP( JC+J-1 ) = CLARND( 5, ISEED ) 511 JC = JC + J 512 200 CONTINUE 513 AP( N*( N+1 ) / 2 ) = SMLNUM*AP( N*( N+1 ) / 2 ) 514 ELSE 515 JC = 1 516 DO 210 J = 1, N 517 CALL CLARNV( 2, ISEED, N-J, AP( JC+1 ) ) 518 CALL CSSCAL( N-J, TSCAL, AP( JC+1 ), 1 ) 519 AP( JC ) = CLARND( 5, ISEED ) 520 JC = JC + N - J + 1 521 210 CONTINUE 522 AP( 1 ) = SMLNUM*AP( 1 ) 523 END IF 524* 525 ELSE IF( IMAT.EQ.13 ) THEN 526* 527* Type 13: Make the first diagonal element in the solve small to 528* cause immediate overflow when dividing by T(j,j). 529* In type 13, the offdiagonal elements are O(1) (CNORM(j) > 1). 530* 531 CALL CLARNV( 2, ISEED, N, B ) 532 IF( UPPER ) THEN 533 JC = 1 534 DO 220 J = 1, N 535 CALL CLARNV( 4, ISEED, J-1, AP( JC ) ) 536 AP( JC+J-1 ) = CLARND( 5, ISEED ) 537 JC = JC + J 538 220 CONTINUE 539 AP( N*( N+1 ) / 2 ) = SMLNUM*AP( N*( N+1 ) / 2 ) 540 ELSE 541 JC = 1 542 DO 230 J = 1, N 543 CALL CLARNV( 4, ISEED, N-J, AP( JC+1 ) ) 544 AP( JC ) = CLARND( 5, ISEED ) 545 JC = JC + N - J + 1 546 230 CONTINUE 547 AP( 1 ) = SMLNUM*AP( 1 ) 548 END IF 549* 550 ELSE IF( IMAT.EQ.14 ) THEN 551* 552* Type 14: T is diagonal with small numbers on the diagonal to 553* make the growth factor underflow, but a small right hand side 554* chosen so that the solution does not overflow. 555* 556 IF( UPPER ) THEN 557 JCOUNT = 1 558 JC = ( N-1 )*N / 2 + 1 559 DO 250 J = N, 1, -1 560 DO 240 I = 1, J - 1 561 AP( JC+I-1 ) = ZERO 562 240 CONTINUE 563 IF( JCOUNT.LE.2 ) THEN 564 AP( JC+J-1 ) = SMLNUM*CLARND( 5, ISEED ) 565 ELSE 566 AP( JC+J-1 ) = CLARND( 5, ISEED ) 567 END IF 568 JCOUNT = JCOUNT + 1 569 IF( JCOUNT.GT.4 ) 570 $ JCOUNT = 1 571 JC = JC - J + 1 572 250 CONTINUE 573 ELSE 574 JCOUNT = 1 575 JC = 1 576 DO 270 J = 1, N 577 DO 260 I = J + 1, N 578 AP( JC+I-J ) = ZERO 579 260 CONTINUE 580 IF( JCOUNT.LE.2 ) THEN 581 AP( JC ) = SMLNUM*CLARND( 5, ISEED ) 582 ELSE 583 AP( JC ) = CLARND( 5, ISEED ) 584 END IF 585 JCOUNT = JCOUNT + 1 586 IF( JCOUNT.GT.4 ) 587 $ JCOUNT = 1 588 JC = JC + N - J + 1 589 270 CONTINUE 590 END IF 591* 592* Set the right hand side alternately zero and small. 593* 594 IF( UPPER ) THEN 595 B( 1 ) = ZERO 596 DO 280 I = N, 2, -2 597 B( I ) = ZERO 598 B( I-1 ) = SMLNUM*CLARND( 5, ISEED ) 599 280 CONTINUE 600 ELSE 601 B( N ) = ZERO 602 DO 290 I = 1, N - 1, 2 603 B( I ) = ZERO 604 B( I+1 ) = SMLNUM*CLARND( 5, ISEED ) 605 290 CONTINUE 606 END IF 607* 608 ELSE IF( IMAT.EQ.15 ) THEN 609* 610* Type 15: Make the diagonal elements small to cause gradual 611* overflow when dividing by T(j,j). To control the amount of 612* scaling needed, the matrix is bidiagonal. 613* 614 TEXP = ONE / MAX( ONE, REAL( N-1 ) ) 615 TSCAL = SMLNUM**TEXP 616 CALL CLARNV( 4, ISEED, N, B ) 617 IF( UPPER ) THEN 618 JC = 1 619 DO 310 J = 1, N 620 DO 300 I = 1, J - 2 621 AP( JC+I-1 ) = ZERO 622 300 CONTINUE 623 IF( J.GT.1 ) 624 $ AP( JC+J-2 ) = CMPLX( -ONE, -ONE ) 625 AP( JC+J-1 ) = TSCAL*CLARND( 5, ISEED ) 626 JC = JC + J 627 310 CONTINUE 628 B( N ) = CMPLX( ONE, ONE ) 629 ELSE 630 JC = 1 631 DO 330 J = 1, N 632 DO 320 I = J + 2, N 633 AP( JC+I-J ) = ZERO 634 320 CONTINUE 635 IF( J.LT.N ) 636 $ AP( JC+1 ) = CMPLX( -ONE, -ONE ) 637 AP( JC ) = TSCAL*CLARND( 5, ISEED ) 638 JC = JC + N - J + 1 639 330 CONTINUE 640 B( 1 ) = CMPLX( ONE, ONE ) 641 END IF 642* 643 ELSE IF( IMAT.EQ.16 ) THEN 644* 645* Type 16: One zero diagonal element. 646* 647 IY = N / 2 + 1 648 IF( UPPER ) THEN 649 JC = 1 650 DO 340 J = 1, N 651 CALL CLARNV( 4, ISEED, J, AP( JC ) ) 652 IF( J.NE.IY ) THEN 653 AP( JC+J-1 ) = CLARND( 5, ISEED )*TWO 654 ELSE 655 AP( JC+J-1 ) = ZERO 656 END IF 657 JC = JC + J 658 340 CONTINUE 659 ELSE 660 JC = 1 661 DO 350 J = 1, N 662 CALL CLARNV( 4, ISEED, N-J+1, AP( JC ) ) 663 IF( J.NE.IY ) THEN 664 AP( JC ) = CLARND( 5, ISEED )*TWO 665 ELSE 666 AP( JC ) = ZERO 667 END IF 668 JC = JC + N - J + 1 669 350 CONTINUE 670 END IF 671 CALL CLARNV( 2, ISEED, N, B ) 672 CALL CSSCAL( N, TWO, B, 1 ) 673* 674 ELSE IF( IMAT.EQ.17 ) THEN 675* 676* Type 17: Make the offdiagonal elements large to cause overflow 677* when adding a column of T. In the non-transposed case, the 678* matrix is constructed to cause overflow when adding a column in 679* every other step. 680* 681 TSCAL = UNFL / ULP 682 TSCAL = ( ONE-ULP ) / TSCAL 683 DO 360 J = 1, N*( N+1 ) / 2 684 AP( J ) = ZERO 685 360 CONTINUE 686 TEXP = ONE 687 IF( UPPER ) THEN 688 JC = ( N-1 )*N / 2 + 1 689 DO 370 J = N, 2, -2 690 AP( JC ) = -TSCAL / REAL( N+1 ) 691 AP( JC+J-1 ) = ONE 692 B( J ) = TEXP*( ONE-ULP ) 693 JC = JC - J + 1 694 AP( JC ) = -( TSCAL / REAL( N+1 ) ) / REAL( N+2 ) 695 AP( JC+J-2 ) = ONE 696 B( J-1 ) = TEXP*REAL( N*N+N-1 ) 697 TEXP = TEXP*TWO 698 JC = JC - J + 2 699 370 CONTINUE 700 B( 1 ) = ( REAL( N+1 ) / REAL( N+2 ) )*TSCAL 701 ELSE 702 JC = 1 703 DO 380 J = 1, N - 1, 2 704 AP( JC+N-J ) = -TSCAL / REAL( N+1 ) 705 AP( JC ) = ONE 706 B( J ) = TEXP*( ONE-ULP ) 707 JC = JC + N - J + 1 708 AP( JC+N-J-1 ) = -( TSCAL / REAL( N+1 ) ) / REAL( N+2 ) 709 AP( JC ) = ONE 710 B( J+1 ) = TEXP*REAL( N*N+N-1 ) 711 TEXP = TEXP*TWO 712 JC = JC + N - J 713 380 CONTINUE 714 B( N ) = ( REAL( N+1 ) / REAL( N+2 ) )*TSCAL 715 END IF 716* 717 ELSE IF( IMAT.EQ.18 ) THEN 718* 719* Type 18: Generate a unit triangular matrix with elements 720* between -1 and 1, and make the right hand side large so that it 721* requires scaling. 722* 723 IF( UPPER ) THEN 724 JC = 1 725 DO 390 J = 1, N 726 CALL CLARNV( 4, ISEED, J-1, AP( JC ) ) 727 AP( JC+J-1 ) = ZERO 728 JC = JC + J 729 390 CONTINUE 730 ELSE 731 JC = 1 732 DO 400 J = 1, N 733 IF( J.LT.N ) 734 $ CALL CLARNV( 4, ISEED, N-J, AP( JC+1 ) ) 735 AP( JC ) = ZERO 736 JC = JC + N - J + 1 737 400 CONTINUE 738 END IF 739* 740* Set the right hand side so that the largest value is BIGNUM. 741* 742 CALL CLARNV( 2, ISEED, N, B ) 743 IY = ICAMAX( N, B, 1 ) 744 BNORM = ABS( B( IY ) ) 745 BSCAL = BIGNUM / MAX( ONE, BNORM ) 746 CALL CSSCAL( N, BSCAL, B, 1 ) 747* 748 ELSE IF( IMAT.EQ.19 ) THEN 749* 750* Type 19: Generate a triangular matrix with elements between 751* BIGNUM/(n-1) and BIGNUM so that at least one of the column 752* norms will exceed BIGNUM. 753* 1/3/91: CLATPS no longer can handle this case 754* 755 TLEFT = BIGNUM / MAX( ONE, REAL( N-1 ) ) 756 TSCAL = BIGNUM*( REAL( N-1 ) / MAX( ONE, REAL( N ) ) ) 757 IF( UPPER ) THEN 758 JC = 1 759 DO 420 J = 1, N 760 CALL CLARNV( 5, ISEED, J, AP( JC ) ) 761 CALL SLARNV( 1, ISEED, J, RWORK ) 762 DO 410 I = 1, J 763 AP( JC+I-1 ) = AP( JC+I-1 )*( TLEFT+RWORK( I )*TSCAL ) 764 410 CONTINUE 765 JC = JC + J 766 420 CONTINUE 767 ELSE 768 JC = 1 769 DO 440 J = 1, N 770 CALL CLARNV( 5, ISEED, N-J+1, AP( JC ) ) 771 CALL SLARNV( 1, ISEED, N-J+1, RWORK ) 772 DO 430 I = J, N 773 AP( JC+I-J ) = AP( JC+I-J )* 774 $ ( TLEFT+RWORK( I-J+1 )*TSCAL ) 775 430 CONTINUE 776 JC = JC + N - J + 1 777 440 CONTINUE 778 END IF 779 CALL CLARNV( 2, ISEED, N, B ) 780 CALL CSSCAL( N, TWO, B, 1 ) 781 END IF 782* 783* Flip the matrix across its counter-diagonal if the transpose will 784* be used. 785* 786 IF( .NOT.LSAME( TRANS, 'N' ) ) THEN 787 IF( UPPER ) THEN 788 JJ = 1 789 JR = N*( N+1 ) / 2 790 DO 460 J = 1, N / 2 791 JL = JJ 792 DO 450 I = J, N - J 793 T = AP( JR-I+J ) 794 AP( JR-I+J ) = AP( JL ) 795 AP( JL ) = T 796 JL = JL + I 797 450 CONTINUE 798 JJ = JJ + J + 1 799 JR = JR - ( N-J+1 ) 800 460 CONTINUE 801 ELSE 802 JL = 1 803 JJ = N*( N+1 ) / 2 804 DO 480 J = 1, N / 2 805 JR = JJ 806 DO 470 I = J, N - J 807 T = AP( JL+I-J ) 808 AP( JL+I-J ) = AP( JR ) 809 AP( JR ) = T 810 JR = JR - I 811 470 CONTINUE 812 JL = JL + N - J + 1 813 JJ = JJ - J - 1 814 480 CONTINUE 815 END IF 816 END IF 817* 818 RETURN 819* 820* End of CLATTP 821* 822 END 823