1*> \brief \b DLATMT 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 DLATMT( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, 12* RANK, KL, KU, PACK, A, LDA, WORK, INFO ) 13* 14* .. Scalar Arguments .. 15* DOUBLE PRECISION COND, DMAX 16* INTEGER INFO, KL, KU, LDA, M, MODE, N, RANK 17* CHARACTER DIST, PACK, SYM 18* .. 19* .. Array Arguments .. 20* DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * ) 21* INTEGER ISEED( 4 ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> DLATMT generates random matrices with specified singular values 31*> (or symmetric/hermitian with specified eigenvalues) 32*> for testing LAPACK programs. 33*> 34*> DLATMT operates by applying the following sequence of 35*> operations: 36*> 37*> Set the diagonal to D, where D may be input or 38*> computed according to MODE, COND, DMAX, and SYM 39*> as described below. 40*> 41*> Generate a matrix with the appropriate band structure, by one 42*> of two methods: 43*> 44*> Method A: 45*> Generate a dense M x N matrix by multiplying D on the left 46*> and the right by random unitary matrices, then: 47*> 48*> Reduce the bandwidth according to KL and KU, using 49*> Householder transformations. 50*> 51*> Method B: 52*> Convert the bandwidth-0 (i.e., diagonal) matrix to a 53*> bandwidth-1 matrix using Givens rotations, "chasing" 54*> out-of-band elements back, much as in QR; then 55*> convert the bandwidth-1 to a bandwidth-2 matrix, etc. 56*> Note that for reasonably small bandwidths (relative to 57*> M and N) this requires less storage, as a dense matrix 58*> is not generated. Also, for symmetric matrices, only 59*> one triangle is generated. 60*> 61*> Method A is chosen if the bandwidth is a large fraction of the 62*> order of the matrix, and LDA is at least M (so a dense 63*> matrix can be stored.) Method B is chosen if the bandwidth 64*> is small (< 1/2 N for symmetric, < .3 N+M for 65*> non-symmetric), or LDA is less than M and not less than the 66*> bandwidth. 67*> 68*> Pack the matrix if desired. Options specified by PACK are: 69*> no packing 70*> zero out upper half (if symmetric) 71*> zero out lower half (if symmetric) 72*> store the upper half columnwise (if symmetric or upper 73*> triangular) 74*> store the lower half columnwise (if symmetric or lower 75*> triangular) 76*> store the lower triangle in banded format (if symmetric 77*> or lower triangular) 78*> store the upper triangle in banded format (if symmetric 79*> or upper triangular) 80*> store the entire matrix in banded format 81*> If Method B is chosen, and band format is specified, then the 82*> matrix will be generated in the band format, so no repacking 83*> will be necessary. 84*> \endverbatim 85* 86* Arguments: 87* ========== 88* 89*> \param[in] M 90*> \verbatim 91*> M is INTEGER 92*> The number of rows of A. Not modified. 93*> \endverbatim 94*> 95*> \param[in] N 96*> \verbatim 97*> N is INTEGER 98*> The number of columns of A. Not modified. 99*> \endverbatim 100*> 101*> \param[in] DIST 102*> \verbatim 103*> DIST is CHARACTER*1 104*> On entry, DIST specifies the type of distribution to be used 105*> to generate the random eigen-/singular values. 106*> 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform ) 107*> 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric ) 108*> 'N' => NORMAL( 0, 1 ) ( 'N' for normal ) 109*> Not modified. 110*> \endverbatim 111*> 112*> \param[in,out] ISEED 113*> \verbatim 114*> ISEED is INTEGER array, dimension ( 4 ) 115*> On entry ISEED specifies the seed of the random number 116*> generator. They should lie between 0 and 4095 inclusive, 117*> and ISEED(4) should be odd. The random number generator 118*> uses a linear congruential sequence limited to small 119*> integers, and so should produce machine independent 120*> random numbers. The values of ISEED are changed on 121*> exit, and can be used in the next call to DLATMT 122*> to continue the same random number sequence. 123*> Changed on exit. 124*> \endverbatim 125*> 126*> \param[in] SYM 127*> \verbatim 128*> SYM is CHARACTER*1 129*> If SYM='S' or 'H', the generated matrix is symmetric, with 130*> eigenvalues specified by D, COND, MODE, and DMAX; they 131*> may be positive, negative, or zero. 132*> If SYM='P', the generated matrix is symmetric, with 133*> eigenvalues (= singular values) specified by D, COND, 134*> MODE, and DMAX; they will not be negative. 135*> If SYM='N', the generated matrix is nonsymmetric, with 136*> singular values specified by D, COND, MODE, and DMAX; 137*> they will not be negative. 138*> Not modified. 139*> \endverbatim 140*> 141*> \param[in,out] D 142*> \verbatim 143*> D is DOUBLE PRECISION array, dimension ( MIN( M , N ) ) 144*> This array is used to specify the singular values or 145*> eigenvalues of A (see SYM, above.) If MODE=0, then D is 146*> assumed to contain the singular/eigenvalues, otherwise 147*> they will be computed according to MODE, COND, and DMAX, 148*> and placed in D. 149*> Modified if MODE is nonzero. 150*> \endverbatim 151*> 152*> \param[in] MODE 153*> \verbatim 154*> MODE is INTEGER 155*> On entry this describes how the singular/eigenvalues are to 156*> be specified: 157*> MODE = 0 means use D as input 158*> 159*> MODE = 1 sets D(1)=1 and D(2:RANK)=1.0/COND 160*> MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND 161*> MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1)) 162*> 163*> MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND) 164*> MODE = 5 sets D to random numbers in the range 165*> ( 1/COND , 1 ) such that their logarithms 166*> are uniformly distributed. 167*> MODE = 6 set D to random numbers from same distribution 168*> as the rest of the matrix. 169*> MODE < 0 has the same meaning as ABS(MODE), except that 170*> the order of the elements of D is reversed. 171*> Thus if MODE is positive, D has entries ranging from 172*> 1 to 1/COND, if negative, from 1/COND to 1, 173*> If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then 174*> the elements of D will also be multiplied by a random 175*> sign (i.e., +1 or -1.) 176*> Not modified. 177*> \endverbatim 178*> 179*> \param[in] COND 180*> \verbatim 181*> COND is DOUBLE PRECISION 182*> On entry, this is used as described under MODE above. 183*> If used, it must be >= 1. Not modified. 184*> \endverbatim 185*> 186*> \param[in] DMAX 187*> \verbatim 188*> DMAX is DOUBLE PRECISION 189*> If MODE is neither -6, 0 nor 6, the contents of D, as 190*> computed according to MODE and COND, will be scaled by 191*> DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or 192*> singular value (which is to say the norm) will be abs(DMAX). 193*> Note that DMAX need not be positive: if DMAX is negative 194*> (or zero), D will be scaled by a negative number (or zero). 195*> Not modified. 196*> \endverbatim 197*> 198*> \param[in] RANK 199*> \verbatim 200*> RANK is INTEGER 201*> The rank of matrix to be generated for modes 1,2,3 only. 202*> D( RANK+1:N ) = 0. 203*> Not modified. 204*> \endverbatim 205*> 206*> \param[in] KL 207*> \verbatim 208*> KL is INTEGER 209*> This specifies the lower bandwidth of the matrix. For 210*> example, KL=0 implies upper triangular, KL=1 implies upper 211*> Hessenberg, and KL being at least M-1 means that the matrix 212*> has full lower bandwidth. KL must equal KU if the matrix 213*> is symmetric. 214*> Not modified. 215*> \endverbatim 216*> 217*> \param[in] KU 218*> \verbatim 219*> KU is INTEGER 220*> This specifies the upper bandwidth of the matrix. For 221*> example, KU=0 implies lower triangular, KU=1 implies lower 222*> Hessenberg, and KU being at least N-1 means that the matrix 223*> has full upper bandwidth. KL must equal KU if the matrix 224*> is symmetric. 225*> Not modified. 226*> \endverbatim 227*> 228*> \param[in] PACK 229*> \verbatim 230*> PACK is CHARACTER*1 231*> This specifies packing of matrix as follows: 232*> 'N' => no packing 233*> 'U' => zero out all subdiagonal entries (if symmetric) 234*> 'L' => zero out all superdiagonal entries (if symmetric) 235*> 'C' => store the upper triangle columnwise 236*> (only if the matrix is symmetric or upper triangular) 237*> 'R' => store the lower triangle columnwise 238*> (only if the matrix is symmetric or lower triangular) 239*> 'B' => store the lower triangle in band storage scheme 240*> (only if matrix symmetric or lower triangular) 241*> 'Q' => store the upper triangle in band storage scheme 242*> (only if matrix symmetric or upper triangular) 243*> 'Z' => store the entire matrix in band storage scheme 244*> (pivoting can be provided for by using this 245*> option to store A in the trailing rows of 246*> the allocated storage) 247*> 248*> Using these options, the various LAPACK packed and banded 249*> storage schemes can be obtained: 250*> GB - use 'Z' 251*> PB, SB or TB - use 'B' or 'Q' 252*> PP, SP or TP - use 'C' or 'R' 253*> 254*> If two calls to DLATMT differ only in the PACK parameter, 255*> they will generate mathematically equivalent matrices. 256*> Not modified. 257*> \endverbatim 258*> 259*> \param[in,out] A 260*> \verbatim 261*> A is DOUBLE PRECISION array, dimension ( LDA, N ) 262*> On exit A is the desired test matrix. A is first generated 263*> in full (unpacked) form, and then packed, if so specified 264*> by PACK. Thus, the first M elements of the first N 265*> columns will always be modified. If PACK specifies a 266*> packed or banded storage scheme, all LDA elements of the 267*> first N columns will be modified; the elements of the 268*> array which do not correspond to elements of the generated 269*> matrix are set to zero. 270*> Modified. 271*> \endverbatim 272*> 273*> \param[in] LDA 274*> \verbatim 275*> LDA is INTEGER 276*> LDA specifies the first dimension of A as declared in the 277*> calling program. If PACK='N', 'U', 'L', 'C', or 'R', then 278*> LDA must be at least M. If PACK='B' or 'Q', then LDA must 279*> be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)). 280*> If PACK='Z', LDA must be large enough to hold the packed 281*> array: MIN( KU, N-1) + MIN( KL, M-1) + 1. 282*> Not modified. 283*> \endverbatim 284*> 285*> \param[out] WORK 286*> \verbatim 287*> WORK is DOUBLE PRECISION array, dimension ( 3*MAX( N , M ) ) 288*> Workspace. 289*> Modified. 290*> \endverbatim 291*> 292*> \param[out] INFO 293*> \verbatim 294*> INFO is INTEGER 295*> Error code. On exit, INFO will be set to one of the 296*> following values: 297*> 0 => normal return 298*> -1 => M negative or unequal to N and SYM='S', 'H', or 'P' 299*> -2 => N negative 300*> -3 => DIST illegal string 301*> -5 => SYM illegal string 302*> -7 => MODE not in range -6 to 6 303*> -8 => COND less than 1.0, and MODE neither -6, 0 nor 6 304*> -10 => KL negative 305*> -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL 306*> -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N'; 307*> or PACK='C' or 'Q' and SYM='N' and KL is not zero; 308*> or PACK='R' or 'B' and SYM='N' and KU is not zero; 309*> or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not 310*> N. 311*> -14 => LDA is less than M, or PACK='Z' and LDA is less than 312*> MIN(KU,N-1) + MIN(KL,M-1) + 1. 313*> 1 => Error return from DLATM7 314*> 2 => Cannot scale to DMAX (max. sing. value is 0) 315*> 3 => Error return from DLAGGE or DLAGSY 316*> \endverbatim 317* 318* Authors: 319* ======== 320* 321*> \author Univ. of Tennessee 322*> \author Univ. of California Berkeley 323*> \author Univ. of Colorado Denver 324*> \author NAG Ltd. 325* 326*> \ingroup double_matgen 327* 328* ===================================================================== 329 SUBROUTINE DLATMT( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, 330 $ RANK, KL, KU, PACK, A, LDA, WORK, INFO ) 331* 332* -- LAPACK computational routine -- 333* -- LAPACK is a software package provided by Univ. of Tennessee, -- 334* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 335* 336* .. Scalar Arguments .. 337 DOUBLE PRECISION COND, DMAX 338 INTEGER INFO, KL, KU, LDA, M, MODE, N, RANK 339 CHARACTER DIST, PACK, SYM 340* .. 341* .. Array Arguments .. 342 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * ) 343 INTEGER ISEED( 4 ) 344* .. 345* 346* ===================================================================== 347* 348* .. Parameters .. 349 DOUBLE PRECISION ZERO 350 PARAMETER ( ZERO = 0.0D0 ) 351 DOUBLE PRECISION ONE 352 PARAMETER ( ONE = 1.0D0 ) 353 DOUBLE PRECISION TWOPI 354 PARAMETER ( TWOPI = 6.28318530717958647692528676655900576839D+0 ) 355* .. 356* .. Local Scalars .. 357 DOUBLE PRECISION ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP 358 INTEGER I, IC, ICOL, IDIST, IENDCH, IINFO, IL, ILDA, 359 $ IOFFG, IOFFST, IPACK, IPACKG, IR, IR1, IR2, 360 $ IROW, IRSIGN, ISKEW, ISYM, ISYMPK, J, JC, JCH, 361 $ JKL, JKU, JR, K, LLB, MINLDA, MNMIN, MR, NC, 362 $ UUB 363 LOGICAL GIVENS, ILEXTR, ILTEMP, TOPDWN 364* .. 365* .. External Functions .. 366 DOUBLE PRECISION DLARND 367 LOGICAL LSAME 368 EXTERNAL DLARND, LSAME 369* .. 370* .. External Subroutines .. 371 EXTERNAL DLATM7, DCOPY, DLAGGE, DLAGSY, DLAROT, 372 $ DLARTG, DLASET, DSCAL, XERBLA 373* .. 374* .. Intrinsic Functions .. 375 INTRINSIC ABS, COS, DBLE, MAX, MIN, MOD, SIN 376* .. 377* .. Executable Statements .. 378* 379* 1) Decode and Test the input parameters. 380* Initialize flags & seed. 381* 382 INFO = 0 383* 384* Quick return if possible 385* 386 IF( M.EQ.0 .OR. N.EQ.0 ) 387 $ RETURN 388* 389* Decode DIST 390* 391 IF( LSAME( DIST, 'U' ) ) THEN 392 IDIST = 1 393 ELSE IF( LSAME( DIST, 'S' ) ) THEN 394 IDIST = 2 395 ELSE IF( LSAME( DIST, 'N' ) ) THEN 396 IDIST = 3 397 ELSE 398 IDIST = -1 399 END IF 400* 401* Decode SYM 402* 403 IF( LSAME( SYM, 'N' ) ) THEN 404 ISYM = 1 405 IRSIGN = 0 406 ELSE IF( LSAME( SYM, 'P' ) ) THEN 407 ISYM = 2 408 IRSIGN = 0 409 ELSE IF( LSAME( SYM, 'S' ) ) THEN 410 ISYM = 2 411 IRSIGN = 1 412 ELSE IF( LSAME( SYM, 'H' ) ) THEN 413 ISYM = 2 414 IRSIGN = 1 415 ELSE 416 ISYM = -1 417 END IF 418* 419* Decode PACK 420* 421 ISYMPK = 0 422 IF( LSAME( PACK, 'N' ) ) THEN 423 IPACK = 0 424 ELSE IF( LSAME( PACK, 'U' ) ) THEN 425 IPACK = 1 426 ISYMPK = 1 427 ELSE IF( LSAME( PACK, 'L' ) ) THEN 428 IPACK = 2 429 ISYMPK = 1 430 ELSE IF( LSAME( PACK, 'C' ) ) THEN 431 IPACK = 3 432 ISYMPK = 2 433 ELSE IF( LSAME( PACK, 'R' ) ) THEN 434 IPACK = 4 435 ISYMPK = 3 436 ELSE IF( LSAME( PACK, 'B' ) ) THEN 437 IPACK = 5 438 ISYMPK = 3 439 ELSE IF( LSAME( PACK, 'Q' ) ) THEN 440 IPACK = 6 441 ISYMPK = 2 442 ELSE IF( LSAME( PACK, 'Z' ) ) THEN 443 IPACK = 7 444 ELSE 445 IPACK = -1 446 END IF 447* 448* Set certain internal parameters 449* 450 MNMIN = MIN( M, N ) 451 LLB = MIN( KL, M-1 ) 452 UUB = MIN( KU, N-1 ) 453 MR = MIN( M, N+LLB ) 454 NC = MIN( N, M+UUB ) 455* 456 IF( IPACK.EQ.5 .OR. IPACK.EQ.6 ) THEN 457 MINLDA = UUB + 1 458 ELSE IF( IPACK.EQ.7 ) THEN 459 MINLDA = LLB + UUB + 1 460 ELSE 461 MINLDA = M 462 END IF 463* 464* Use Givens rotation method if bandwidth small enough, 465* or if LDA is too small to store the matrix unpacked. 466* 467 GIVENS = .FALSE. 468 IF( ISYM.EQ.1 ) THEN 469 IF( DBLE( LLB+UUB ).LT.0.3D0*DBLE( MAX( 1, MR+NC ) ) ) 470 $ GIVENS = .TRUE. 471 ELSE 472 IF( 2*LLB.LT.M ) 473 $ GIVENS = .TRUE. 474 END IF 475 IF( LDA.LT.M .AND. LDA.GE.MINLDA ) 476 $ GIVENS = .TRUE. 477* 478* Set INFO if an error 479* 480 IF( M.LT.0 ) THEN 481 INFO = -1 482 ELSE IF( M.NE.N .AND. ISYM.NE.1 ) THEN 483 INFO = -1 484 ELSE IF( N.LT.0 ) THEN 485 INFO = -2 486 ELSE IF( IDIST.EQ.-1 ) THEN 487 INFO = -3 488 ELSE IF( ISYM.EQ.-1 ) THEN 489 INFO = -5 490 ELSE IF( ABS( MODE ).GT.6 ) THEN 491 INFO = -7 492 ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE ) 493 $ THEN 494 INFO = -8 495 ELSE IF( KL.LT.0 ) THEN 496 INFO = -10 497 ELSE IF( KU.LT.0 .OR. ( ISYM.NE.1 .AND. KL.NE.KU ) ) THEN 498 INFO = -11 499 ELSE IF( IPACK.EQ.-1 .OR. ( ISYMPK.EQ.1 .AND. ISYM.EQ.1 ) .OR. 500 $ ( ISYMPK.EQ.2 .AND. ISYM.EQ.1 .AND. KL.GT.0 ) .OR. 501 $ ( ISYMPK.EQ.3 .AND. ISYM.EQ.1 .AND. KU.GT.0 ) .OR. 502 $ ( ISYMPK.NE.0 .AND. M.NE.N ) ) THEN 503 INFO = -12 504 ELSE IF( LDA.LT.MAX( 1, MINLDA ) ) THEN 505 INFO = -14 506 END IF 507* 508 IF( INFO.NE.0 ) THEN 509 CALL XERBLA( 'DLATMT', -INFO ) 510 RETURN 511 END IF 512* 513* Initialize random number generator 514* 515 DO 100 I = 1, 4 516 ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 ) 517 100 CONTINUE 518* 519 IF( MOD( ISEED( 4 ), 2 ).NE.1 ) 520 $ ISEED( 4 ) = ISEED( 4 ) + 1 521* 522* 2) Set up D if indicated. 523* 524* Compute D according to COND and MODE 525* 526 CALL DLATM7( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, RANK, 527 $ IINFO ) 528 IF( IINFO.NE.0 ) THEN 529 INFO = 1 530 RETURN 531 END IF 532* 533* Choose Top-Down if D is (apparently) increasing, 534* Bottom-Up if D is (apparently) decreasing. 535* 536 IF( ABS( D( 1 ) ).LE.ABS( D( RANK ) ) ) THEN 537 TOPDWN = .TRUE. 538 ELSE 539 TOPDWN = .FALSE. 540 END IF 541* 542 IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN 543* 544* Scale by DMAX 545* 546 TEMP = ABS( D( 1 ) ) 547 DO 110 I = 2, RANK 548 TEMP = MAX( TEMP, ABS( D( I ) ) ) 549 110 CONTINUE 550* 551 IF( TEMP.GT.ZERO ) THEN 552 ALPHA = DMAX / TEMP 553 ELSE 554 INFO = 2 555 RETURN 556 END IF 557* 558 CALL DSCAL( RANK, ALPHA, D, 1 ) 559* 560 END IF 561* 562* 3) Generate Banded Matrix using Givens rotations. 563* Also the special case of UUB=LLB=0 564* 565* Compute Addressing constants to cover all 566* storage formats. Whether GE, SY, GB, or SB, 567* upper or lower triangle or both, 568* the (i,j)-th element is in 569* A( i - ISKEW*j + IOFFST, j ) 570* 571 IF( IPACK.GT.4 ) THEN 572 ILDA = LDA - 1 573 ISKEW = 1 574 IF( IPACK.GT.5 ) THEN 575 IOFFST = UUB + 1 576 ELSE 577 IOFFST = 1 578 END IF 579 ELSE 580 ILDA = LDA 581 ISKEW = 0 582 IOFFST = 0 583 END IF 584* 585* IPACKG is the format that the matrix is generated in. If this is 586* different from IPACK, then the matrix must be repacked at the 587* end. It also signals how to compute the norm, for scaling. 588* 589 IPACKG = 0 590 CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA ) 591* 592* Diagonal Matrix -- We are done, unless it 593* is to be stored SP/PP/TP (PACK='R' or 'C') 594* 595 IF( LLB.EQ.0 .AND. UUB.EQ.0 ) THEN 596 CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 ) 597 IF( IPACK.LE.2 .OR. IPACK.GE.5 ) 598 $ IPACKG = IPACK 599* 600 ELSE IF( GIVENS ) THEN 601* 602* Check whether to use Givens rotations, 603* Householder transformations, or nothing. 604* 605 IF( ISYM.EQ.1 ) THEN 606* 607* Non-symmetric -- A = U D V 608* 609 IF( IPACK.GT.4 ) THEN 610 IPACKG = IPACK 611 ELSE 612 IPACKG = 0 613 END IF 614* 615 CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 ) 616* 617 IF( TOPDWN ) THEN 618 JKL = 0 619 DO 140 JKU = 1, UUB 620* 621* Transform from bandwidth JKL, JKU-1 to JKL, JKU 622* 623* Last row actually rotated is M 624* Last column actually rotated is MIN( M+JKU, N ) 625* 626 DO 130 JR = 1, MIN( M+JKU, N ) + JKL - 1 627 EXTRA = ZERO 628 ANGLE = TWOPI*DLARND( 1, ISEED ) 629 C = COS( ANGLE ) 630 S = SIN( ANGLE ) 631 ICOL = MAX( 1, JR-JKL ) 632 IF( JR.LT.M ) THEN 633 IL = MIN( N, JR+JKU ) + 1 - ICOL 634 CALL DLAROT( .TRUE., JR.GT.JKL, .FALSE., IL, C, 635 $ S, A( JR-ISKEW*ICOL+IOFFST, ICOL ), 636 $ ILDA, EXTRA, DUMMY ) 637 END IF 638* 639* Chase "EXTRA" back up 640* 641 IR = JR 642 IC = ICOL 643 DO 120 JCH = JR - JKL, 1, -JKL - JKU 644 IF( IR.LT.M ) THEN 645 CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST, 646 $ IC+1 ), EXTRA, C, S, DUMMY ) 647 END IF 648 IROW = MAX( 1, JCH-JKU ) 649 IL = IR + 2 - IROW 650 TEMP = ZERO 651 ILTEMP = JCH.GT.JKU 652 CALL DLAROT( .FALSE., ILTEMP, .TRUE., IL, C, -S, 653 $ A( IROW-ISKEW*IC+IOFFST, IC ), 654 $ ILDA, TEMP, EXTRA ) 655 IF( ILTEMP ) THEN 656 CALL DLARTG( A( IROW+1-ISKEW*( IC+1 )+IOFFST, 657 $ IC+1 ), TEMP, C, S, DUMMY ) 658 ICOL = MAX( 1, JCH-JKU-JKL ) 659 IL = IC + 2 - ICOL 660 EXTRA = ZERO 661 CALL DLAROT( .TRUE., JCH.GT.JKU+JKL, .TRUE., 662 $ IL, C, -S, A( IROW-ISKEW*ICOL+ 663 $ IOFFST, ICOL ), ILDA, EXTRA, 664 $ TEMP ) 665 IC = ICOL 666 IR = IROW 667 END IF 668 120 CONTINUE 669 130 CONTINUE 670 140 CONTINUE 671* 672 JKU = UUB 673 DO 170 JKL = 1, LLB 674* 675* Transform from bandwidth JKL-1, JKU to JKL, JKU 676* 677 DO 160 JC = 1, MIN( N+JKL, M ) + JKU - 1 678 EXTRA = ZERO 679 ANGLE = TWOPI*DLARND( 1, ISEED ) 680 C = COS( ANGLE ) 681 S = SIN( ANGLE ) 682 IROW = MAX( 1, JC-JKU ) 683 IF( JC.LT.N ) THEN 684 IL = MIN( M, JC+JKL ) + 1 - IROW 685 CALL DLAROT( .FALSE., JC.GT.JKU, .FALSE., IL, C, 686 $ S, A( IROW-ISKEW*JC+IOFFST, JC ), 687 $ ILDA, EXTRA, DUMMY ) 688 END IF 689* 690* Chase "EXTRA" back up 691* 692 IC = JC 693 IR = IROW 694 DO 150 JCH = JC - JKU, 1, -JKL - JKU 695 IF( IC.LT.N ) THEN 696 CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST, 697 $ IC+1 ), EXTRA, C, S, DUMMY ) 698 END IF 699 ICOL = MAX( 1, JCH-JKL ) 700 IL = IC + 2 - ICOL 701 TEMP = ZERO 702 ILTEMP = JCH.GT.JKL 703 CALL DLAROT( .TRUE., ILTEMP, .TRUE., IL, C, -S, 704 $ A( IR-ISKEW*ICOL+IOFFST, ICOL ), 705 $ ILDA, TEMP, EXTRA ) 706 IF( ILTEMP ) THEN 707 CALL DLARTG( A( IR+1-ISKEW*( ICOL+1 )+IOFFST, 708 $ ICOL+1 ), TEMP, C, S, DUMMY ) 709 IROW = MAX( 1, JCH-JKL-JKU ) 710 IL = IR + 2 - IROW 711 EXTRA = ZERO 712 CALL DLAROT( .FALSE., JCH.GT.JKL+JKU, .TRUE., 713 $ IL, C, -S, A( IROW-ISKEW*ICOL+ 714 $ IOFFST, ICOL ), ILDA, EXTRA, 715 $ TEMP ) 716 IC = ICOL 717 IR = IROW 718 END IF 719 150 CONTINUE 720 160 CONTINUE 721 170 CONTINUE 722* 723 ELSE 724* 725* Bottom-Up -- Start at the bottom right. 726* 727 JKL = 0 728 DO 200 JKU = 1, UUB 729* 730* Transform from bandwidth JKL, JKU-1 to JKL, JKU 731* 732* First row actually rotated is M 733* First column actually rotated is MIN( M+JKU, N ) 734* 735 IENDCH = MIN( M, N+JKL ) - 1 736 DO 190 JC = MIN( M+JKU, N ) - 1, 1 - JKL, -1 737 EXTRA = ZERO 738 ANGLE = TWOPI*DLARND( 1, ISEED ) 739 C = COS( ANGLE ) 740 S = SIN( ANGLE ) 741 IROW = MAX( 1, JC-JKU+1 ) 742 IF( JC.GT.0 ) THEN 743 IL = MIN( M, JC+JKL+1 ) + 1 - IROW 744 CALL DLAROT( .FALSE., .FALSE., JC+JKL.LT.M, IL, 745 $ C, S, A( IROW-ISKEW*JC+IOFFST, 746 $ JC ), ILDA, DUMMY, EXTRA ) 747 END IF 748* 749* Chase "EXTRA" back down 750* 751 IC = JC 752 DO 180 JCH = JC + JKL, IENDCH, JKL + JKU 753 ILEXTR = IC.GT.0 754 IF( ILEXTR ) THEN 755 CALL DLARTG( A( JCH-ISKEW*IC+IOFFST, IC ), 756 $ EXTRA, C, S, DUMMY ) 757 END IF 758 IC = MAX( 1, IC ) 759 ICOL = MIN( N-1, JCH+JKU ) 760 ILTEMP = JCH + JKU.LT.N 761 TEMP = ZERO 762 CALL DLAROT( .TRUE., ILEXTR, ILTEMP, ICOL+2-IC, 763 $ C, S, A( JCH-ISKEW*IC+IOFFST, IC ), 764 $ ILDA, EXTRA, TEMP ) 765 IF( ILTEMP ) THEN 766 CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFST, 767 $ ICOL ), TEMP, C, S, DUMMY ) 768 IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH 769 EXTRA = ZERO 770 CALL DLAROT( .FALSE., .TRUE., 771 $ JCH+JKL+JKU.LE.IENDCH, IL, C, S, 772 $ A( JCH-ISKEW*ICOL+IOFFST, 773 $ ICOL ), ILDA, TEMP, EXTRA ) 774 IC = ICOL 775 END IF 776 180 CONTINUE 777 190 CONTINUE 778 200 CONTINUE 779* 780 JKU = UUB 781 DO 230 JKL = 1, LLB 782* 783* Transform from bandwidth JKL-1, JKU to JKL, JKU 784* 785* First row actually rotated is MIN( N+JKL, M ) 786* First column actually rotated is N 787* 788 IENDCH = MIN( N, M+JKU ) - 1 789 DO 220 JR = MIN( N+JKL, M ) - 1, 1 - JKU, -1 790 EXTRA = ZERO 791 ANGLE = TWOPI*DLARND( 1, ISEED ) 792 C = COS( ANGLE ) 793 S = SIN( ANGLE ) 794 ICOL = MAX( 1, JR-JKL+1 ) 795 IF( JR.GT.0 ) THEN 796 IL = MIN( N, JR+JKU+1 ) + 1 - ICOL 797 CALL DLAROT( .TRUE., .FALSE., JR+JKU.LT.N, IL, 798 $ C, S, A( JR-ISKEW*ICOL+IOFFST, 799 $ ICOL ), ILDA, DUMMY, EXTRA ) 800 END IF 801* 802* Chase "EXTRA" back down 803* 804 IR = JR 805 DO 210 JCH = JR + JKU, IENDCH, JKL + JKU 806 ILEXTR = IR.GT.0 807 IF( ILEXTR ) THEN 808 CALL DLARTG( A( IR-ISKEW*JCH+IOFFST, JCH ), 809 $ EXTRA, C, S, DUMMY ) 810 END IF 811 IR = MAX( 1, IR ) 812 IROW = MIN( M-1, JCH+JKL ) 813 ILTEMP = JCH + JKL.LT.M 814 TEMP = ZERO 815 CALL DLAROT( .FALSE., ILEXTR, ILTEMP, IROW+2-IR, 816 $ C, S, A( IR-ISKEW*JCH+IOFFST, 817 $ JCH ), ILDA, EXTRA, TEMP ) 818 IF( ILTEMP ) THEN 819 CALL DLARTG( A( IROW-ISKEW*JCH+IOFFST, JCH ), 820 $ TEMP, C, S, DUMMY ) 821 IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH 822 EXTRA = ZERO 823 CALL DLAROT( .TRUE., .TRUE., 824 $ JCH+JKL+JKU.LE.IENDCH, IL, C, S, 825 $ A( IROW-ISKEW*JCH+IOFFST, JCH ), 826 $ ILDA, TEMP, EXTRA ) 827 IR = IROW 828 END IF 829 210 CONTINUE 830 220 CONTINUE 831 230 CONTINUE 832 END IF 833* 834 ELSE 835* 836* Symmetric -- A = U D U' 837* 838 IPACKG = IPACK 839 IOFFG = IOFFST 840* 841 IF( TOPDWN ) THEN 842* 843* Top-Down -- Generate Upper triangle only 844* 845 IF( IPACK.GE.5 ) THEN 846 IPACKG = 6 847 IOFFG = UUB + 1 848 ELSE 849 IPACKG = 1 850 END IF 851 CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 ) 852* 853 DO 260 K = 1, UUB 854 DO 250 JC = 1, N - 1 855 IROW = MAX( 1, JC-K ) 856 IL = MIN( JC+1, K+2 ) 857 EXTRA = ZERO 858 TEMP = A( JC-ISKEW*( JC+1 )+IOFFG, JC+1 ) 859 ANGLE = TWOPI*DLARND( 1, ISEED ) 860 C = COS( ANGLE ) 861 S = SIN( ANGLE ) 862 CALL DLAROT( .FALSE., JC.GT.K, .TRUE., IL, C, S, 863 $ A( IROW-ISKEW*JC+IOFFG, JC ), ILDA, 864 $ EXTRA, TEMP ) 865 CALL DLAROT( .TRUE., .TRUE., .FALSE., 866 $ MIN( K, N-JC )+1, C, S, 867 $ A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA, 868 $ TEMP, DUMMY ) 869* 870* Chase EXTRA back up the matrix 871* 872 ICOL = JC 873 DO 240 JCH = JC - K, 1, -K 874 CALL DLARTG( A( JCH+1-ISKEW*( ICOL+1 )+IOFFG, 875 $ ICOL+1 ), EXTRA, C, S, DUMMY ) 876 TEMP = A( JCH-ISKEW*( JCH+1 )+IOFFG, JCH+1 ) 877 CALL DLAROT( .TRUE., .TRUE., .TRUE., K+2, C, -S, 878 $ A( ( 1-ISKEW )*JCH+IOFFG, JCH ), 879 $ ILDA, TEMP, EXTRA ) 880 IROW = MAX( 1, JCH-K ) 881 IL = MIN( JCH+1, K+2 ) 882 EXTRA = ZERO 883 CALL DLAROT( .FALSE., JCH.GT.K, .TRUE., IL, C, 884 $ -S, A( IROW-ISKEW*JCH+IOFFG, JCH ), 885 $ ILDA, EXTRA, TEMP ) 886 ICOL = JCH 887 240 CONTINUE 888 250 CONTINUE 889 260 CONTINUE 890* 891* If we need lower triangle, copy from upper. Note that 892* the order of copying is chosen to work for 'q' -> 'b' 893* 894 IF( IPACK.NE.IPACKG .AND. IPACK.NE.3 ) THEN 895 DO 280 JC = 1, N 896 IROW = IOFFST - ISKEW*JC 897 DO 270 JR = JC, MIN( N, JC+UUB ) 898 A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR ) 899 270 CONTINUE 900 280 CONTINUE 901 IF( IPACK.EQ.5 ) THEN 902 DO 300 JC = N - UUB + 1, N 903 DO 290 JR = N + 2 - JC, UUB + 1 904 A( JR, JC ) = ZERO 905 290 CONTINUE 906 300 CONTINUE 907 END IF 908 IF( IPACKG.EQ.6 ) THEN 909 IPACKG = IPACK 910 ELSE 911 IPACKG = 0 912 END IF 913 END IF 914 ELSE 915* 916* Bottom-Up -- Generate Lower triangle only 917* 918 IF( IPACK.GE.5 ) THEN 919 IPACKG = 5 920 IF( IPACK.EQ.6 ) 921 $ IOFFG = 1 922 ELSE 923 IPACKG = 2 924 END IF 925 CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 ) 926* 927 DO 330 K = 1, UUB 928 DO 320 JC = N - 1, 1, -1 929 IL = MIN( N+1-JC, K+2 ) 930 EXTRA = ZERO 931 TEMP = A( 1+( 1-ISKEW )*JC+IOFFG, JC ) 932 ANGLE = TWOPI*DLARND( 1, ISEED ) 933 C = COS( ANGLE ) 934 S = -SIN( ANGLE ) 935 CALL DLAROT( .FALSE., .TRUE., N-JC.GT.K, IL, C, S, 936 $ A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA, 937 $ TEMP, EXTRA ) 938 ICOL = MAX( 1, JC-K+1 ) 939 CALL DLAROT( .TRUE., .FALSE., .TRUE., JC+2-ICOL, C, 940 $ S, A( JC-ISKEW*ICOL+IOFFG, ICOL ), 941 $ ILDA, DUMMY, TEMP ) 942* 943* Chase EXTRA back down the matrix 944* 945 ICOL = JC 946 DO 310 JCH = JC + K, N - 1, K 947 CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFG, ICOL ), 948 $ EXTRA, C, S, DUMMY ) 949 TEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH ) 950 CALL DLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S, 951 $ A( JCH-ISKEW*ICOL+IOFFG, ICOL ), 952 $ ILDA, EXTRA, TEMP ) 953 IL = MIN( N+1-JCH, K+2 ) 954 EXTRA = ZERO 955 CALL DLAROT( .FALSE., .TRUE., N-JCH.GT.K, IL, C, 956 $ S, A( ( 1-ISKEW )*JCH+IOFFG, JCH ), 957 $ ILDA, TEMP, EXTRA ) 958 ICOL = JCH 959 310 CONTINUE 960 320 CONTINUE 961 330 CONTINUE 962* 963* If we need upper triangle, copy from lower. Note that 964* the order of copying is chosen to work for 'b' -> 'q' 965* 966 IF( IPACK.NE.IPACKG .AND. IPACK.NE.4 ) THEN 967 DO 350 JC = N, 1, -1 968 IROW = IOFFST - ISKEW*JC 969 DO 340 JR = JC, MAX( 1, JC-UUB ), -1 970 A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR ) 971 340 CONTINUE 972 350 CONTINUE 973 IF( IPACK.EQ.6 ) THEN 974 DO 370 JC = 1, UUB 975 DO 360 JR = 1, UUB + 1 - JC 976 A( JR, JC ) = ZERO 977 360 CONTINUE 978 370 CONTINUE 979 END IF 980 IF( IPACKG.EQ.5 ) THEN 981 IPACKG = IPACK 982 ELSE 983 IPACKG = 0 984 END IF 985 END IF 986 END IF 987 END IF 988* 989 ELSE 990* 991* 4) Generate Banded Matrix by first 992* Rotating by random Unitary matrices, 993* then reducing the bandwidth using Householder 994* transformations. 995* 996* Note: we should get here only if LDA .ge. N 997* 998 IF( ISYM.EQ.1 ) THEN 999* 1000* Non-symmetric -- A = U D V 1001* 1002 CALL DLAGGE( MR, NC, LLB, UUB, D, A, LDA, ISEED, WORK, 1003 $ IINFO ) 1004 ELSE 1005* 1006* Symmetric -- A = U D U' 1007* 1008 CALL DLAGSY( M, LLB, D, A, LDA, ISEED, WORK, IINFO ) 1009* 1010 END IF 1011 IF( IINFO.NE.0 ) THEN 1012 INFO = 3 1013 RETURN 1014 END IF 1015 END IF 1016* 1017* 5) Pack the matrix 1018* 1019 IF( IPACK.NE.IPACKG ) THEN 1020 IF( IPACK.EQ.1 ) THEN 1021* 1022* 'U' -- Upper triangular, not packed 1023* 1024 DO 390 J = 1, M 1025 DO 380 I = J + 1, M 1026 A( I, J ) = ZERO 1027 380 CONTINUE 1028 390 CONTINUE 1029* 1030 ELSE IF( IPACK.EQ.2 ) THEN 1031* 1032* 'L' -- Lower triangular, not packed 1033* 1034 DO 410 J = 2, M 1035 DO 400 I = 1, J - 1 1036 A( I, J ) = ZERO 1037 400 CONTINUE 1038 410 CONTINUE 1039* 1040 ELSE IF( IPACK.EQ.3 ) THEN 1041* 1042* 'C' -- Upper triangle packed Columnwise. 1043* 1044 ICOL = 1 1045 IROW = 0 1046 DO 430 J = 1, M 1047 DO 420 I = 1, J 1048 IROW = IROW + 1 1049 IF( IROW.GT.LDA ) THEN 1050 IROW = 1 1051 ICOL = ICOL + 1 1052 END IF 1053 A( IROW, ICOL ) = A( I, J ) 1054 420 CONTINUE 1055 430 CONTINUE 1056* 1057 ELSE IF( IPACK.EQ.4 ) THEN 1058* 1059* 'R' -- Lower triangle packed Columnwise. 1060* 1061 ICOL = 1 1062 IROW = 0 1063 DO 450 J = 1, M 1064 DO 440 I = J, M 1065 IROW = IROW + 1 1066 IF( IROW.GT.LDA ) THEN 1067 IROW = 1 1068 ICOL = ICOL + 1 1069 END IF 1070 A( IROW, ICOL ) = A( I, J ) 1071 440 CONTINUE 1072 450 CONTINUE 1073* 1074 ELSE IF( IPACK.GE.5 ) THEN 1075* 1076* 'B' -- The lower triangle is packed as a band matrix. 1077* 'Q' -- The upper triangle is packed as a band matrix. 1078* 'Z' -- The whole matrix is packed as a band matrix. 1079* 1080 IF( IPACK.EQ.5 ) 1081 $ UUB = 0 1082 IF( IPACK.EQ.6 ) 1083 $ LLB = 0 1084* 1085 DO 470 J = 1, UUB 1086 DO 460 I = MIN( J+LLB, M ), 1, -1 1087 A( I-J+UUB+1, J ) = A( I, J ) 1088 460 CONTINUE 1089 470 CONTINUE 1090* 1091 DO 490 J = UUB + 2, N 1092 DO 480 I = J - UUB, MIN( J+LLB, M ) 1093 A( I-J+UUB+1, J ) = A( I, J ) 1094 480 CONTINUE 1095 490 CONTINUE 1096 END IF 1097* 1098* If packed, zero out extraneous elements. 1099* 1100* Symmetric/Triangular Packed -- 1101* zero out everything after A(IROW,ICOL) 1102* 1103 IF( IPACK.EQ.3 .OR. IPACK.EQ.4 ) THEN 1104 DO 510 JC = ICOL, M 1105 DO 500 JR = IROW + 1, LDA 1106 A( JR, JC ) = ZERO 1107 500 CONTINUE 1108 IROW = 0 1109 510 CONTINUE 1110* 1111 ELSE IF( IPACK.GE.5 ) THEN 1112* 1113* Packed Band -- 1114* 1st row is now in A( UUB+2-j, j), zero above it 1115* m-th row is now in A( M+UUB-j,j), zero below it 1116* last non-zero diagonal is now in A( UUB+LLB+1,j ), 1117* zero below it, too. 1118* 1119 IR1 = UUB + LLB + 2 1120 IR2 = UUB + M + 2 1121 DO 540 JC = 1, N 1122 DO 520 JR = 1, UUB + 1 - JC 1123 A( JR, JC ) = ZERO 1124 520 CONTINUE 1125 DO 530 JR = MAX( 1, MIN( IR1, IR2-JC ) ), LDA 1126 A( JR, JC ) = ZERO 1127 530 CONTINUE 1128 540 CONTINUE 1129 END IF 1130 END IF 1131* 1132 RETURN 1133* 1134* End of DLATMT 1135* 1136 END 1137