1*> \brief \b DLATME 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 DLATME( N, DIST, ISEED, D, MODE, COND, DMAX, EI, 12* RSIGN, 13* UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, 14* A, 15* LDA, WORK, INFO ) 16* 17* .. Scalar Arguments .. 18* CHARACTER DIST, RSIGN, SIM, UPPER 19* INTEGER INFO, KL, KU, LDA, MODE, MODES, N 20* DOUBLE PRECISION ANORM, COND, CONDS, DMAX 21* .. 22* .. Array Arguments .. 23* CHARACTER EI( * ) 24* INTEGER ISEED( 4 ) 25* DOUBLE PRECISION A( LDA, * ), D( * ), DS( * ), WORK( * ) 26* .. 27* 28* 29*> \par Purpose: 30* ============= 31*> 32*> \verbatim 33*> 34*> DLATME generates random non-symmetric square matrices with 35*> specified eigenvalues for testing LAPACK programs. 36*> 37*> DLATME operates by applying the following sequence of 38*> operations: 39*> 40*> 1. Set the diagonal to D, where D may be input or 41*> computed according to MODE, COND, DMAX, and RSIGN 42*> as described below. 43*> 44*> 2. If complex conjugate pairs are desired (MODE=0 and EI(1)='R', 45*> or MODE=5), certain pairs of adjacent elements of D are 46*> interpreted as the real and complex parts of a complex 47*> conjugate pair; A thus becomes block diagonal, with 1x1 48*> and 2x2 blocks. 49*> 50*> 3. If UPPER='T', the upper triangle of A is set to random values 51*> out of distribution DIST. 52*> 53*> 4. If SIM='T', A is multiplied on the left by a random matrix 54*> X, whose singular values are specified by DS, MODES, and 55*> CONDS, and on the right by X inverse. 56*> 57*> 5. If KL < N-1, the lower bandwidth is reduced to KL using 58*> Householder transformations. If KU < N-1, the upper 59*> bandwidth is reduced to KU. 60*> 61*> 6. If ANORM is not negative, the matrix is scaled to have 62*> maximum-element-norm ANORM. 63*> 64*> (Note: since the matrix cannot be reduced beyond Hessenberg form, 65*> no packing options are available.) 66*> \endverbatim 67* 68* Arguments: 69* ========== 70* 71*> \param[in] N 72*> \verbatim 73*> N is INTEGER 74*> The number of columns (or rows) of A. Not modified. 75*> \endverbatim 76*> 77*> \param[in] DIST 78*> \verbatim 79*> DIST is CHARACTER*1 80*> On entry, DIST specifies the type of distribution to be used 81*> to generate the random eigen-/singular values, and for the 82*> upper triangle (see UPPER). 83*> 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform ) 84*> 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric ) 85*> 'N' => NORMAL( 0, 1 ) ( 'N' for normal ) 86*> Not modified. 87*> \endverbatim 88*> 89*> \param[in,out] ISEED 90*> \verbatim 91*> ISEED is INTEGER array, dimension ( 4 ) 92*> On entry ISEED specifies the seed of the random number 93*> generator. They should lie between 0 and 4095 inclusive, 94*> and ISEED(4) should be odd. The random number generator 95*> uses a linear congruential sequence limited to small 96*> integers, and so should produce machine independent 97*> random numbers. The values of ISEED are changed on 98*> exit, and can be used in the next call to DLATME 99*> to continue the same random number sequence. 100*> Changed on exit. 101*> \endverbatim 102*> 103*> \param[in,out] D 104*> \verbatim 105*> D is DOUBLE PRECISION array, dimension ( N ) 106*> This array is used to specify the eigenvalues of A. If 107*> MODE=0, then D is assumed to contain the eigenvalues (but 108*> see the description of EI), otherwise they will be 109*> computed according to MODE, COND, DMAX, and RSIGN and 110*> placed in D. 111*> Modified if MODE is nonzero. 112*> \endverbatim 113*> 114*> \param[in] MODE 115*> \verbatim 116*> MODE is INTEGER 117*> On entry this describes how the eigenvalues are to 118*> be specified: 119*> MODE = 0 means use D (with EI) as input 120*> MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND 121*> MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND 122*> MODE = 3 sets D(I)=COND**(-(I-1)/(N-1)) 123*> MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND) 124*> MODE = 5 sets D to random numbers in the range 125*> ( 1/COND , 1 ) such that their logarithms 126*> are uniformly distributed. Each odd-even pair 127*> of elements will be either used as two real 128*> eigenvalues or as the real and imaginary part 129*> of a complex conjugate pair of eigenvalues; 130*> the choice of which is done is random, with 131*> 50-50 probability, for each pair. 132*> MODE = 6 set D to random numbers from same distribution 133*> as the rest of the matrix. 134*> MODE < 0 has the same meaning as ABS(MODE), except that 135*> the order of the elements of D is reversed. 136*> Thus if MODE is between 1 and 4, D has entries ranging 137*> from 1 to 1/COND, if between -1 and -4, D has entries 138*> ranging from 1/COND to 1, 139*> Not modified. 140*> \endverbatim 141*> 142*> \param[in] COND 143*> \verbatim 144*> COND is DOUBLE PRECISION 145*> On entry, this is used as described under MODE above. 146*> If used, it must be >= 1. Not modified. 147*> \endverbatim 148*> 149*> \param[in] DMAX 150*> \verbatim 151*> DMAX is DOUBLE PRECISION 152*> If MODE is neither -6, 0 nor 6, the contents of D, as 153*> computed according to MODE and COND, will be scaled by 154*> DMAX / max(abs(D(i))). Note that DMAX need not be 155*> positive: if DMAX is negative (or zero), D will be 156*> scaled by a negative number (or zero). 157*> Not modified. 158*> \endverbatim 159*> 160*> \param[in] EI 161*> \verbatim 162*> EI is CHARACTER*1 array, dimension ( N ) 163*> If MODE is 0, and EI(1) is not ' ' (space character), 164*> this array specifies which elements of D (on input) are 165*> real eigenvalues and which are the real and imaginary parts 166*> of a complex conjugate pair of eigenvalues. The elements 167*> of EI may then only have the values 'R' and 'I'. If 168*> EI(j)='R' and EI(j+1)='I', then the j-th eigenvalue is 169*> CMPLX( D(j) , D(j+1) ), and the (j+1)-th is the complex 170*> conjugate thereof. If EI(j)=EI(j+1)='R', then the j-th 171*> eigenvalue is D(j) (i.e., real). EI(1) may not be 'I', 172*> nor may two adjacent elements of EI both have the value 'I'. 173*> If MODE is not 0, then EI is ignored. If MODE is 0 and 174*> EI(1)=' ', then the eigenvalues will all be real. 175*> Not modified. 176*> \endverbatim 177*> 178*> \param[in] RSIGN 179*> \verbatim 180*> RSIGN is CHARACTER*1 181*> If MODE is not 0, 6, or -6, and RSIGN='T', then the 182*> elements of D, as computed according to MODE and COND, will 183*> be multiplied by a random sign (+1 or -1). If RSIGN='F', 184*> they will not be. RSIGN may only have the values 'T' or 185*> 'F'. 186*> Not modified. 187*> \endverbatim 188*> 189*> \param[in] UPPER 190*> \verbatim 191*> UPPER is CHARACTER*1 192*> If UPPER='T', then the elements of A above the diagonal 193*> (and above the 2x2 diagonal blocks, if A has complex 194*> eigenvalues) will be set to random numbers out of DIST. 195*> If UPPER='F', they will not. UPPER may only have the 196*> values 'T' or 'F'. 197*> Not modified. 198*> \endverbatim 199*> 200*> \param[in] SIM 201*> \verbatim 202*> SIM is CHARACTER*1 203*> If SIM='T', then A will be operated on by a "similarity 204*> transform", i.e., multiplied on the left by a matrix X and 205*> on the right by X inverse. X = U S V, where U and V are 206*> random unitary matrices and S is a (diagonal) matrix of 207*> singular values specified by DS, MODES, and CONDS. If 208*> SIM='F', then A will not be transformed. 209*> Not modified. 210*> \endverbatim 211*> 212*> \param[in,out] DS 213*> \verbatim 214*> DS is DOUBLE PRECISION array, dimension ( N ) 215*> This array is used to specify the singular values of X, 216*> in the same way that D specifies the eigenvalues of A. 217*> If MODE=0, the DS contains the singular values, which 218*> may not be zero. 219*> Modified if MODE is nonzero. 220*> \endverbatim 221*> 222*> \param[in] MODES 223*> \verbatim 224*> MODES is INTEGER 225*> \endverbatim 226*> 227*> \param[in] CONDS 228*> \verbatim 229*> CONDS is DOUBLE PRECISION 230*> Same as MODE and COND, but for specifying the diagonal 231*> of S. MODES=-6 and +6 are not allowed (since they would 232*> result in randomly ill-conditioned eigenvalues.) 233*> \endverbatim 234*> 235*> \param[in] KL 236*> \verbatim 237*> KL is INTEGER 238*> This specifies the lower bandwidth of the matrix. KL=1 239*> specifies upper Hessenberg form. If KL is at least N-1, 240*> then A will have full lower bandwidth. KL must be at 241*> least 1. 242*> Not modified. 243*> \endverbatim 244*> 245*> \param[in] KU 246*> \verbatim 247*> KU is INTEGER 248*> This specifies the upper bandwidth of the matrix. KU=1 249*> specifies lower Hessenberg form. If KU is at least N-1, 250*> then A will have full upper bandwidth; if KU and KL 251*> are both at least N-1, then A will be dense. Only one of 252*> KU and KL may be less than N-1. KU must be at least 1. 253*> Not modified. 254*> \endverbatim 255*> 256*> \param[in] ANORM 257*> \verbatim 258*> ANORM is DOUBLE PRECISION 259*> If ANORM is not negative, then A will be scaled by a non- 260*> negative real number to make the maximum-element-norm of A 261*> to be ANORM. 262*> Not modified. 263*> \endverbatim 264*> 265*> \param[out] A 266*> \verbatim 267*> A is DOUBLE PRECISION array, dimension ( LDA, N ) 268*> On exit A is the desired test matrix. 269*> Modified. 270*> \endverbatim 271*> 272*> \param[in] LDA 273*> \verbatim 274*> LDA is INTEGER 275*> LDA specifies the first dimension of A as declared in the 276*> calling program. LDA must be at least N. 277*> Not modified. 278*> \endverbatim 279*> 280*> \param[out] WORK 281*> \verbatim 282*> WORK is DOUBLE PRECISION array, dimension ( 3*N ) 283*> Workspace. 284*> Modified. 285*> \endverbatim 286*> 287*> \param[out] INFO 288*> \verbatim 289*> INFO is INTEGER 290*> Error code. On exit, INFO will be set to one of the 291*> following values: 292*> 0 => normal return 293*> -1 => N negative 294*> -2 => DIST illegal string 295*> -5 => MODE not in range -6 to 6 296*> -6 => COND less than 1.0, and MODE neither -6, 0 nor 6 297*> -8 => EI(1) is not ' ' or 'R', EI(j) is not 'R' or 'I', or 298*> two adjacent elements of EI are 'I'. 299*> -9 => RSIGN is not 'T' or 'F' 300*> -10 => UPPER is not 'T' or 'F' 301*> -11 => SIM is not 'T' or 'F' 302*> -12 => MODES=0 and DS has a zero singular value. 303*> -13 => MODES is not in the range -5 to 5. 304*> -14 => MODES is nonzero and CONDS is less than 1. 305*> -15 => KL is less than 1. 306*> -16 => KU is less than 1, or KL and KU are both less than 307*> N-1. 308*> -19 => LDA is less than N. 309*> 1 => Error return from DLATM1 (computing D) 310*> 2 => Cannot scale to DMAX (max. eigenvalue is 0) 311*> 3 => Error return from DLATM1 (computing DS) 312*> 4 => Error return from DLARGE 313*> 5 => Zero singular value from DLATM1. 314*> \endverbatim 315* 316* Authors: 317* ======== 318* 319*> \author Univ. of Tennessee 320*> \author Univ. of California Berkeley 321*> \author Univ. of Colorado Denver 322*> \author NAG Ltd. 323* 324*> \date December 2016 325* 326*> \ingroup double_matgen 327* 328* ===================================================================== 329 SUBROUTINE DLATME( N, DIST, ISEED, D, MODE, COND, DMAX, EI, 330 $ RSIGN, 331 $ UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, 332 $ A, 333 $ LDA, WORK, INFO ) 334* 335* -- LAPACK computational routine (version 3.7.0) -- 336* -- LAPACK is a software package provided by Univ. of Tennessee, -- 337* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 338* December 2016 339* 340* .. Scalar Arguments .. 341 CHARACTER DIST, RSIGN, SIM, UPPER 342 INTEGER INFO, KL, KU, LDA, MODE, MODES, N 343 DOUBLE PRECISION ANORM, COND, CONDS, DMAX 344* .. 345* .. Array Arguments .. 346 CHARACTER EI( * ) 347 INTEGER ISEED( 4 ) 348 DOUBLE PRECISION A( LDA, * ), D( * ), DS( * ), WORK( * ) 349* .. 350* 351* ===================================================================== 352* 353* .. Parameters .. 354 DOUBLE PRECISION ZERO 355 PARAMETER ( ZERO = 0.0D0 ) 356 DOUBLE PRECISION ONE 357 PARAMETER ( ONE = 1.0D0 ) 358 DOUBLE PRECISION HALF 359 PARAMETER ( HALF = 1.0D0 / 2.0D0 ) 360* .. 361* .. Local Scalars .. 362 LOGICAL BADEI, BADS, USEEI 363 INTEGER I, IC, ICOLS, IDIST, IINFO, IR, IROWS, IRSIGN, 364 $ ISIM, IUPPER, J, JC, JCR, JR 365 DOUBLE PRECISION ALPHA, TAU, TEMP, XNORMS 366* .. 367* .. Local Arrays .. 368 DOUBLE PRECISION TEMPA( 1 ) 369* .. 370* .. External Functions .. 371 LOGICAL LSAME 372 DOUBLE PRECISION DLANGE, DLARAN 373 EXTERNAL LSAME, DLANGE, DLARAN 374* .. 375* .. External Subroutines .. 376 EXTERNAL DCOPY, DGEMV, DGER, DLARFG, DLARGE, DLARNV, 377 $ DLASET, DLATM1, DSCAL, XERBLA 378* .. 379* .. Intrinsic Functions .. 380 INTRINSIC ABS, MAX, MOD 381* .. 382* .. Executable Statements .. 383* 384* 1) Decode and Test the input parameters. 385* Initialize flags & seed. 386* 387 INFO = 0 388* 389* Quick return if possible 390* 391 IF( N.EQ.0 ) 392 $ RETURN 393* 394* Decode DIST 395* 396 IF( LSAME( DIST, 'U' ) ) THEN 397 IDIST = 1 398 ELSE IF( LSAME( DIST, 'S' ) ) THEN 399 IDIST = 2 400 ELSE IF( LSAME( DIST, 'N' ) ) THEN 401 IDIST = 3 402 ELSE 403 IDIST = -1 404 END IF 405* 406* Check EI 407* 408 USEEI = .TRUE. 409 BADEI = .FALSE. 410 IF( LSAME( EI( 1 ), ' ' ) .OR. MODE.NE.0 ) THEN 411 USEEI = .FALSE. 412 ELSE 413 IF( LSAME( EI( 1 ), 'R' ) ) THEN 414 DO 10 J = 2, N 415 IF( LSAME( EI( J ), 'I' ) ) THEN 416 IF( LSAME( EI( J-1 ), 'I' ) ) 417 $ BADEI = .TRUE. 418 ELSE 419 IF( .NOT.LSAME( EI( J ), 'R' ) ) 420 $ BADEI = .TRUE. 421 END IF 422 10 CONTINUE 423 ELSE 424 BADEI = .TRUE. 425 END IF 426 END IF 427* 428* Decode RSIGN 429* 430 IF( LSAME( RSIGN, 'T' ) ) THEN 431 IRSIGN = 1 432 ELSE IF( LSAME( RSIGN, 'F' ) ) THEN 433 IRSIGN = 0 434 ELSE 435 IRSIGN = -1 436 END IF 437* 438* Decode UPPER 439* 440 IF( LSAME( UPPER, 'T' ) ) THEN 441 IUPPER = 1 442 ELSE IF( LSAME( UPPER, 'F' ) ) THEN 443 IUPPER = 0 444 ELSE 445 IUPPER = -1 446 END IF 447* 448* Decode SIM 449* 450 IF( LSAME( SIM, 'T' ) ) THEN 451 ISIM = 1 452 ELSE IF( LSAME( SIM, 'F' ) ) THEN 453 ISIM = 0 454 ELSE 455 ISIM = -1 456 END IF 457* 458* Check DS, if MODES=0 and ISIM=1 459* 460 BADS = .FALSE. 461 IF( MODES.EQ.0 .AND. ISIM.EQ.1 ) THEN 462 DO 20 J = 1, N 463 IF( DS( J ).EQ.ZERO ) 464 $ BADS = .TRUE. 465 20 CONTINUE 466 END IF 467* 468* Set INFO if an error 469* 470 IF( N.LT.0 ) THEN 471 INFO = -1 472 ELSE IF( IDIST.EQ.-1 ) THEN 473 INFO = -2 474 ELSE IF( ABS( MODE ).GT.6 ) THEN 475 INFO = -5 476 ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE ) 477 $ THEN 478 INFO = -6 479 ELSE IF( BADEI ) THEN 480 INFO = -8 481 ELSE IF( IRSIGN.EQ.-1 ) THEN 482 INFO = -9 483 ELSE IF( IUPPER.EQ.-1 ) THEN 484 INFO = -10 485 ELSE IF( ISIM.EQ.-1 ) THEN 486 INFO = -11 487 ELSE IF( BADS ) THEN 488 INFO = -12 489 ELSE IF( ISIM.EQ.1 .AND. ABS( MODES ).GT.5 ) THEN 490 INFO = -13 491 ELSE IF( ISIM.EQ.1 .AND. MODES.NE.0 .AND. CONDS.LT.ONE ) THEN 492 INFO = -14 493 ELSE IF( KL.LT.1 ) THEN 494 INFO = -15 495 ELSE IF( KU.LT.1 .OR. ( KU.LT.N-1 .AND. KL.LT.N-1 ) ) THEN 496 INFO = -16 497 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 498 INFO = -19 499 END IF 500* 501 IF( INFO.NE.0 ) THEN 502 CALL XERBLA( 'DLATME', -INFO ) 503 RETURN 504 END IF 505* 506* Initialize random number generator 507* 508 DO 30 I = 1, 4 509 ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 ) 510 30 CONTINUE 511* 512 IF( MOD( ISEED( 4 ), 2 ).NE.1 ) 513 $ ISEED( 4 ) = ISEED( 4 ) + 1 514* 515* 2) Set up diagonal of A 516* 517* Compute D according to COND and MODE 518* 519 CALL DLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, N, IINFO ) 520 IF( IINFO.NE.0 ) THEN 521 INFO = 1 522 RETURN 523 END IF 524 IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN 525* 526* Scale by DMAX 527* 528 TEMP = ABS( D( 1 ) ) 529 DO 40 I = 2, N 530 TEMP = MAX( TEMP, ABS( D( I ) ) ) 531 40 CONTINUE 532* 533 IF( TEMP.GT.ZERO ) THEN 534 ALPHA = DMAX / TEMP 535 ELSE IF( DMAX.NE.ZERO ) THEN 536 INFO = 2 537 RETURN 538 ELSE 539 ALPHA = ZERO 540 END IF 541* 542 CALL DSCAL( N, ALPHA, D, 1 ) 543* 544 END IF 545* 546 CALL DLASET( 'Full', N, N, ZERO, ZERO, A, LDA ) 547 CALL DCOPY( N, D, 1, A, LDA+1 ) 548* 549* Set up complex conjugate pairs 550* 551 IF( MODE.EQ.0 ) THEN 552 IF( USEEI ) THEN 553 DO 50 J = 2, N 554 IF( LSAME( EI( J ), 'I' ) ) THEN 555 A( J-1, J ) = A( J, J ) 556 A( J, J-1 ) = -A( J, J ) 557 A( J, J ) = A( J-1, J-1 ) 558 END IF 559 50 CONTINUE 560 END IF 561* 562 ELSE IF( ABS( MODE ).EQ.5 ) THEN 563* 564 DO 60 J = 2, N, 2 565 IF( DLARAN( ISEED ).GT.HALF ) THEN 566 A( J-1, J ) = A( J, J ) 567 A( J, J-1 ) = -A( J, J ) 568 A( J, J ) = A( J-1, J-1 ) 569 END IF 570 60 CONTINUE 571 END IF 572* 573* 3) If UPPER='T', set upper triangle of A to random numbers. 574* (but don't modify the corners of 2x2 blocks.) 575* 576 IF( IUPPER.NE.0 ) THEN 577 DO 70 JC = 2, N 578 IF( A( JC-1, JC ).NE.ZERO ) THEN 579 JR = JC - 2 580 ELSE 581 JR = JC - 1 582 END IF 583 CALL DLARNV( IDIST, ISEED, JR, A( 1, JC ) ) 584 70 CONTINUE 585 END IF 586* 587* 4) If SIM='T', apply similarity transformation. 588* 589* -1 590* Transform is X A X , where X = U S V, thus 591* 592* it is U S V A V' (1/S) U' 593* 594 IF( ISIM.NE.0 ) THEN 595* 596* Compute S (singular values of the eigenvector matrix) 597* according to CONDS and MODES 598* 599 CALL DLATM1( MODES, CONDS, 0, 0, ISEED, DS, N, IINFO ) 600 IF( IINFO.NE.0 ) THEN 601 INFO = 3 602 RETURN 603 END IF 604* 605* Multiply by V and V' 606* 607 CALL DLARGE( N, A, LDA, ISEED, WORK, IINFO ) 608 IF( IINFO.NE.0 ) THEN 609 INFO = 4 610 RETURN 611 END IF 612* 613* Multiply by S and (1/S) 614* 615 DO 80 J = 1, N 616 CALL DSCAL( N, DS( J ), A( J, 1 ), LDA ) 617 IF( DS( J ).NE.ZERO ) THEN 618 CALL DSCAL( N, ONE / DS( J ), A( 1, J ), 1 ) 619 ELSE 620 INFO = 5 621 RETURN 622 END IF 623 80 CONTINUE 624* 625* Multiply by U and U' 626* 627 CALL DLARGE( N, A, LDA, ISEED, WORK, IINFO ) 628 IF( IINFO.NE.0 ) THEN 629 INFO = 4 630 RETURN 631 END IF 632 END IF 633* 634* 5) Reduce the bandwidth. 635* 636 IF( KL.LT.N-1 ) THEN 637* 638* Reduce bandwidth -- kill column 639* 640 DO 90 JCR = KL + 1, N - 1 641 IC = JCR - KL 642 IROWS = N + 1 - JCR 643 ICOLS = N + KL - JCR 644* 645 CALL DCOPY( IROWS, A( JCR, IC ), 1, WORK, 1 ) 646 XNORMS = WORK( 1 ) 647 CALL DLARFG( IROWS, XNORMS, WORK( 2 ), 1, TAU ) 648 WORK( 1 ) = ONE 649* 650 CALL DGEMV( 'T', IROWS, ICOLS, ONE, A( JCR, IC+1 ), LDA, 651 $ WORK, 1, ZERO, WORK( IROWS+1 ), 1 ) 652 CALL DGER( IROWS, ICOLS, -TAU, WORK, 1, WORK( IROWS+1 ), 1, 653 $ A( JCR, IC+1 ), LDA ) 654* 655 CALL DGEMV( 'N', N, IROWS, ONE, A( 1, JCR ), LDA, WORK, 1, 656 $ ZERO, WORK( IROWS+1 ), 1 ) 657 CALL DGER( N, IROWS, -TAU, WORK( IROWS+1 ), 1, WORK, 1, 658 $ A( 1, JCR ), LDA ) 659* 660 A( JCR, IC ) = XNORMS 661 CALL DLASET( 'Full', IROWS-1, 1, ZERO, ZERO, A( JCR+1, IC ), 662 $ LDA ) 663 90 CONTINUE 664 ELSE IF( KU.LT.N-1 ) THEN 665* 666* Reduce upper bandwidth -- kill a row at a time. 667* 668 DO 100 JCR = KU + 1, N - 1 669 IR = JCR - KU 670 IROWS = N + KU - JCR 671 ICOLS = N + 1 - JCR 672* 673 CALL DCOPY( ICOLS, A( IR, JCR ), LDA, WORK, 1 ) 674 XNORMS = WORK( 1 ) 675 CALL DLARFG( ICOLS, XNORMS, WORK( 2 ), 1, TAU ) 676 WORK( 1 ) = ONE 677* 678 CALL DGEMV( 'N', IROWS, ICOLS, ONE, A( IR+1, JCR ), LDA, 679 $ WORK, 1, ZERO, WORK( ICOLS+1 ), 1 ) 680 CALL DGER( IROWS, ICOLS, -TAU, WORK( ICOLS+1 ), 1, WORK, 1, 681 $ A( IR+1, JCR ), LDA ) 682* 683 CALL DGEMV( 'C', ICOLS, N, ONE, A( JCR, 1 ), LDA, WORK, 1, 684 $ ZERO, WORK( ICOLS+1 ), 1 ) 685 CALL DGER( ICOLS, N, -TAU, WORK, 1, WORK( ICOLS+1 ), 1, 686 $ A( JCR, 1 ), LDA ) 687* 688 A( IR, JCR ) = XNORMS 689 CALL DLASET( 'Full', 1, ICOLS-1, ZERO, ZERO, A( IR, JCR+1 ), 690 $ LDA ) 691 100 CONTINUE 692 END IF 693* 694* Scale the matrix to have norm ANORM 695* 696 IF( ANORM.GE.ZERO ) THEN 697 TEMP = DLANGE( 'M', N, N, A, LDA, TEMPA ) 698 IF( TEMP.GT.ZERO ) THEN 699 ALPHA = ANORM / TEMP 700 DO 110 J = 1, N 701 CALL DSCAL( N, ALPHA, A( 1, J ), 1 ) 702 110 CONTINUE 703 END IF 704 END IF 705* 706 RETURN 707* 708* End of DLATME 709* 710 END 711