1 SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) 2* 3* -- LAPACK driver routine (version 1.1) -- 4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5* Courant Institute, Argonne National Lab, and Rice University 6* March 31, 1993 7* 8* .. Scalar Arguments .. 9 CHARACTER JOBZ, UPLO 10 INTEGER INFO, LDA, LWORK, N 11* .. 12* .. Array Arguments .. 13 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ) 14* .. 15* 16* Purpose 17* ======= 18* 19* DSYEV computes all eigenvalues and, optionally, eigenvectors of a 20* real symmetric matrix A. 21* 22* Arguments 23* ========= 24* 25* JOBZ (input) CHARACTER*1 26* = 'N': Compute eigenvalues only; 27* = 'V': Compute eigenvalues and eigenvectors. 28* 29* UPLO (input) CHARACTER*1 30* = 'U': Upper triangle of A is stored; 31* = 'L': Lower triangle of A is stored. 32* 33* N (input) INTEGER 34* The order of the matrix A. N >= 0. 35* 36* A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 37* On entry, the symmetric matrix A. If UPLO = 'U', the 38* leading N-by-N upper triangular part of A contains the 39* upper triangular part of the matrix A. If UPLO = 'L', 40* the leading N-by-N lower triangular part of A contains 41* the lower triangular part of the matrix A. 42* On exit, if JOBZ = 'V', then if INFO = 0, A contains the 43* orthonormal eigenvectors of the matrix A. 44* If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 45* or the upper triangle (if UPLO='U') of A, including the 46* diagonal, is destroyed. 47* 48* LDA (input) INTEGER 49* The leading dimension of the array A. LDA >= max(1,N). 50* 51* W (output) DOUBLE PRECISION array, dimension (N) 52* If INFO = 0, the eigenvalues in ascending order. 53* 54* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 55* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 56* 57* LWORK (input) INTEGER 58* The length of the array WORK. LWORK >= max(1,3*N-1). 59* For optimal efficiency, LWORK >= (NB+2)*N, 60* where NB is the blocksize for DSYTRD returned by ILAENV. 61* 62* INFO (output) INTEGER 63* = 0: successful exit 64* < 0: if INFO = -i, the i-th argument had an illegal value 65* > 0: if INFO = i, the algorithm failed to converge; i 66* off-diagonal elements of an intermediate tridiagonal 67* form did not converge to zero. 68* 69* ===================================================================== 70* 71* .. Parameters .. 72 DOUBLE PRECISION ZERO, ONE 73 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 74* .. 75* .. Local Scalars .. 76 LOGICAL LOWER, WANTZ 77 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE, J, 78 $ LLWORK, LOPT 79 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 80 $ SMLNUM 81* .. 82* .. External Functions .. 83 LOGICAL LSAME 84 DOUBLE PRECISION DLAMCH, DLANSY 85 EXTERNAL LSAME, DLAMCH, DLANSY 86* .. 87* .. External Subroutines .. 88 EXTERNAL DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD, XERBLA 89* .. 90* .. Intrinsic Functions .. 91 INTRINSIC MAX, SQRT 92* .. 93* .. Executable Statements .. 94* 95* Test the input parameters. 96* 97 WANTZ = LSAME( JOBZ, 'V' ) 98 LOWER = LSAME( UPLO, 'L' ) 99* 100 INFO = 0 101 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 102 INFO = -1 103 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 104 INFO = -2 105 ELSE IF( N.LT.0 ) THEN 106 INFO = -3 107 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 108 INFO = -5 109 ELSE IF( LWORK.LT.MAX( 1, 3*N-1 ) ) THEN 110 INFO = -8 111 END IF 112* 113 IF( INFO.NE.0 ) THEN 114 CALL XERBLA( 'DSYEV ', -INFO ) 115 RETURN 116 END IF 117* 118* Quick return if possible 119* 120 IF( N.EQ.0 ) THEN 121 WORK( 1 ) = 1 122 RETURN 123 END IF 124* 125 IF( N.EQ.1 ) THEN 126 W( 1 ) = A( 1, 1 ) 127 WORK( 1 ) = 3 128 IF( WANTZ ) 129 $ A( 1, 1 ) = ONE 130 RETURN 131 END IF 132* 133* Get machine constants. 134* 135 SAFMIN = DLAMCH( 'Safe minimum' ) 136 EPS = DLAMCH( 'Precision' ) 137 SMLNUM = SAFMIN / EPS 138 BIGNUM = ONE / SMLNUM 139 RMIN = SQRT( SMLNUM ) 140 RMAX = SQRT( BIGNUM ) 141* 142* Scale matrix to allowable range, if necessary. 143* 144 ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK ) 145 ISCALE = 0 146 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 147 ISCALE = 1 148 SIGMA = RMIN / ANRM 149 ELSE IF( ANRM.GT.RMAX ) THEN 150 ISCALE = 1 151 SIGMA = RMAX / ANRM 152 END IF 153 IF( ISCALE.EQ.1 ) THEN 154 IF( LOWER ) THEN 155 DO 10 J = 1, N 156 CALL DSCAL( N-J+1, SIGMA, A( J, J ), 1 ) 157 10 CONTINUE 158 ELSE 159 DO 20 J = 1, N 160 CALL DSCAL( J, SIGMA, A( 1, J ), 1 ) 161 20 CONTINUE 162 END IF 163 END IF 164* 165* Call DSYTRD to reduce symmetric matrix to tridiagonal form. 166* 167 INDE = 1 168 INDTAU = INDE + N 169 INDWRK = INDTAU + N 170 LLWORK = LWORK - INDWRK + 1 171 CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ), 172 $ WORK( INDWRK ), LLWORK, IINFO ) 173 LOPT = 2*N + WORK( INDWRK ) 174* 175* For eigenvalues only, call DSTERF. For eigenvectors, first call 176* DORGTR to generate the orthogonal matrix, then call DSTEQR. 177* 178 IF( .NOT.WANTZ ) THEN 179 CALL DSTERF( N, W, WORK( INDE ), INFO ) 180 ELSE 181 CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ), 182 $ LLWORK, IINFO ) 183 CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ), 184 $ INFO ) 185 END IF 186* 187* If matrix was scaled, then rescale eigenvalues appropriately. 188* 189 IF( ISCALE.EQ.1 ) THEN 190 IF( INFO.EQ.0 ) THEN 191 IMAX = N 192 ELSE 193 IMAX = INFO - 1 194 END IF 195 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 196 END IF 197* 198* Set WORK(1) to optimal workspace size. 199* 200 WORK( 1 ) = MAX( 3*N-1, LOPT ) 201* 202 RETURN 203* 204* End of DSYEV 205* 206 END 207 SUBROUTINE DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO ) 208* 209* -- LAPACK routine (version 1.1) -- 210* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 211* Courant Institute, Argonne National Lab, and Rice University 212* March 31, 1993 213* 214* .. Scalar Arguments .. 215 CHARACTER COMPZ 216 INTEGER INFO, LDZ, N 217* .. 218* .. Array Arguments .. 219 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * ) 220* .. 221* 222* Purpose 223* ======= 224* 225* DSTEQR computes all eigenvalues and, optionally, eigenvectors of a 226* symmetric tridiagonal matrix using the implicit QL or QR method. 227* The eigenvectors of a full or band symmetric matrix can also be found 228* if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to 229* tridiagonal form. 230* 231* Arguments 232* ========= 233* 234* COMPZ (input) CHARACTER*1 235* = 'N': Compute eigenvalues only. 236* = 'V': Compute eigenvalues and eigenvectors of the original 237* symmetric matrix. On entry, Z must contain the 238* orthogonal matrix used to reduce the original matrix 239* to tridiagonal form. 240* = 'I': Compute eigenvalues and eigenvectors of the 241* tridiagonal matrix. Z is initialized to the identity 242* matrix. 243* 244* N (input) INTEGER 245* The order of the matrix. N >= 0. 246* 247* D (input/output) DOUBLE PRECISION array, dimension (N) 248* On entry, the diagonal elements of the tridiagonal matrix. 249* On exit, if INFO = 0, the eigenvalues in ascending order. 250* 251* E (input/output) DOUBLE PRECISION array, dimension (N-1) 252* On entry, the (n-1) subdiagonal elements of the tridiagonal 253* matrix. 254* On exit, E has been destroyed. 255* 256* Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) 257* On entry, if COMPZ = 'V', then Z contains the orthogonal 258* matrix used in the reduction to tridiagonal form. 259* On exit, if COMPZ = 'V', Z contains the orthonormal 260* eigenvectors of the original symmetric matrix, and if 261* COMPZ = 'I', Z contains the orthonormal eigenvectors of 262* the symmetric tridiagonal matrix. If an error exit is 263* made, Z contains the eigenvectors associated with the 264* stored eigenvalues. 265* If COMPZ = 'N', then Z is not referenced. 266* 267* LDZ (input) INTEGER 268* The leading dimension of the array Z. LDZ >= 1, and if 269* eigenvectors are desired, then LDZ >= max(1,N). 270* 271* WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2)) 272* If COMPZ = 'N', then WORK is not referenced. 273* 274* INFO (output) INTEGER 275* = 0: successful exit 276* < 0: if INFO = -i, the i-th argument had an illegal value 277* > 0: the algorithm has failed to find all the eigenvalues in 278* a total of 30*N iterations; if INFO = i, then i 279* elements of E have not converged to zero; on exit, D 280* and E contain the elements of a symmetric tridiagonal 281* matrix which is orthogonally similar to the original 282* matrix. 283* 284* ===================================================================== 285* 286* .. Parameters .. 287 DOUBLE PRECISION ZERO, ONE, TWO 288 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 289 INTEGER MAXIT 290 PARAMETER ( MAXIT = 30 ) 291* .. 292* .. Local Scalars .. 293 INTEGER I, ICOMPZ, II, J, JTOT, K, L, L1, LEND, LENDM1, 294 $ LENDP1, LM1, M, MM, MM1, NM1, NMAXIT 295 DOUBLE PRECISION B, C, EPS, F, G, P, R, RT1, RT2, S, TST 296* .. 297* .. External Functions .. 298 LOGICAL LSAME 299 DOUBLE PRECISION DLAMCH, DLAPY2 300 EXTERNAL LSAME, DLAMCH, DLAPY2 301* .. 302* .. External Subroutines .. 303 EXTERNAL DLAE2, DLAEV2, DLARTG, DLASR, DLAZRO, DSWAP, 304 $ XERBLA 305* .. 306* .. Intrinsic Functions .. 307 INTRINSIC ABS, MAX, SIGN 308* .. 309* .. Executable Statements .. 310* 311* Test the input parameters. 312* 313 INFO = 0 314* 315 IF( LSAME( COMPZ, 'N' ) ) THEN 316 ICOMPZ = 0 317 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 318 ICOMPZ = 1 319 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 320 ICOMPZ = 2 321 ELSE 322 ICOMPZ = -1 323 END IF 324 IF( ICOMPZ.LT.0 ) THEN 325 INFO = -1 326 ELSE IF( N.LT.0 ) THEN 327 INFO = -2 328 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, 329 $ N ) ) ) THEN 330 INFO = -6 331 END IF 332 IF( INFO.NE.0 ) THEN 333 CALL XERBLA( 'DSTEQR', -INFO ) 334 RETURN 335 END IF 336* 337* Quick return if possible 338* 339 IF( N.EQ.0 ) 340 $ RETURN 341* 342 IF( N.EQ.1 ) THEN 343 IF( ICOMPZ.GT.0 ) 344 $ Z( 1, 1 ) = ONE 345 RETURN 346 END IF 347* 348* Determine the unit roundoff for this environment. 349* 350 EPS = DLAMCH( 'E' ) 351* 352* Compute the eigenvalues and eigenvectors of the tridiagonal 353* matrix. 354* 355 IF( ICOMPZ.EQ.2 ) 356 $ CALL DLAZRO( N, N, ZERO, ONE, Z, LDZ ) 357* 358 NMAXIT = N*MAXIT 359 JTOT = 0 360* 361* Determine where the matrix splits and choose QL or QR iteration 362* for each block, according to whether top or bottom diagonal 363* element is smaller. 364* 365 L1 = 1 366 NM1 = N - 1 367* 368 10 CONTINUE 369 IF( L1.GT.N ) 370 $ GO TO 160 371 IF( L1.GT.1 ) 372 $ E( L1-1 ) = ZERO 373 IF( L1.LE.NM1 ) THEN 374 DO 20 M = L1, NM1 375 TST = ABS( E( M ) ) 376 IF( TST.LE.EPS*( ABS( D( M ) )+ABS( D( M+1 ) ) ) ) 377 $ GO TO 30 378 20 CONTINUE 379 END IF 380 M = N 381* 382 30 CONTINUE 383 L = L1 384 LEND = M 385 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 386 L = LEND 387 LEND = L1 388 END IF 389 L1 = M + 1 390* 391 IF( LEND.GE.L ) THEN 392* 393* QL Iteration 394* 395* Look for small subdiagonal element. 396* 397 40 CONTINUE 398 IF( L.NE.LEND ) THEN 399 LENDM1 = LEND - 1 400 DO 50 M = L, LENDM1 401 TST = ABS( E( M ) ) 402 IF( TST.LE.EPS*( ABS( D( M ) )+ABS( D( M+1 ) ) ) ) 403 $ GO TO 60 404 50 CONTINUE 405 END IF 406* 407 M = LEND 408* 409 60 CONTINUE 410 IF( M.LT.LEND ) 411 $ E( M ) = ZERO 412 P = D( L ) 413 IF( M.EQ.L ) 414 $ GO TO 80 415* 416* If remaining matrix is 2-by-2, use DLAE2 or DLAEV2 417* to compute its eigensystem. 418* 419 IF( M.EQ.L+1 ) THEN 420 IF( ICOMPZ.GT.0 ) THEN 421 CALL DLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S ) 422 WORK( L ) = C 423 WORK( N-1+L ) = S 424 CALL DLASR( 'R', 'V', 'B', N, 2, WORK( L ), 425 $ WORK( N-1+L ), Z( 1, L ), LDZ ) 426 ELSE 427 CALL DLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 ) 428 END IF 429 D( L ) = RT1 430 D( L+1 ) = RT2 431 E( L ) = ZERO 432 L = L + 2 433 IF( L.LE.LEND ) 434 $ GO TO 40 435 GO TO 10 436 END IF 437* 438 IF( JTOT.EQ.NMAXIT ) 439 $ GO TO 140 440 JTOT = JTOT + 1 441* 442* Form shift. 443* 444 G = ( D( L+1 )-P ) / ( TWO*E( L ) ) 445 R = DLAPY2( G, ONE ) 446 G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) ) 447* 448 S = ONE 449 C = ONE 450 P = ZERO 451* 452* Inner loop 453* 454 MM1 = M - 1 455 DO 70 I = MM1, L, -1 456 F = S*E( I ) 457 B = C*E( I ) 458 CALL DLARTG( G, F, C, S, R ) 459 IF( I.NE.M-1 ) 460 $ E( I+1 ) = R 461 G = D( I+1 ) - P 462 R = ( D( I )-G )*S + TWO*C*B 463 P = S*R 464 D( I+1 ) = G + P 465 G = C*R - B 466* 467* If eigenvectors are desired, then save rotations. 468* 469 IF( ICOMPZ.GT.0 ) THEN 470 WORK( I ) = C 471 WORK( N-1+I ) = -S 472 END IF 473* 474 70 CONTINUE 475* 476* If eigenvectors are desired, then apply saved rotations. 477* 478 IF( ICOMPZ.GT.0 ) THEN 479 MM = M - L + 1 480 CALL DLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ), 481 $ Z( 1, L ), LDZ ) 482 END IF 483* 484 D( L ) = D( L ) - P 485 E( L ) = G 486 GO TO 40 487* 488* Eigenvalue found. 489* 490 80 CONTINUE 491 D( L ) = P 492* 493 L = L + 1 494 IF( L.LE.LEND ) 495 $ GO TO 40 496 GO TO 10 497* 498 ELSE 499* 500* QR Iteration 501* 502* Look for small superdiagonal element. 503* 504 90 CONTINUE 505 IF( L.NE.LEND ) THEN 506 LENDP1 = LEND + 1 507 DO 100 M = L, LENDP1, -1 508 TST = ABS( E( M-1 ) ) 509 IF( TST.LE.EPS*( ABS( D( M ) )+ABS( D( M-1 ) ) ) ) 510 $ GO TO 110 511 100 CONTINUE 512 END IF 513* 514 M = LEND 515* 516 110 CONTINUE 517 IF( M.GT.LEND ) 518 $ E( M-1 ) = ZERO 519 P = D( L ) 520 IF( M.EQ.L ) 521 $ GO TO 130 522* 523* If remaining matrix is 2-by-2, use DLAE2 or DLAEV2 524* to compute its eigensystem. 525* 526 IF( M.EQ.L-1 ) THEN 527 IF( ICOMPZ.GT.0 ) THEN 528 CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S ) 529 WORK( M ) = C 530 WORK( N-1+M ) = S 531 CALL DLASR( 'R', 'V', 'F', N, 2, WORK( M ), 532 $ WORK( N-1+M ), Z( 1, L-1 ), LDZ ) 533 ELSE 534 CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 ) 535 END IF 536 D( L-1 ) = RT1 537 D( L ) = RT2 538 E( L-1 ) = ZERO 539 L = L - 2 540 IF( L.GE.LEND ) 541 $ GO TO 90 542 GO TO 10 543 END IF 544* 545 IF( JTOT.EQ.NMAXIT ) 546 $ GO TO 140 547 JTOT = JTOT + 1 548* 549* Form shift. 550* 551 G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ) 552 R = DLAPY2( G, ONE ) 553 G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) ) 554* 555 S = ONE 556 C = ONE 557 P = ZERO 558* 559* Inner loop 560* 561 LM1 = L - 1 562 DO 120 I = M, LM1 563 F = S*E( I ) 564 B = C*E( I ) 565 CALL DLARTG( G, F, C, S, R ) 566 IF( I.NE.M ) 567 $ E( I-1 ) = R 568 G = D( I ) - P 569 R = ( D( I+1 )-G )*S + TWO*C*B 570 P = S*R 571 D( I ) = G + P 572 G = C*R - B 573* 574* If eigenvectors are desired, then save rotations. 575* 576 IF( ICOMPZ.GT.0 ) THEN 577 WORK( I ) = C 578 WORK( N-1+I ) = S 579 END IF 580* 581 120 CONTINUE 582* 583* If eigenvectors are desired, then apply saved rotations. 584* 585 IF( ICOMPZ.GT.0 ) THEN 586 MM = L - M + 1 587 CALL DLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ), 588 $ Z( 1, M ), LDZ ) 589 END IF 590* 591 D( L ) = D( L ) - P 592 E( LM1 ) = G 593 GO TO 90 594* 595* Eigenvalue found. 596* 597 130 CONTINUE 598 D( L ) = P 599* 600 L = L - 1 601 IF( L.GE.LEND ) 602 $ GO TO 90 603 GO TO 10 604* 605 END IF 606* 607* Set error -- no convergence to an eigenvalue after a total 608* of N*MAXIT iterations. 609* 610 140 CONTINUE 611 DO 150 I = 1, N - 1 612 IF( E( I ).NE.ZERO ) 613 $ INFO = INFO + 1 614 150 CONTINUE 615 RETURN 616* 617* Order eigenvalues and eigenvectors. 618* 619 160 CONTINUE 620 DO 180 II = 2, N 621 I = II - 1 622 K = I 623 P = D( I ) 624 DO 170 J = II, N 625 IF( D( J ).LT.P ) THEN 626 K = J 627 P = D( J ) 628 END IF 629 170 CONTINUE 630 IF( K.NE.I ) THEN 631 D( K ) = D( I ) 632 D( I ) = P 633 IF( ICOMPZ.GT.0 ) 634 $ CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) 635 END IF 636 180 CONTINUE 637* 638 RETURN 639* 640* End of DSTEQR 641* 642 END 643 SUBROUTINE DLARTG( F, G, CS, SN, R ) 644* 645* -- LAPACK auxiliary routine (version 1.1) -- 646* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 647* Courant Institute, Argonne National Lab, and Rice University 648* October 31, 1992 649* 650* .. Scalar Arguments .. 651 DOUBLE PRECISION CS, F, G, R, SN 652* .. 653* 654* Purpose 655* ======= 656* 657* DLARTG generate a plane rotation so that 658* 659* [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1. 660* [ -SN CS ] [ G ] [ 0 ] 661* 662* This is a faster version of the BLAS1 routine DROTG, except for 663* the following differences: 664* F and G are unchanged on return. 665* If G=0, then CS=1 and SN=0. 666* If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any 667* floating point operations (saves work in DBDSQR when 668* there are zeros on the diagonal). 669* 670* Arguments 671* ========= 672* 673* F (input) DOUBLE PRECISION 674* The first component of vector to be rotated. 675* 676* G (input) DOUBLE PRECISION 677* The second component of vector to be rotated. 678* 679* CS (output) DOUBLE PRECISION 680* The cosine of the rotation. 681* 682* SN (output) DOUBLE PRECISION 683* The sine of the rotation. 684* 685* R (output) DOUBLE PRECISION 686* The nonzero component of the rotated vector. 687* 688* ===================================================================== 689* 690* .. Parameters .. 691 DOUBLE PRECISION ZERO 692 PARAMETER ( ZERO = 0.0D0 ) 693 DOUBLE PRECISION ONE 694 PARAMETER ( ONE = 1.0D0 ) 695* .. 696* .. Local Scalars .. 697 DOUBLE PRECISION T, TT 698* .. 699* .. Intrinsic Functions .. 700 INTRINSIC ABS, SQRT 701* .. 702* .. Executable Statements .. 703* 704 IF( G.EQ.ZERO ) THEN 705 CS = ONE 706 SN = ZERO 707 R = F 708 ELSE IF( F.EQ.ZERO ) THEN 709 CS = ZERO 710 SN = ONE 711 R = G 712 ELSE 713 IF( ABS( F ).GT.ABS( G ) ) THEN 714 T = G / F 715 TT = SQRT( ONE+T*T ) 716 CS = ONE / TT 717 SN = T*CS 718 R = F*TT 719 ELSE 720 T = F / G 721 TT = SQRT( ONE+T*T ) 722 SN = ONE / TT 723 CS = T*SN 724 R = G*TT 725 END IF 726 END IF 727 RETURN 728* 729* End of DLARTG 730* 731 END 732 SUBROUTINE DLASR( SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA ) 733* 734* -- LAPACK auxiliary routine (version 1.1) -- 735* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 736* Courant Institute, Argonne National Lab, and Rice University 737* October 31, 1992 738* 739* .. Scalar Arguments .. 740 CHARACTER DIRECT, PIVOT, SIDE 741 INTEGER LDA, M, N 742* .. 743* .. Array Arguments .. 744 DOUBLE PRECISION A( LDA, * ), C( * ), S( * ) 745* .. 746* 747* Purpose 748* ======= 749* 750* DLASR performs the transformation 751* 752* A := P*A, when SIDE = 'L' or 'l' ( Left-hand side ) 753* 754* A := A*P', when SIDE = 'R' or 'r' ( Right-hand side ) 755* 756* where A is an m by n real matrix and P is an orthogonal matrix, 757* consisting of a sequence of plane rotations determined by the 758* parameters PIVOT and DIRECT as follows ( z = m when SIDE = 'L' or 'l' 759* and z = n when SIDE = 'R' or 'r' ): 760* 761* When DIRECT = 'F' or 'f' ( Forward sequence ) then 762* 763* P = P( z - 1 )*...*P( 2 )*P( 1 ), 764* 765* and when DIRECT = 'B' or 'b' ( Backward sequence ) then 766* 767* P = P( 1 )*P( 2 )*...*P( z - 1 ), 768* 769* where P( k ) is a plane rotation matrix for the following planes: 770* 771* when PIVOT = 'V' or 'v' ( Variable pivot ), 772* the plane ( k, k + 1 ) 773* 774* when PIVOT = 'T' or 't' ( Top pivot ), 775* the plane ( 1, k + 1 ) 776* 777* when PIVOT = 'B' or 'b' ( Bottom pivot ), 778* the plane ( k, z ) 779* 780* c( k ) and s( k ) must contain the cosine and sine that define the 781* matrix P( k ). The two by two plane rotation part of the matrix 782* P( k ), R( k ), is assumed to be of the form 783* 784* R( k ) = ( c( k ) s( k ) ). 785* ( -s( k ) c( k ) ) 786* 787* This version vectorises across rows of the array A when SIDE = 'L'. 788* 789* Arguments 790* ========= 791* 792* SIDE (input) CHARACTER*1 793* Specifies whether the plane rotation matrix P is applied to 794* A on the left or the right. 795* = 'L': Left, compute A := P*A 796* = 'R': Right, compute A:= A*P' 797* 798* DIRECT (input) CHARACTER*1 799* Specifies whether P is a forward or backward sequence of 800* plane rotations. 801* = 'F': Forward, P = P( z - 1 )*...*P( 2 )*P( 1 ) 802* = 'B': Backward, P = P( 1 )*P( 2 )*...*P( z - 1 ) 803* 804* PIVOT (input) CHARACTER*1 805* Specifies the plane for which P(k) is a plane rotation 806* matrix. 807* = 'V': Variable pivot, the plane (k,k+1) 808* = 'T': Top pivot, the plane (1,k+1) 809* = 'B': Bottom pivot, the plane (k,z) 810* 811* M (input) INTEGER 812* The number of rows of the matrix A. If m <= 1, an immediate 813* return is effected. 814* 815* N (input) INTEGER 816* The number of columns of the matrix A. If n <= 1, an 817* immediate return is effected. 818* 819* C, S (input) DOUBLE PRECISION arrays, dimension 820* (M-1) if SIDE = 'L' 821* (N-1) if SIDE = 'R' 822* c(k) and s(k) contain the cosine and sine that define the 823* matrix P(k). The two by two plane rotation part of the 824* matrix P(k), R(k), is assumed to be of the form 825* R( k ) = ( c( k ) s( k ) ). 826* ( -s( k ) c( k ) ) 827* 828* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 829* The m by n matrix A. On exit, A is overwritten by P*A if 830* SIDE = 'R' or by A*P' if SIDE = 'L'. 831* 832* LDA (input) INTEGER 833* The leading dimension of the array A. LDA >= max(1,M). 834* 835* ===================================================================== 836* 837* .. Parameters .. 838 DOUBLE PRECISION ONE, ZERO 839 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 840* .. 841* .. Local Scalars .. 842 INTEGER I, INFO, J 843 DOUBLE PRECISION CTEMP, STEMP, TEMP 844* .. 845* .. External Functions .. 846 LOGICAL LSAME 847 EXTERNAL LSAME 848* .. 849* .. External Subroutines .. 850 EXTERNAL XERBLA 851* .. 852* .. Intrinsic Functions .. 853 INTRINSIC MAX 854* .. 855* .. Executable Statements .. 856* 857* Test the input parameters 858* 859 INFO = 0 860 IF( .NOT.( LSAME( SIDE, 'L' ) .OR. LSAME( SIDE, 'R' ) ) ) THEN 861 INFO = 1 862 ELSE IF( .NOT.( LSAME( PIVOT, 'V' ) .OR. LSAME( PIVOT, 863 $ 'T' ) .OR. LSAME( PIVOT, 'B' ) ) ) THEN 864 INFO = 2 865 ELSE IF( .NOT.( LSAME( DIRECT, 'F' ) .OR. LSAME( DIRECT, 'B' ) ) ) 866 $ THEN 867 INFO = 3 868 ELSE IF( M.LT.0 ) THEN 869 INFO = 4 870 ELSE IF( N.LT.0 ) THEN 871 INFO = 5 872 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 873 INFO = 9 874 END IF 875 IF( INFO.NE.0 ) THEN 876 CALL XERBLA( 'DLASR ', INFO ) 877 RETURN 878 END IF 879* 880* Quick return if possible 881* 882 IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) 883 $ RETURN 884 IF( LSAME( SIDE, 'L' ) ) THEN 885* 886* Form P * A 887* 888 IF( LSAME( PIVOT, 'V' ) ) THEN 889 IF( LSAME( DIRECT, 'F' ) ) THEN 890 DO 20 J = 1, M - 1 891 CTEMP = C( J ) 892 STEMP = S( J ) 893 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 894 DO 10 I = 1, N 895 TEMP = A( J+1, I ) 896 A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I ) 897 A( J, I ) = STEMP*TEMP + CTEMP*A( J, I ) 898 10 CONTINUE 899 END IF 900 20 CONTINUE 901 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 902 DO 40 J = M - 1, 1, -1 903 CTEMP = C( J ) 904 STEMP = S( J ) 905 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 906 DO 30 I = 1, N 907 TEMP = A( J+1, I ) 908 A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I ) 909 A( J, I ) = STEMP*TEMP + CTEMP*A( J, I ) 910 30 CONTINUE 911 END IF 912 40 CONTINUE 913 END IF 914 ELSE IF( LSAME( PIVOT, 'T' ) ) THEN 915 IF( LSAME( DIRECT, 'F' ) ) THEN 916 DO 60 J = 2, M 917 CTEMP = C( J-1 ) 918 STEMP = S( J-1 ) 919 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 920 DO 50 I = 1, N 921 TEMP = A( J, I ) 922 A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I ) 923 A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I ) 924 50 CONTINUE 925 END IF 926 60 CONTINUE 927 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 928 DO 80 J = M, 2, -1 929 CTEMP = C( J-1 ) 930 STEMP = S( J-1 ) 931 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 932 DO 70 I = 1, N 933 TEMP = A( J, I ) 934 A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I ) 935 A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I ) 936 70 CONTINUE 937 END IF 938 80 CONTINUE 939 END IF 940 ELSE IF( LSAME( PIVOT, 'B' ) ) THEN 941 IF( LSAME( DIRECT, 'F' ) ) THEN 942 DO 100 J = 1, M - 1 943 CTEMP = C( J ) 944 STEMP = S( J ) 945 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 946 DO 90 I = 1, N 947 TEMP = A( J, I ) 948 A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP 949 A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP 950 90 CONTINUE 951 END IF 952 100 CONTINUE 953 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 954 DO 120 J = M - 1, 1, -1 955 CTEMP = C( J ) 956 STEMP = S( J ) 957 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 958 DO 110 I = 1, N 959 TEMP = A( J, I ) 960 A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP 961 A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP 962 110 CONTINUE 963 END IF 964 120 CONTINUE 965 END IF 966 END IF 967 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 968* 969* Form A * P' 970* 971 IF( LSAME( PIVOT, 'V' ) ) THEN 972 IF( LSAME( DIRECT, 'F' ) ) THEN 973 DO 140 J = 1, N - 1 974 CTEMP = C( J ) 975 STEMP = S( J ) 976 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 977 DO 130 I = 1, M 978 TEMP = A( I, J+1 ) 979 A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J ) 980 A( I, J ) = STEMP*TEMP + CTEMP*A( I, J ) 981 130 CONTINUE 982 END IF 983 140 CONTINUE 984 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 985 DO 160 J = N - 1, 1, -1 986 CTEMP = C( J ) 987 STEMP = S( J ) 988 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 989 DO 150 I = 1, M 990 TEMP = A( I, J+1 ) 991 A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J ) 992 A( I, J ) = STEMP*TEMP + CTEMP*A( I, J ) 993 150 CONTINUE 994 END IF 995 160 CONTINUE 996 END IF 997 ELSE IF( LSAME( PIVOT, 'T' ) ) THEN 998 IF( LSAME( DIRECT, 'F' ) ) THEN 999 DO 180 J = 2, N 1000 CTEMP = C( J-1 ) 1001 STEMP = S( J-1 ) 1002 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 1003 DO 170 I = 1, M 1004 TEMP = A( I, J ) 1005 A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 ) 1006 A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 ) 1007 170 CONTINUE 1008 END IF 1009 180 CONTINUE 1010 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 1011 DO 200 J = N, 2, -1 1012 CTEMP = C( J-1 ) 1013 STEMP = S( J-1 ) 1014 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 1015 DO 190 I = 1, M 1016 TEMP = A( I, J ) 1017 A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 ) 1018 A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 ) 1019 190 CONTINUE 1020 END IF 1021 200 CONTINUE 1022 END IF 1023 ELSE IF( LSAME( PIVOT, 'B' ) ) THEN 1024 IF( LSAME( DIRECT, 'F' ) ) THEN 1025 DO 220 J = 1, N - 1 1026 CTEMP = C( J ) 1027 STEMP = S( J ) 1028 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 1029 DO 210 I = 1, M 1030 TEMP = A( I, J ) 1031 A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP 1032 A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP 1033 210 CONTINUE 1034 END IF 1035 220 CONTINUE 1036 ELSE IF( LSAME( DIRECT, 'B' ) ) THEN 1037 DO 240 J = N - 1, 1, -1 1038 CTEMP = C( J ) 1039 STEMP = S( J ) 1040 IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN 1041 DO 230 I = 1, M 1042 TEMP = A( I, J ) 1043 A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP 1044 A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP 1045 230 CONTINUE 1046 END IF 1047 240 CONTINUE 1048 END IF 1049 END IF 1050 END IF 1051* 1052 RETURN 1053* 1054* End of DLASR 1055* 1056 END 1057 SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) 1058* 1059* -- LAPACK auxiliary routine (version 1.1) -- 1060* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 1061* Courant Institute, Argonne National Lab, and Rice University 1062* October 31, 1992 1063* 1064* .. Scalar Arguments .. 1065 DOUBLE PRECISION A, B, C, CS1, RT1, RT2, SN1 1066* .. 1067* 1068* Purpose 1069* ======= 1070* 1071* DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix 1072* [ A B ] 1073* [ B C ]. 1074* On return, RT1 is the eigenvalue of larger absolute value, RT2 is the 1075* eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right 1076* eigenvector for RT1, giving the decomposition 1077* 1078* [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] 1079* [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. 1080* 1081* Arguments 1082* ========= 1083* 1084* A (input) DOUBLE PRECISION 1085* The (1,1) entry of the 2-by-2 matrix. 1086* 1087* B (input) DOUBLE PRECISION 1088* The (1,2) entry and the conjugate of the (2,1) entry of the 1089* 2-by-2 matrix. 1090* 1091* C (input) DOUBLE PRECISION 1092* The (2,2) entry of the 2-by-2 matrix. 1093* 1094* RT1 (output) DOUBLE PRECISION 1095* The eigenvalue of larger absolute value. 1096* 1097* RT2 (output) DOUBLE PRECISION 1098* The eigenvalue of smaller absolute value. 1099* 1100* CS1 (output) DOUBLE PRECISION 1101* SN1 (output) DOUBLE PRECISION 1102* The vector (CS1, SN1) is a unit right eigenvector for RT1. 1103* 1104* Further Details 1105* =============== 1106* 1107* RT1 is accurate to a few ulps barring over/underflow. 1108* 1109* RT2 may be inaccurate if there is massive cancellation in the 1110* determinant A*C-B*B; higher precision or correctly rounded or 1111* correctly truncated arithmetic would be needed to compute RT2 1112* accurately in all cases. 1113* 1114* CS1 and SN1 are accurate to a few ulps barring over/underflow. 1115* 1116* Overflow is possible only if RT1 is within a factor of 5 of overflow. 1117* Underflow is harmless if the input data is 0 or exceeds 1118* underflow_threshold / macheps. 1119* 1120* ===================================================================== 1121* 1122* .. Parameters .. 1123 DOUBLE PRECISION ONE 1124 PARAMETER ( ONE = 1.0D0 ) 1125 DOUBLE PRECISION TWO 1126 PARAMETER ( TWO = 2.0D0 ) 1127 DOUBLE PRECISION ZERO 1128 PARAMETER ( ZERO = 0.0D0 ) 1129 DOUBLE PRECISION HALF 1130 PARAMETER ( HALF = 0.5D0 ) 1131* .. 1132* .. Local Scalars .. 1133 INTEGER SGN1, SGN2 1134 DOUBLE PRECISION AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM, 1135 $ TB, TN 1136* .. 1137* .. Intrinsic Functions .. 1138 INTRINSIC ABS, SQRT 1139* .. 1140* .. Executable Statements .. 1141* 1142* Compute the eigenvalues 1143* 1144 SM = A + C 1145 DF = A - C 1146 ADF = ABS( DF ) 1147 TB = B + B 1148 AB = ABS( TB ) 1149 IF( ABS( A ).GT.ABS( C ) ) THEN 1150 ACMX = A 1151 ACMN = C 1152 ELSE 1153 ACMX = C 1154 ACMN = A 1155 END IF 1156 IF( ADF.GT.AB ) THEN 1157 RT = ADF*SQRT( ONE+( AB / ADF )**2 ) 1158 ELSE IF( ADF.LT.AB ) THEN 1159 RT = AB*SQRT( ONE+( ADF / AB )**2 ) 1160 ELSE 1161* 1162* Includes case AB=ADF=0 1163* 1164 RT = AB*SQRT( TWO ) 1165 END IF 1166 IF( SM.LT.ZERO ) THEN 1167 RT1 = HALF*( SM-RT ) 1168 SGN1 = -1 1169* 1170* Order of execution important. 1171* To get fully accurate smaller eigenvalue, 1172* next line needs to be executed in higher precision. 1173* 1174 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 1175 ELSE IF( SM.GT.ZERO ) THEN 1176 RT1 = HALF*( SM+RT ) 1177 SGN1 = 1 1178* 1179* Order of execution important. 1180* To get fully accurate smaller eigenvalue, 1181* next line needs to be executed in higher precision. 1182* 1183 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 1184 ELSE 1185* 1186* Includes case RT1 = RT2 = 0 1187* 1188 RT1 = HALF*RT 1189 RT2 = -HALF*RT 1190 SGN1 = 1 1191 END IF 1192* 1193* Compute the eigenvector 1194* 1195 IF( DF.GE.ZERO ) THEN 1196 CS = DF + RT 1197 SGN2 = 1 1198 ELSE 1199 CS = DF - RT 1200 SGN2 = -1 1201 END IF 1202 ACS = ABS( CS ) 1203 IF( ACS.GT.AB ) THEN 1204 CT = -TB / CS 1205 SN1 = ONE / SQRT( ONE+CT*CT ) 1206 CS1 = CT*SN1 1207 ELSE 1208 IF( AB.EQ.ZERO ) THEN 1209 CS1 = ONE 1210 SN1 = ZERO 1211 ELSE 1212 TN = -CS / TB 1213 CS1 = ONE / SQRT( ONE+TN*TN ) 1214 SN1 = TN*CS1 1215 END IF 1216 END IF 1217 IF( SGN1.EQ.SGN2 ) THEN 1218 TN = CS1 1219 CS1 = -SN1 1220 SN1 = TN 1221 END IF 1222 RETURN 1223* 1224* End of DLAEV2 1225* 1226 END 1227 SUBROUTINE DLAZRO( M, N, ALPHA, BETA, A, LDA ) 1228* 1229* -- LAPACK auxiliary routine (version 1.1) -- 1230* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 1231* Courant Institute, Argonne National Lab, and Rice University 1232* October 31, 1992 1233* 1234* .. Scalar Arguments .. 1235 INTEGER LDA, M, N 1236 DOUBLE PRECISION ALPHA, BETA 1237* .. 1238* .. Array Arguments .. 1239 DOUBLE PRECISION A( LDA, * ) 1240* .. 1241* 1242* Purpose 1243* ======= 1244* 1245* DLAZRO initializes a 2-D array A to BETA on the diagonal and 1246* ALPHA on the offdiagonals. 1247* 1248* Arguments 1249* ========= 1250* 1251* M (input) INTEGER 1252* The number of rows of the matrix A. M >= 0. 1253* 1254* N (input) INTEGER 1255* The number of columns of the matrix A. N >= 0. 1256* 1257* ALPHA (input) DOUBLE PRECISION 1258* The constant to which the offdiagonal elements are to be set. 1259* 1260* BETA (input) DOUBLE PRECISION 1261* The constant to which the diagonal elements are to be set. 1262* 1263* A (output) DOUBLE PRECISION array, dimension (LDA,N) 1264* On exit, the leading m by n submatrix of A is set such that 1265* A(i,j) = ALPHA, 1 <= i <= m, 1 <= j <= n, i <> j 1266* A(i,i) = BETA, 1 <= i <= min(m,n). 1267* 1268* LDA (input) INTEGER 1269* The leading dimension of the array A. LDA >= max(1,M). 1270* 1271* ===================================================================== 1272* 1273* .. Local Scalars .. 1274 INTEGER I, J 1275* .. 1276* .. Intrinsic Functions .. 1277 INTRINSIC MIN 1278* .. 1279* .. Executable Statements .. 1280* 1281 DO 20 J = 1, N 1282 DO 10 I = 1, M 1283 A( I, J ) = ALPHA 1284 10 CONTINUE 1285 20 CONTINUE 1286* 1287 DO 30 I = 1, MIN( M, N ) 1288 A( I, I ) = BETA 1289 30 CONTINUE 1290* 1291 RETURN 1292* 1293* End of DLAZRO 1294* 1295 END 1296 SUBROUTINE DORGTR( UPLO, N, A, LDA, TAU, WORK, LWORK, INFO ) 1297* 1298* -- LAPACK routine (version 1.1) -- 1299* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 1300* Courant Institute, Argonne National Lab, and Rice University 1301* March 31, 1993 1302* 1303* .. Scalar Arguments .. 1304 CHARACTER UPLO 1305 INTEGER INFO, LDA, LWORK, N 1306* .. 1307* .. Array Arguments .. 1308 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( LWORK ) 1309* .. 1310* 1311* Purpose 1312* ======= 1313* 1314* DORGTR generates a real orthogonal matrix Q which is defined as the 1315* product of n-1 elementary reflectors of order N, as returned by 1316* DSYTRD: 1317* 1318* if UPLO = 'U', Q = H(n-1) . . . H(2) H(1), 1319* 1320* if UPLO = 'L', Q = H(1) H(2) . . . H(n-1). 1321* 1322* Arguments 1323* ========= 1324* 1325* UPLO (input) CHARACTER*1 1326* = 'U': Upper triangle of A contains elementary reflectors 1327* from DSYTRD; 1328* = 'L': Lower triangle of A contains elementary reflectors 1329* from DSYTRD. 1330* 1331* N (input) INTEGER 1332* The order of the matrix Q. N >= 0. 1333* 1334* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 1335* On entry, the vectors which define the elementary reflectors, 1336* as returned by DSYTRD. 1337* On exit, the N-by-N orthogonal matrix Q. 1338* 1339* LDA (input) INTEGER 1340* The leading dimension of the array A. LDA >= max(1,N). 1341* 1342* TAU (input) DOUBLE PRECISION array, dimension (N-1) 1343* TAU(i) must contain the scalar factor of the elementary 1344* reflector H(i), as returned by DSYTRD. 1345* 1346* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 1347* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 1348* 1349* LWORK (input) INTEGER 1350* The dimension of the array WORK. LWORK >= max(1,N-1). 1351* For optimum performance LWORK >= (N-1)*NB, where NB is 1352* the optimal blocksize. 1353* 1354* INFO (output) INTEGER 1355* = 0: successful exit 1356* < 0: if INFO = -i, the i-th argument had an illegal value 1357* 1358* ===================================================================== 1359* 1360* .. Parameters .. 1361 DOUBLE PRECISION ZERO, ONE 1362 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 1363* .. 1364* .. Local Scalars .. 1365 LOGICAL UPPER 1366 INTEGER I, IINFO, J 1367* .. 1368* .. External Functions .. 1369 LOGICAL LSAME 1370 EXTERNAL LSAME 1371* .. 1372* .. External Subroutines .. 1373 EXTERNAL DORGQL, DORGQR, XERBLA 1374* .. 1375* .. Intrinsic Functions .. 1376 INTRINSIC MAX 1377* .. 1378* .. Executable Statements .. 1379* 1380* Test the input arguments 1381* 1382 INFO = 0 1383 UPPER = LSAME( UPLO, 'U' ) 1384 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 1385 INFO = -1 1386 ELSE IF( N.LT.0 ) THEN 1387 INFO = -2 1388 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 1389 INFO = -4 1390 ELSE IF( LWORK.LT.MAX( 1, N-1 ) ) THEN 1391 INFO = -7 1392 END IF 1393 IF( INFO.NE.0 ) THEN 1394 CALL XERBLA( 'DORGTR', -INFO ) 1395 RETURN 1396 END IF 1397* 1398* Quick return if possible 1399* 1400 IF( N.EQ.0 ) THEN 1401 WORK( 1 ) = 1 1402 RETURN 1403 END IF 1404* 1405 IF( UPPER ) THEN 1406* 1407* Q was determined by a call to DSYTRD with UPLO = 'U' 1408* 1409* Shift the vectors which define the elementary reflectors one 1410* column to the left, and set the last row and column of Q to 1411* those of the unit matrix 1412* 1413 DO 20 J = 1, N - 1 1414 DO 10 I = 1, J - 1 1415 A( I, J ) = A( I, J+1 ) 1416 10 CONTINUE 1417 A( N, J ) = ZERO 1418 20 CONTINUE 1419 DO 30 I = 1, N - 1 1420 A( I, N ) = ZERO 1421 30 CONTINUE 1422 A( N, N ) = ONE 1423* 1424* Generate Q(1:n-1,1:n-1) 1425* 1426 CALL DORGQL( N-1, N-1, N-1, A, LDA, TAU, WORK, LWORK, IINFO ) 1427* 1428 ELSE 1429* 1430* Q was determined by a call to DSYTRD with UPLO = 'L'. 1431* 1432* Shift the vectors which define the elementary reflectors one 1433* column to the right, and set the first row and column of Q to 1434* those of the unit matrix 1435* 1436 DO 50 J = N, 2, -1 1437 A( 1, J ) = ZERO 1438 DO 40 I = J + 1, N 1439 A( I, J ) = A( I, J-1 ) 1440 40 CONTINUE 1441 50 CONTINUE 1442 A( 1, 1 ) = ONE 1443 DO 60 I = 2, N 1444 A( I, 1 ) = ZERO 1445 60 CONTINUE 1446 IF( N.GT.1 ) THEN 1447* 1448* Generate Q(2:n,2:n) 1449* 1450 CALL DORGQR( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK, 1451 $ LWORK, IINFO ) 1452 END IF 1453 END IF 1454 RETURN 1455* 1456* End of DORGTR 1457* 1458 END 1459 SUBROUTINE DORGQR( M, N, K, A, LDA, TAU, WORK, LWORK, INFO ) 1460* 1461* -- LAPACK routine (version 1.1) -- 1462* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 1463* Courant Institute, Argonne National Lab, and Rice University 1464* March 31, 1993 1465* 1466* .. Scalar Arguments .. 1467 INTEGER INFO, K, LDA, LWORK, M, N 1468* .. 1469* .. Array Arguments .. 1470 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( LWORK ) 1471* .. 1472* 1473* Purpose 1474* ======= 1475* 1476* DORGQR generates an M-by-N real matrix Q with orthonormal columns, 1477* which is defined as the first N columns of a product of K elementary 1478* reflectors of order M 1479* 1480* Q = H(1) H(2) . . . H(k) 1481* 1482* as returned by DGEQRF. 1483* 1484* Arguments 1485* ========= 1486* 1487* M (input) INTEGER 1488* The number of rows of the matrix Q. M >= 0. 1489* 1490* N (input) INTEGER 1491* The number of columns of the matrix Q. M >= N >= 0. 1492* 1493* K (input) INTEGER 1494* The number of elementary reflectors whose product defines the 1495* matrix Q. N >= K >= 0. 1496* 1497* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 1498* On entry, the i-th column must contain the vector which 1499* defines the elementary reflector H(i), for i = 1,2,...,k, as 1500* returned by DGEQRF in the first k columns of its array 1501* argument A. 1502* On exit, the M-by-N matrix Q. 1503* 1504* LDA (input) INTEGER 1505* The first dimension of the array A. LDA >= max(1,M). 1506* 1507* TAU (input) DOUBLE PRECISION array, dimension (K) 1508* TAU(i) must contain the scalar factor of the elementary 1509* reflector H(i), as returned by DGEQRF. 1510* 1511* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 1512* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 1513* 1514* LWORK (input) INTEGER 1515* The dimension of the array WORK. LWORK >= max(1,N). 1516* For optimum performance LWORK >= N*NB, where NB is the 1517* optimal blocksize. 1518* 1519* INFO (output) INTEGER 1520* = 0: successful exit 1521* < 0: if INFO = -i, the i-th argument has an illegal value 1522* 1523* ===================================================================== 1524* 1525* .. Parameters .. 1526 DOUBLE PRECISION ZERO 1527 PARAMETER ( ZERO = 0.0D+0 ) 1528* .. 1529* .. Local Scalars .. 1530 INTEGER I, IB, IINFO, IWS, J, KI, KK, L, LDWORK, NB, 1531 $ NBMIN, NX 1532* .. 1533* .. External Subroutines .. 1534 EXTERNAL DLARFB, DLARFT, DORG2R, XERBLA 1535* .. 1536* .. Intrinsic Functions .. 1537 INTRINSIC MAX, MIN 1538* .. 1539* .. External Functions .. 1540 INTEGER ILAENV 1541 EXTERNAL ILAENV 1542* .. 1543* .. Executable Statements .. 1544* 1545* Test the input arguments 1546* 1547 INFO = 0 1548 IF( M.LT.0 ) THEN 1549 INFO = -1 1550 ELSE IF( N.LT.0 .OR. N.GT.M ) THEN 1551 INFO = -2 1552 ELSE IF( K.LT.0 .OR. K.GT.N ) THEN 1553 INFO = -3 1554 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 1555 INFO = -5 1556 ELSE IF( LWORK.LT.MAX( 1, N ) ) THEN 1557 INFO = -8 1558 END IF 1559 IF( INFO.NE.0 ) THEN 1560 CALL XERBLA( 'DORGQR', -INFO ) 1561 RETURN 1562 END IF 1563* 1564* Quick return if possible 1565* 1566 IF( N.LE.0 ) THEN 1567 WORK( 1 ) = 1 1568 RETURN 1569 END IF 1570* 1571* Determine the block size. 1572* 1573 NB = ILAENV( 1, 'DORGQR', ' ', M, N, K, -1 ) 1574 NBMIN = 2 1575 NX = 0 1576 IWS = N 1577 IF( NB.GT.1 .AND. NB.LT.K ) THEN 1578* 1579* Determine when to cross over from blocked to unblocked code. 1580* 1581 NX = MAX( 0, ILAENV( 3, 'DORGQR', ' ', M, N, K, -1 ) ) 1582 IF( NX.LT.K ) THEN 1583* 1584* Determine if workspace is large enough for blocked code. 1585* 1586 LDWORK = N 1587 IWS = LDWORK*NB 1588 IF( LWORK.LT.IWS ) THEN 1589* 1590* Not enough workspace to use optimal NB: reduce NB and 1591* determine the minimum value of NB. 1592* 1593 NB = LWORK / LDWORK 1594 NBMIN = MAX( 2, ILAENV( 2, 'DORGQR', ' ', M, N, K, -1 ) ) 1595 END IF 1596 END IF 1597 END IF 1598* 1599 IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN 1600* 1601* Use blocked code after the last block. 1602* The first kk columns are handled by the block method. 1603* 1604 KI = ( ( K-NX-1 ) / NB )*NB 1605 KK = MIN( K, KI+NB ) 1606* 1607* Set A(1:kk,kk+1:n) to zero. 1608* 1609 DO 20 J = KK + 1, N 1610 DO 10 I = 1, KK 1611 A( I, J ) = ZERO 1612 10 CONTINUE 1613 20 CONTINUE 1614 ELSE 1615 KK = 0 1616 END IF 1617* 1618* Use unblocked code for the last or only block. 1619* 1620 IF( KK.LT.N ) 1621 $ CALL DORG2R( M-KK, N-KK, K-KK, A( KK+1, KK+1 ), LDA, 1622 $ TAU( KK+1 ), WORK, IINFO ) 1623* 1624 IF( KK.GT.0 ) THEN 1625* 1626* Use blocked code 1627* 1628 DO 50 I = KI + 1, 1, -NB 1629 IB = MIN( NB, K-I+1 ) 1630 IF( I+IB.LE.N ) THEN 1631* 1632* Form the triangular factor of the block reflector 1633* H = H(i) H(i+1) . . . H(i+ib-1) 1634* 1635 CALL DLARFT( 'Forward', 'Columnwise', M-I+1, IB, 1636 $ A( I, I ), LDA, TAU( I ), WORK, LDWORK ) 1637* 1638* Apply H to A(i:m,i+ib:n) from the left 1639* 1640 CALL DLARFB( 'Left', 'No transpose', 'Forward', 1641 $ 'Columnwise', M-I+1, N-I-IB+1, IB, 1642 $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ), 1643 $ LDA, WORK( IB+1 ), LDWORK ) 1644 END IF 1645* 1646* Apply H to rows i:m of current block 1647* 1648 CALL DORG2R( M-I+1, IB, IB, A( I, I ), LDA, TAU( I ), WORK, 1649 $ IINFO ) 1650* 1651* Set rows 1:i-1 of current block to zero 1652* 1653 DO 40 J = I, I + IB - 1 1654 DO 30 L = 1, I - 1 1655 A( L, J ) = ZERO 1656 30 CONTINUE 1657 40 CONTINUE 1658 50 CONTINUE 1659 END IF 1660* 1661 WORK( 1 ) = IWS 1662 RETURN 1663* 1664* End of DORGQR 1665* 1666 END 1667 SUBROUTINE DORG2R( M, N, K, A, LDA, TAU, WORK, INFO ) 1668* 1669* -- LAPACK routine (version 1.1) -- 1670* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 1671* Courant Institute, Argonne National Lab, and Rice University 1672* February 29, 1992 1673* 1674* .. Scalar Arguments .. 1675 INTEGER INFO, K, LDA, M, N 1676* .. 1677* .. Array Arguments .. 1678 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) 1679* .. 1680* 1681* Purpose 1682* ======= 1683* 1684* DORG2R generates an m by n real matrix Q with orthonormal columns, 1685* which is defined as the first n columns of a product of k elementary 1686* reflectors of order m 1687* 1688* Q = H(1) H(2) . . . H(k) 1689* 1690* as returned by DGEQRF. 1691* 1692* Arguments 1693* ========= 1694* 1695* M (input) INTEGER 1696* The number of rows of the matrix Q. M >= 0. 1697* 1698* N (input) INTEGER 1699* The number of columns of the matrix Q. M >= N >= 0. 1700* 1701* K (input) INTEGER 1702* The number of elementary reflectors whose product defines the 1703* matrix Q. N >= K >= 0. 1704* 1705* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 1706* On entry, the i-th column must contain the vector which 1707* defines the elementary reflector H(i), for i = 1,2,...,k, as 1708* returned by DGEQRF in the first k columns of its array 1709* argument A. 1710* On exit, the m-by-n matrix Q. 1711* 1712* LDA (input) INTEGER 1713* The first dimension of the array A. LDA >= max(1,M). 1714* 1715* TAU (input) DOUBLE PRECISION array, dimension (K) 1716* TAU(i) must contain the scalar factor of the elementary 1717* reflector H(i), as returned by DGEQRF. 1718* 1719* WORK (workspace) DOUBLE PRECISION array, dimension (N) 1720* 1721* INFO (output) INTEGER 1722* = 0: successful exit 1723* < 0: if INFO = -i, the i-th argument has an illegal value 1724* 1725* ===================================================================== 1726* 1727* .. Parameters .. 1728 DOUBLE PRECISION ONE, ZERO 1729 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 1730* .. 1731* .. Local Scalars .. 1732 INTEGER I, J, L 1733* .. 1734* .. External Subroutines .. 1735 EXTERNAL DLARF, DSCAL, XERBLA 1736* .. 1737* .. Intrinsic Functions .. 1738 INTRINSIC MAX 1739* .. 1740* .. Executable Statements .. 1741* 1742* Test the input arguments 1743* 1744 INFO = 0 1745 IF( M.LT.0 ) THEN 1746 INFO = -1 1747 ELSE IF( N.LT.0 .OR. N.GT.M ) THEN 1748 INFO = -2 1749 ELSE IF( K.LT.0 .OR. K.GT.N ) THEN 1750 INFO = -3 1751 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 1752 INFO = -5 1753 END IF 1754 IF( INFO.NE.0 ) THEN 1755 CALL XERBLA( 'DORG2R', -INFO ) 1756 RETURN 1757 END IF 1758* 1759* Quick return if possible 1760* 1761 IF( N.LE.0 ) 1762 $ RETURN 1763* 1764* Initialise columns k+1:n to columns of the unit matrix 1765* 1766 DO 20 J = K + 1, N 1767 DO 10 L = 1, M 1768 A( L, J ) = ZERO 1769 10 CONTINUE 1770 A( J, J ) = ONE 1771 20 CONTINUE 1772* 1773 DO 40 I = K, 1, -1 1774* 1775* Apply H(i) to A(i:m,i:n) from the left 1776* 1777 IF( I.LT.N ) THEN 1778 A( I, I ) = ONE 1779 CALL DLARF( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), 1780 $ A( I, I+1 ), LDA, WORK ) 1781 END IF 1782 IF( I.LT.M ) 1783 $ CALL DSCAL( M-I, -TAU( I ), A( I+1, I ), 1 ) 1784 A( I, I ) = ONE - TAU( I ) 1785* 1786* Set A(1:i-1,i) to zero 1787* 1788 DO 30 L = 1, I - 1 1789 A( L, I ) = ZERO 1790 30 CONTINUE 1791 40 CONTINUE 1792 RETURN 1793* 1794* End of DORG2R 1795* 1796 END 1797 SUBROUTINE DORGQL( M, N, K, A, LDA, TAU, WORK, LWORK, INFO ) 1798* 1799* -- LAPACK routine (version 1.1) -- 1800* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 1801* Courant Institute, Argonne National Lab, and Rice University 1802* March 31, 1993 1803* 1804* .. Scalar Arguments .. 1805 INTEGER INFO, K, LDA, LWORK, M, N 1806* .. 1807* .. Array Arguments .. 1808 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( LWORK ) 1809* .. 1810* 1811* Purpose 1812* ======= 1813* 1814* DORGQL generates an M-by-N real matrix Q with orthonormal columns, 1815* which is defined as the last N columns of a product of K elementary 1816* reflectors of order M 1817* 1818* Q = H(k) . . . H(2) H(1) 1819* 1820* as returned by DGEQLF. 1821* 1822* Arguments 1823* ========= 1824* 1825* M (input) INTEGER 1826* The number of rows of the matrix Q. M >= 0. 1827* 1828* N (input) INTEGER 1829* The number of columns of the matrix Q. M >= N >= 0. 1830* 1831* K (input) INTEGER 1832* The number of elementary reflectors whose product defines the 1833* matrix Q. N >= K >= 0. 1834* 1835* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 1836* On entry, the (n-k+i)-th column must contain the vector which 1837* defines the elementary reflector H(i), for i = 1,2,...,k, as 1838* returned by DGEQLF in the last k columns of its array 1839* argument A. 1840* On exit, the M-by-N matrix Q. 1841* 1842* LDA (input) INTEGER 1843* The first dimension of the array A. LDA >= max(1,M). 1844* 1845* TAU (input) DOUBLE PRECISION array, dimension (K) 1846* TAU(i) must contain the scalar factor of the elementary 1847* reflector H(i), as returned by DGEQLF. 1848* 1849* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 1850* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 1851* 1852* LWORK (input) INTEGER 1853* The dimension of the array WORK. LWORK >= max(1,N). 1854* For optimum performance LWORK >= N*NB, where NB is the 1855* optimal blocksize. 1856* 1857* INFO (output) INTEGER 1858* = 0: successful exit 1859* < 0: if INFO = -i, the i-th argument has an illegal value 1860* 1861* ===================================================================== 1862* 1863* .. Parameters .. 1864 DOUBLE PRECISION ZERO 1865 PARAMETER ( ZERO = 0.0D+0 ) 1866* .. 1867* .. Local Scalars .. 1868 INTEGER I, IB, IINFO, IWS, J, KK, L, LDWORK, NB, NBMIN, 1869 $ NX 1870* .. 1871* .. External Subroutines .. 1872 EXTERNAL DLARFB, DLARFT, DORG2L, XERBLA 1873* .. 1874* .. Intrinsic Functions .. 1875 INTRINSIC MAX, MIN 1876* .. 1877* .. External Functions .. 1878 INTEGER ILAENV 1879 EXTERNAL ILAENV 1880* .. 1881* .. Executable Statements .. 1882* 1883* Test the input arguments 1884* 1885 INFO = 0 1886 IF( M.LT.0 ) THEN 1887 INFO = -1 1888 ELSE IF( N.LT.0 .OR. N.GT.M ) THEN 1889 INFO = -2 1890 ELSE IF( K.LT.0 .OR. K.GT.N ) THEN 1891 INFO = -3 1892 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 1893 INFO = -5 1894 ELSE IF( LWORK.LT.MAX( 1, N ) ) THEN 1895 INFO = -8 1896 END IF 1897 IF( INFO.NE.0 ) THEN 1898 CALL XERBLA( 'DORGQL', -INFO ) 1899 RETURN 1900 END IF 1901* 1902* Quick return if possible 1903* 1904 IF( N.LE.0 ) THEN 1905 WORK( 1 ) = 1 1906 RETURN 1907 END IF 1908* 1909* Determine the block size. 1910* 1911 NB = ILAENV( 1, 'DORGQL', ' ', M, N, K, -1 ) 1912 NBMIN = 2 1913 NX = 0 1914 IWS = N 1915 IF( NB.GT.1 .AND. NB.LT.K ) THEN 1916* 1917* Determine when to cross over from blocked to unblocked code. 1918* 1919 NX = MAX( 0, ILAENV( 3, 'DORGQL', ' ', M, N, K, -1 ) ) 1920 IF( NX.LT.K ) THEN 1921* 1922* Determine if workspace is large enough for blocked code. 1923* 1924 LDWORK = N 1925 IWS = LDWORK*NB 1926 IF( LWORK.LT.IWS ) THEN 1927* 1928* Not enough workspace to use optimal NB: reduce NB and 1929* determine the minimum value of NB. 1930* 1931 NB = LWORK / LDWORK 1932 NBMIN = MAX( 2, ILAENV( 2, 'DORGQL', ' ', M, N, K, -1 ) ) 1933 END IF 1934 END IF 1935 END IF 1936* 1937 IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN 1938* 1939* Use blocked code after the first block. 1940* The last kk columns are handled by the block method. 1941* 1942 KK = MIN( K, ( ( K-NX+NB-1 ) / NB )*NB ) 1943* 1944* Set A(m-kk+1:m,1:n-kk) to zero. 1945* 1946 DO 20 J = 1, N - KK 1947 DO 10 I = M - KK + 1, M 1948 A( I, J ) = ZERO 1949 10 CONTINUE 1950 20 CONTINUE 1951 ELSE 1952 KK = 0 1953 END IF 1954* 1955* Use unblocked code for the first or only block. 1956* 1957 CALL DORG2L( M-KK, N-KK, K-KK, A, LDA, TAU, WORK, IINFO ) 1958* 1959 IF( KK.GT.0 ) THEN 1960* 1961* Use blocked code 1962* 1963 DO 50 I = K - KK + 1, K, NB 1964 IB = MIN( NB, K-I+1 ) 1965 IF( N-K+I.GT.1 ) THEN 1966* 1967* Form the triangular factor of the block reflector 1968* H = H(i+ib-1) . . . H(i+1) H(i) 1969* 1970 CALL DLARFT( 'Backward', 'Columnwise', M-K+I+IB-1, IB, 1971 $ A( 1, N-K+I ), LDA, TAU( I ), WORK, LDWORK ) 1972* 1973* Apply H to A(1:m-k+i+ib-1,1:n-k+i-1) from the left 1974* 1975 CALL DLARFB( 'Left', 'No transpose', 'Backward', 1976 $ 'Columnwise', M-K+I+IB-1, N-K+I-1, IB, 1977 $ A( 1, N-K+I ), LDA, WORK, LDWORK, A, LDA, 1978 $ WORK( IB+1 ), LDWORK ) 1979 END IF 1980* 1981* Apply H to rows 1:m-k+i+ib-1 of current block 1982* 1983 CALL DORG2L( M-K+I+IB-1, IB, IB, A( 1, N-K+I ), LDA, 1984 $ TAU( I ), WORK, IINFO ) 1985* 1986* Set rows m-k+i+ib:m of current block to zero 1987* 1988 DO 40 J = N - K + I, N - K + I + IB - 1 1989 DO 30 L = M - K + I + IB, M 1990 A( L, J ) = ZERO 1991 30 CONTINUE 1992 40 CONTINUE 1993 50 CONTINUE 1994 END IF 1995* 1996 WORK( 1 ) = IWS 1997 RETURN 1998* 1999* End of DORGQL 2000* 2001 END 2002 SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, 2003 $ T, LDT, C, LDC, WORK, LDWORK ) 2004* 2005* -- LAPACK auxiliary routine (version 1.1) -- 2006* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 2007* Courant Institute, Argonne National Lab, and Rice University 2008* February 29, 1992 2009* 2010* .. Scalar Arguments .. 2011 CHARACTER DIRECT, SIDE, STOREV, TRANS 2012 INTEGER K, LDC, LDT, LDV, LDWORK, M, N 2013* .. 2014* .. Array Arguments .. 2015 DOUBLE PRECISION C( LDC, * ), T( LDT, * ), V( LDV, * ), 2016 $ WORK( LDWORK, * ) 2017* .. 2018* 2019* Purpose 2020* ======= 2021* 2022* DLARFB applies a real block reflector H or its transpose H' to a 2023* real m by n matrix C, from either the left or the right. 2024* 2025* Arguments 2026* ========= 2027* 2028* SIDE (input) CHARACTER*1 2029* = 'L': apply H or H' from the Left 2030* = 'R': apply H or H' from the Right 2031* 2032* TRANS (input) CHARACTER*1 2033* = 'N': apply H (No transpose) 2034* = 'T': apply H' (Transpose) 2035* 2036* DIRECT (input) CHARACTER*1 2037* Indicates how H is formed from a product of elementary 2038* reflectors 2039* = 'F': H = H(1) H(2) . . . H(k) (Forward) 2040* = 'B': H = H(k) . . . H(2) H(1) (Backward) 2041* 2042* STOREV (input) CHARACTER*1 2043* Indicates how the vectors which define the elementary 2044* reflectors are stored: 2045* = 'C': Columnwise 2046* = 'R': Rowwise 2047* 2048* M (input) INTEGER 2049* The number of rows of the matrix C. 2050* 2051* N (input) INTEGER 2052* The number of columns of the matrix C. 2053* 2054* K (input) INTEGER 2055* The order of the matrix T (= the number of elementary 2056* reflectors whose product defines the block reflector). 2057* 2058* V (input) DOUBLE PRECISION array, dimension 2059* (LDV,K) if STOREV = 'C' 2060* (LDV,M) if STOREV = 'R' and SIDE = 'L' 2061* (LDV,N) if STOREV = 'R' and SIDE = 'R' 2062* The matrix V. See further details. 2063* 2064* LDV (input) INTEGER 2065* The leading dimension of the array V. 2066* If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); 2067* if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); 2068* if STOREV = 'R', LDV >= K. 2069* 2070* T (input) DOUBLE PRECISION array, dimension (LDT,K) 2071* The triangular k by k matrix T in the representation of the 2072* block reflector. 2073* 2074* LDT (input) INTEGER 2075* The leading dimension of the array T. LDT >= K. 2076* 2077* C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 2078* On entry, the m by n matrix C. 2079* On exit, C is overwritten by H*C or H'*C or C*H or C*H'. 2080* 2081* LDC (input) INTEGER 2082* The leading dimension of the array C. LDA >= max(1,M). 2083* 2084* WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,K) 2085* 2086* LDWORK (input) INTEGER 2087* The leading dimension of the array WORK. 2088* If SIDE = 'L', LDWORK >= max(1,N); 2089* if SIDE = 'R', LDWORK >= max(1,M). 2090* 2091* ===================================================================== 2092* 2093* .. Parameters .. 2094 DOUBLE PRECISION ONE 2095 PARAMETER ( ONE = 1.0D+0 ) 2096* .. 2097* .. Local Scalars .. 2098 CHARACTER TRANST 2099 INTEGER I, J 2100* .. 2101* .. External Functions .. 2102 LOGICAL LSAME 2103 EXTERNAL LSAME 2104* .. 2105* .. External Subroutines .. 2106 EXTERNAL DCOPY, DGEMM, DTRMM 2107* .. 2108* .. Executable Statements .. 2109* 2110* Quick return if possible 2111* 2112 IF( M.LE.0 .OR. N.LE.0 ) 2113 $ RETURN 2114* 2115 IF( LSAME( TRANS, 'N' ) ) THEN 2116 TRANST = 'T' 2117 ELSE 2118 TRANST = 'N' 2119 END IF 2120* 2121 IF( LSAME( STOREV, 'C' ) ) THEN 2122* 2123 IF( LSAME( DIRECT, 'F' ) ) THEN 2124* 2125* Let V = ( V1 ) (first K rows) 2126* ( V2 ) 2127* where V1 is unit lower triangular. 2128* 2129 IF( LSAME( SIDE, 'L' ) ) THEN 2130* 2131* Form H * C or H' * C where C = ( C1 ) 2132* ( C2 ) 2133* 2134* W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) 2135* 2136* W := C1' 2137* 2138 DO 10 J = 1, K 2139 CALL DCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 2140 10 CONTINUE 2141* 2142* W := W * V1 2143* 2144 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N, 2145 $ K, ONE, V, LDV, WORK, LDWORK ) 2146 IF( M.GT.K ) THEN 2147* 2148* W := W + C2'*V2 2149* 2150 CALL DGEMM( 'Transpose', 'No transpose', N, K, M-K, 2151 $ ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV, 2152 $ ONE, WORK, LDWORK ) 2153 END IF 2154* 2155* W := W * T' or W * T 2156* 2157 CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K, 2158 $ ONE, T, LDT, WORK, LDWORK ) 2159* 2160* C := C - V * W' 2161* 2162 IF( M.GT.K ) THEN 2163* 2164* C2 := C2 - V2 * W' 2165* 2166 CALL DGEMM( 'No transpose', 'Transpose', M-K, N, K, 2167 $ -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE, 2168 $ C( K+1, 1 ), LDC ) 2169 END IF 2170* 2171* W := W * V1' 2172* 2173 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', N, K, 2174 $ ONE, V, LDV, WORK, LDWORK ) 2175* 2176* C1 := C1 - W' 2177* 2178 DO 30 J = 1, K 2179 DO 20 I = 1, N 2180 C( J, I ) = C( J, I ) - WORK( I, J ) 2181 20 CONTINUE 2182 30 CONTINUE 2183* 2184 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 2185* 2186* Form C * H or C * H' where C = ( C1 C2 ) 2187* 2188* W := C * V = (C1*V1 + C2*V2) (stored in WORK) 2189* 2190* W := C1 2191* 2192 DO 40 J = 1, K 2193 CALL DCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 ) 2194 40 CONTINUE 2195* 2196* W := W * V1 2197* 2198 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M, 2199 $ K, ONE, V, LDV, WORK, LDWORK ) 2200 IF( N.GT.K ) THEN 2201* 2202* W := W + C2 * V2 2203* 2204 CALL DGEMM( 'No transpose', 'No transpose', M, K, N-K, 2205 $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV, 2206 $ ONE, WORK, LDWORK ) 2207 END IF 2208* 2209* W := W * T or W * T' 2210* 2211 CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K, 2212 $ ONE, T, LDT, WORK, LDWORK ) 2213* 2214* C := C - W * V' 2215* 2216 IF( N.GT.K ) THEN 2217* 2218* C2 := C2 - W * V2' 2219* 2220 CALL DGEMM( 'No transpose', 'Transpose', M, N-K, K, 2221 $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE, 2222 $ C( 1, K+1 ), LDC ) 2223 END IF 2224* 2225* W := W * V1' 2226* 2227 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', M, K, 2228 $ ONE, V, LDV, WORK, LDWORK ) 2229* 2230* C1 := C1 - W 2231* 2232 DO 60 J = 1, K 2233 DO 50 I = 1, M 2234 C( I, J ) = C( I, J ) - WORK( I, J ) 2235 50 CONTINUE 2236 60 CONTINUE 2237 END IF 2238* 2239 ELSE 2240* 2241* Let V = ( V1 ) 2242* ( V2 ) (last K rows) 2243* where V2 is unit upper triangular. 2244* 2245 IF( LSAME( SIDE, 'L' ) ) THEN 2246* 2247* Form H * C or H' * C where C = ( C1 ) 2248* ( C2 ) 2249* 2250* W := C' * V = (C1'*V1 + C2'*V2) (stored in WORK) 2251* 2252* W := C2' 2253* 2254 DO 70 J = 1, K 2255 CALL DCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 ) 2256 70 CONTINUE 2257* 2258* W := W * V2 2259* 2260 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N, 2261 $ K, ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK ) 2262 IF( M.GT.K ) THEN 2263* 2264* W := W + C1'*V1 2265* 2266 CALL DGEMM( 'Transpose', 'No transpose', N, K, M-K, 2267 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) 2268 END IF 2269* 2270* W := W * T' or W * T 2271* 2272 CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, 2273 $ ONE, T, LDT, WORK, LDWORK ) 2274* 2275* C := C - V * W' 2276* 2277 IF( M.GT.K ) THEN 2278* 2279* C1 := C1 - V1 * W' 2280* 2281 CALL DGEMM( 'No transpose', 'Transpose', M-K, N, K, 2282 $ -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC ) 2283 END IF 2284* 2285* W := W * V2' 2286* 2287 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', N, K, 2288 $ ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK ) 2289* 2290* C2 := C2 - W' 2291* 2292 DO 90 J = 1, K 2293 DO 80 I = 1, N 2294 C( M-K+J, I ) = C( M-K+J, I ) - WORK( I, J ) 2295 80 CONTINUE 2296 90 CONTINUE 2297* 2298 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 2299* 2300* Form C * H or C * H' where C = ( C1 C2 ) 2301* 2302* W := C * V = (C1*V1 + C2*V2) (stored in WORK) 2303* 2304* W := C2 2305* 2306 DO 100 J = 1, K 2307 CALL DCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 ) 2308 100 CONTINUE 2309* 2310* W := W * V2 2311* 2312 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M, 2313 $ K, ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK ) 2314 IF( N.GT.K ) THEN 2315* 2316* W := W + C1 * V1 2317* 2318 CALL DGEMM( 'No transpose', 'No transpose', M, K, N-K, 2319 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) 2320 END IF 2321* 2322* W := W * T or W * T' 2323* 2324 CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, 2325 $ ONE, T, LDT, WORK, LDWORK ) 2326* 2327* C := C - W * V' 2328* 2329 IF( N.GT.K ) THEN 2330* 2331* C1 := C1 - W * V1' 2332* 2333 CALL DGEMM( 'No transpose', 'Transpose', M, N-K, K, 2334 $ -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC ) 2335 END IF 2336* 2337* W := W * V2' 2338* 2339 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', M, K, 2340 $ ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK ) 2341* 2342* C2 := C2 - W 2343* 2344 DO 120 J = 1, K 2345 DO 110 I = 1, M 2346 C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J ) 2347 110 CONTINUE 2348 120 CONTINUE 2349 END IF 2350 END IF 2351* 2352 ELSE IF( LSAME( STOREV, 'R' ) ) THEN 2353* 2354 IF( LSAME( DIRECT, 'F' ) ) THEN 2355* 2356* Let V = ( V1 V2 ) (V1: first K columns) 2357* where V1 is unit upper triangular. 2358* 2359 IF( LSAME( SIDE, 'L' ) ) THEN 2360* 2361* Form H * C or H' * C where C = ( C1 ) 2362* ( C2 ) 2363* 2364* W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) 2365* 2366* W := C1' 2367* 2368 DO 130 J = 1, K 2369 CALL DCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 2370 130 CONTINUE 2371* 2372* W := W * V1' 2373* 2374 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', N, K, 2375 $ ONE, V, LDV, WORK, LDWORK ) 2376 IF( M.GT.K ) THEN 2377* 2378* W := W + C2'*V2' 2379* 2380 CALL DGEMM( 'Transpose', 'Transpose', N, K, M-K, ONE, 2381 $ C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, ONE, 2382 $ WORK, LDWORK ) 2383 END IF 2384* 2385* W := W * T' or W * T 2386* 2387 CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K, 2388 $ ONE, T, LDT, WORK, LDWORK ) 2389* 2390* C := C - V' * W' 2391* 2392 IF( M.GT.K ) THEN 2393* 2394* C2 := C2 - V2' * W' 2395* 2396 CALL DGEMM( 'Transpose', 'Transpose', M-K, N, K, -ONE, 2397 $ V( 1, K+1 ), LDV, WORK, LDWORK, ONE, 2398 $ C( K+1, 1 ), LDC ) 2399 END IF 2400* 2401* W := W * V1 2402* 2403 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N, 2404 $ K, ONE, V, LDV, WORK, LDWORK ) 2405* 2406* C1 := C1 - W' 2407* 2408 DO 150 J = 1, K 2409 DO 140 I = 1, N 2410 C( J, I ) = C( J, I ) - WORK( I, J ) 2411 140 CONTINUE 2412 150 CONTINUE 2413* 2414 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 2415* 2416* Form C * H or C * H' where C = ( C1 C2 ) 2417* 2418* W := C * V' = (C1*V1' + C2*V2') (stored in WORK) 2419* 2420* W := C1 2421* 2422 DO 160 J = 1, K 2423 CALL DCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 ) 2424 160 CONTINUE 2425* 2426* W := W * V1' 2427* 2428 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', M, K, 2429 $ ONE, V, LDV, WORK, LDWORK ) 2430 IF( N.GT.K ) THEN 2431* 2432* W := W + C2 * V2' 2433* 2434 CALL DGEMM( 'No transpose', 'Transpose', M, K, N-K, 2435 $ ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV, 2436 $ ONE, WORK, LDWORK ) 2437 END IF 2438* 2439* W := W * T or W * T' 2440* 2441 CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K, 2442 $ ONE, T, LDT, WORK, LDWORK ) 2443* 2444* C := C - W * V 2445* 2446 IF( N.GT.K ) THEN 2447* 2448* C2 := C2 - W * V2 2449* 2450 CALL DGEMM( 'No transpose', 'No transpose', M, N-K, K, 2451 $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, ONE, 2452 $ C( 1, K+1 ), LDC ) 2453 END IF 2454* 2455* W := W * V1 2456* 2457 CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M, 2458 $ K, ONE, V, LDV, WORK, LDWORK ) 2459* 2460* C1 := C1 - W 2461* 2462 DO 180 J = 1, K 2463 DO 170 I = 1, M 2464 C( I, J ) = C( I, J ) - WORK( I, J ) 2465 170 CONTINUE 2466 180 CONTINUE 2467* 2468 END IF 2469* 2470 ELSE 2471* 2472* Let V = ( V1 V2 ) (V2: last K columns) 2473* where V2 is unit lower triangular. 2474* 2475 IF( LSAME( SIDE, 'L' ) ) THEN 2476* 2477* Form H * C or H' * C where C = ( C1 ) 2478* ( C2 ) 2479* 2480* W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) 2481* 2482* W := C2' 2483* 2484 DO 190 J = 1, K 2485 CALL DCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 ) 2486 190 CONTINUE 2487* 2488* W := W * V2' 2489* 2490 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', N, K, 2491 $ ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK ) 2492 IF( M.GT.K ) THEN 2493* 2494* W := W + C1'*V1' 2495* 2496 CALL DGEMM( 'Transpose', 'Transpose', N, K, M-K, ONE, 2497 $ C, LDC, V, LDV, ONE, WORK, LDWORK ) 2498 END IF 2499* 2500* W := W * T' or W * T 2501* 2502 CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, 2503 $ ONE, T, LDT, WORK, LDWORK ) 2504* 2505* C := C - V' * W' 2506* 2507 IF( M.GT.K ) THEN 2508* 2509* C1 := C1 - V1' * W' 2510* 2511 CALL DGEMM( 'Transpose', 'Transpose', M-K, N, K, -ONE, 2512 $ V, LDV, WORK, LDWORK, ONE, C, LDC ) 2513 END IF 2514* 2515* W := W * V2 2516* 2517 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N, 2518 $ K, ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK ) 2519* 2520* C2 := C2 - W' 2521* 2522 DO 210 J = 1, K 2523 DO 200 I = 1, N 2524 C( M-K+J, I ) = C( M-K+J, I ) - WORK( I, J ) 2525 200 CONTINUE 2526 210 CONTINUE 2527* 2528 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 2529* 2530* Form C * H or C * H' where C = ( C1 C2 ) 2531* 2532* W := C * V' = (C1*V1' + C2*V2') (stored in WORK) 2533* 2534* W := C2 2535* 2536 DO 220 J = 1, K 2537 CALL DCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 ) 2538 220 CONTINUE 2539* 2540* W := W * V2' 2541* 2542 CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', M, K, 2543 $ ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK ) 2544 IF( N.GT.K ) THEN 2545* 2546* W := W + C1 * V1' 2547* 2548 CALL DGEMM( 'No transpose', 'Transpose', M, K, N-K, 2549 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) 2550 END IF 2551* 2552* W := W * T or W * T' 2553* 2554 CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, 2555 $ ONE, T, LDT, WORK, LDWORK ) 2556* 2557* C := C - W * V 2558* 2559 IF( N.GT.K ) THEN 2560* 2561* C1 := C1 - W * V1 2562* 2563 CALL DGEMM( 'No transpose', 'No transpose', M, N-K, K, 2564 $ -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC ) 2565 END IF 2566* 2567* W := W * V2 2568* 2569 CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M, 2570 $ K, ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK ) 2571* 2572* C1 := C1 - W 2573* 2574 DO 240 J = 1, K 2575 DO 230 I = 1, M 2576 C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J ) 2577 230 CONTINUE 2578 240 CONTINUE 2579* 2580 END IF 2581* 2582 END IF 2583 END IF 2584* 2585 RETURN 2586* 2587* End of DLARFB 2588* 2589 END 2590 SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) 2591* 2592* -- LAPACK auxiliary routine (version 1.1) -- 2593* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 2594* Courant Institute, Argonne National Lab, and Rice University 2595* February 29, 1992 2596* 2597* .. Scalar Arguments .. 2598 CHARACTER DIRECT, STOREV 2599 INTEGER K, LDT, LDV, N 2600* .. 2601* .. Array Arguments .. 2602 DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) 2603* .. 2604* 2605* Purpose 2606* ======= 2607* 2608* DLARFT forms the triangular factor T of a real block reflector H 2609* of order n, which is defined as a product of k elementary reflectors. 2610* 2611* If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; 2612* 2613* If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. 2614* 2615* If STOREV = 'C', the vector which defines the elementary reflector 2616* H(i) is stored in the i-th column of the array V, and 2617* 2618* H = I - V * T * V' 2619* 2620* If STOREV = 'R', the vector which defines the elementary reflector 2621* H(i) is stored in the i-th row of the array V, and 2622* 2623* H = I - V' * T * V 2624* 2625* Arguments 2626* ========= 2627* 2628* DIRECT (input) CHARACTER*1 2629* Specifies the order in which the elementary reflectors are 2630* multiplied to form the block reflector: 2631* = 'F': H = H(1) H(2) . . . H(k) (Forward) 2632* = 'B': H = H(k) . . . H(2) H(1) (Backward) 2633* 2634* STOREV (input) CHARACTER*1 2635* Specifies how the vectors which define the elementary 2636* reflectors are stored (see also Further Details): 2637* = 'C': columnwise 2638* = 'R': rowwise 2639* 2640* N (input) INTEGER 2641* The order of the block reflector H. N >= 0. 2642* 2643* K (input) INTEGER 2644* The order of the triangular factor T (= the number of 2645* elementary reflectors). K >= 1. 2646* 2647* V (input/output) DOUBLE PRECISION array, dimension 2648* (LDV,K) if STOREV = 'C' 2649* (LDV,N) if STOREV = 'R' 2650* The matrix V. See further details. 2651* 2652* LDV (input) INTEGER 2653* The leading dimension of the array V. 2654* If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. 2655* 2656* TAU (input) DOUBLE PRECISION array, dimension (K) 2657* TAU(i) must contain the scalar factor of the elementary 2658* reflector H(i). 2659* 2660* T (output) DOUBLE PRECISION array, dimension (LDT,K) 2661* The k by k triangular factor T of the block reflector. 2662* If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is 2663* lower triangular. The rest of the array is not used. 2664* 2665* LDT (input) INTEGER 2666* The leading dimension of the array T. LDT >= K. 2667* 2668* Further Details 2669* =============== 2670* 2671* The shape of the matrix V and the storage of the vectors which define 2672* the H(i) is best illustrated by the following example with n = 5 and 2673* k = 3. The elements equal to 1 are not stored; the corresponding 2674* array elements are modified but restored on exit. The rest of the 2675* array is not used. 2676* 2677* DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': 2678* 2679* V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) 2680* ( v1 1 ) ( 1 v2 v2 v2 ) 2681* ( v1 v2 1 ) ( 1 v3 v3 ) 2682* ( v1 v2 v3 ) 2683* ( v1 v2 v3 ) 2684* 2685* DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': 2686* 2687* V = ( v1 v2 v3 ) V = ( v1 v1 1 ) 2688* ( v1 v2 v3 ) ( v2 v2 v2 1 ) 2689* ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) 2690* ( 1 v3 ) 2691* ( 1 ) 2692* 2693* ===================================================================== 2694* 2695* .. Parameters .. 2696 DOUBLE PRECISION ONE, ZERO 2697 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 2698* .. 2699* .. Local Scalars .. 2700 INTEGER I, J 2701 DOUBLE PRECISION VII 2702* .. 2703* .. External Subroutines .. 2704 EXTERNAL DGEMV, DTRMV 2705* .. 2706* .. External Functions .. 2707 LOGICAL LSAME 2708 EXTERNAL LSAME 2709* .. 2710* .. Executable Statements .. 2711* 2712* Quick return if possible 2713* 2714 IF( N.EQ.0 ) 2715 $ RETURN 2716* 2717 IF( LSAME( DIRECT, 'F' ) ) THEN 2718 DO 20 I = 1, K 2719 IF( TAU( I ).EQ.ZERO ) THEN 2720* 2721* H(i) = I 2722* 2723 DO 10 J = 1, I 2724 T( J, I ) = ZERO 2725 10 CONTINUE 2726 ELSE 2727* 2728* general case 2729* 2730 VII = V( I, I ) 2731 V( I, I ) = ONE 2732 IF( LSAME( STOREV, 'C' ) ) THEN 2733* 2734* T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i) 2735* 2736 CALL DGEMV( 'Transpose', N-I+1, I-1, -TAU( I ), 2737 $ V( I, 1 ), LDV, V( I, I ), 1, ZERO, 2738 $ T( 1, I ), 1 ) 2739 ELSE 2740* 2741* T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)' 2742* 2743 CALL DGEMV( 'No transpose', I-1, N-I+1, -TAU( I ), 2744 $ V( 1, I ), LDV, V( I, I ), LDV, ZERO, 2745 $ T( 1, I ), 1 ) 2746 END IF 2747 V( I, I ) = VII 2748* 2749* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) 2750* 2751 CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, 2752 $ LDT, T( 1, I ), 1 ) 2753 T( I, I ) = TAU( I ) 2754 END IF 2755 20 CONTINUE 2756 ELSE 2757 DO 40 I = K, 1, -1 2758 IF( TAU( I ).EQ.ZERO ) THEN 2759* 2760* H(i) = I 2761* 2762 DO 30 J = I, K 2763 T( J, I ) = ZERO 2764 30 CONTINUE 2765 ELSE 2766* 2767* general case 2768* 2769 IF( I.LT.K ) THEN 2770 IF( LSAME( STOREV, 'C' ) ) THEN 2771 VII = V( N-K+I, I ) 2772 V( N-K+I, I ) = ONE 2773* 2774* T(i+1:k,i) := 2775* - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i) 2776* 2777 CALL DGEMV( 'Transpose', N-K+I, K-I, -TAU( I ), 2778 $ V( 1, I+1 ), LDV, V( 1, I ), 1, ZERO, 2779 $ T( I+1, I ), 1 ) 2780 V( N-K+I, I ) = VII 2781 ELSE 2782 VII = V( I, N-K+I ) 2783 V( I, N-K+I ) = ONE 2784* 2785* T(i+1:k,i) := 2786* - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)' 2787* 2788 CALL DGEMV( 'No transpose', K-I, N-K+I, -TAU( I ), 2789 $ V( I+1, 1 ), LDV, V( I, 1 ), LDV, ZERO, 2790 $ T( I+1, I ), 1 ) 2791 V( I, N-K+I ) = VII 2792 END IF 2793* 2794* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) 2795* 2796 CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, 2797 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) 2798 END IF 2799 T( I, I ) = TAU( I ) 2800 END IF 2801 40 CONTINUE 2802 END IF 2803 RETURN 2804* 2805* End of DLARFT 2806* 2807 END 2808 SUBROUTINE DORG2L( M, N, K, A, LDA, TAU, WORK, INFO ) 2809* 2810* -- LAPACK routine (version 1.1) -- 2811* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 2812* Courant Institute, Argonne National Lab, and Rice University 2813* February 29, 1992 2814* 2815* .. Scalar Arguments .. 2816 INTEGER INFO, K, LDA, M, N 2817* .. 2818* .. Array Arguments .. 2819 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) 2820* .. 2821* 2822* Purpose 2823* ======= 2824* 2825* DORG2L generates an m by n real matrix Q with orthonormal columns, 2826* which is defined as the last n columns of a product of k elementary 2827* reflectors of order m 2828* 2829* Q = H(k) . . . H(2) H(1) 2830* 2831* as returned by DGEQLF. 2832* 2833* Arguments 2834* ========= 2835* 2836* M (input) INTEGER 2837* The number of rows of the matrix Q. M >= 0. 2838* 2839* N (input) INTEGER 2840* The number of columns of the matrix Q. M >= N >= 0. 2841* 2842* K (input) INTEGER 2843* The number of elementary reflectors whose product defines the 2844* matrix Q. N >= K >= 0. 2845* 2846* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 2847* On entry, the (n-k+i)-th column must contain the vector which 2848* defines the elementary reflector H(i), for i = 1,2,...,k, as 2849* returned by DGEQLF in the last k columns of its array 2850* argument A. 2851* On exit, the m by n matrix Q. 2852* 2853* LDA (input) INTEGER 2854* The first dimension of the array A. LDA >= max(1,M). 2855* 2856* TAU (input) DOUBLE PRECISION array, dimension (K) 2857* TAU(i) must contain the scalar factor of the elementary 2858* reflector H(i), as returned by DGEQLF. 2859* 2860* WORK (workspace) DOUBLE PRECISION array, dimension (N) 2861* 2862* INFO (output) INTEGER 2863* = 0: successful exit 2864* < 0: if INFO = -i, the i-th argument has an illegal value 2865* 2866* ===================================================================== 2867* 2868* .. Parameters .. 2869 DOUBLE PRECISION ONE, ZERO 2870 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 2871* .. 2872* .. Local Scalars .. 2873 INTEGER I, II, J, L 2874* .. 2875* .. External Subroutines .. 2876 EXTERNAL DLARF, DSCAL, XERBLA 2877* .. 2878* .. Intrinsic Functions .. 2879 INTRINSIC MAX 2880* .. 2881* .. Executable Statements .. 2882* 2883* Test the input arguments 2884* 2885 INFO = 0 2886 IF( M.LT.0 ) THEN 2887 INFO = -1 2888 ELSE IF( N.LT.0 .OR. N.GT.M ) THEN 2889 INFO = -2 2890 ELSE IF( K.LT.0 .OR. K.GT.N ) THEN 2891 INFO = -3 2892 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 2893 INFO = -5 2894 END IF 2895 IF( INFO.NE.0 ) THEN 2896 CALL XERBLA( 'DORG2L', -INFO ) 2897 RETURN 2898 END IF 2899* 2900* Quick return if possible 2901* 2902 IF( N.LE.0 ) 2903 $ RETURN 2904* 2905* Initialise columns 1:n-k to columns of the unit matrix 2906* 2907 DO 20 J = 1, N - K 2908 DO 10 L = 1, M 2909 A( L, J ) = ZERO 2910 10 CONTINUE 2911 A( M-N+J, J ) = ONE 2912 20 CONTINUE 2913* 2914 DO 40 I = 1, K 2915 II = N - K + I 2916* 2917* Apply H(i) to A(1:m-k+i,1:n-k+i) from the left 2918* 2919 A( M-N+II, II ) = ONE 2920 CALL DLARF( 'Left', M-N+II, II-1, A( 1, II ), 1, TAU( I ), A, 2921 $ LDA, WORK ) 2922 CALL DSCAL( M-N+II-1, -TAU( I ), A( 1, II ), 1 ) 2923 A( M-N+II, II ) = ONE - TAU( I ) 2924* 2925* Set A(m-k+i+1:m,n-k+i) to zero 2926* 2927 DO 30 L = M - N + II + 1, M 2928 A( L, II ) = ZERO 2929 30 CONTINUE 2930 40 CONTINUE 2931 RETURN 2932* 2933* End of DORG2L 2934* 2935 END 2936 SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) 2937* 2938* -- LAPACK auxiliary routine (version 1.1) -- 2939* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 2940* Courant Institute, Argonne National Lab, and Rice University 2941* February 29, 1992 2942* 2943* .. Scalar Arguments .. 2944 CHARACTER SIDE 2945 INTEGER INCV, LDC, M, N 2946 DOUBLE PRECISION TAU 2947* .. 2948* .. Array Arguments .. 2949 DOUBLE PRECISION C( LDC, * ), V( * ), WORK( * ) 2950* .. 2951* 2952* Purpose 2953* ======= 2954* 2955* DLARF applies a real elementary reflector H to a real m by n matrix 2956* C, from either the left or the right. H is represented in the form 2957* 2958* H = I - tau * v * v' 2959* 2960* where tau is a real scalar and v is a real vector. 2961* 2962* If tau = 0, then H is taken to be the unit matrix. 2963* 2964* Arguments 2965* ========= 2966* 2967* SIDE (input) CHARACTER*1 2968* = 'L': form H * C 2969* = 'R': form C * H 2970* 2971* M (input) INTEGER 2972* The number of rows of the matrix C. 2973* 2974* N (input) INTEGER 2975* The number of columns of the matrix C. 2976* 2977* V (input) DOUBLE PRECISION array, dimension 2978* (1 + (M-1)*abs(INCV)) if SIDE = 'L' 2979* or (1 + (N-1)*abs(INCV)) if SIDE = 'R' 2980* The vector v in the representation of H. V is not used if 2981* TAU = 0. 2982* 2983* INCV (input) INTEGER 2984* The increment between elements of v. INCV <> 0. 2985* 2986* TAU (input) DOUBLE PRECISION 2987* The value tau in the representation of H. 2988* 2989* C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 2990* On entry, the m by n matrix C. 2991* On exit, C is overwritten by the matrix H * C if SIDE = 'L', 2992* or C * H if SIDE = 'R'. 2993* 2994* LDC (input) INTEGER 2995* The leading dimension of the array C. LDC >= max(1,M). 2996* 2997* WORK (workspace) DOUBLE PRECISION array, dimension 2998* (N) if SIDE = 'L' 2999* or (M) if SIDE = 'R' 3000* 3001* ===================================================================== 3002* 3003* .. Parameters .. 3004 DOUBLE PRECISION ONE, ZERO 3005 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 3006* .. 3007* .. External Subroutines .. 3008 EXTERNAL DGEMV, DGER 3009* .. 3010* .. External Functions .. 3011 LOGICAL LSAME 3012 EXTERNAL LSAME 3013* .. 3014* .. Executable Statements .. 3015* 3016 IF( LSAME( SIDE, 'L' ) ) THEN 3017* 3018* Form H * C 3019* 3020 IF( TAU.NE.ZERO ) THEN 3021* 3022* w := C' * v 3023* 3024 CALL DGEMV( 'Transpose', M, N, ONE, C, LDC, V, INCV, ZERO, 3025 $ WORK, 1 ) 3026* 3027* C := C - v * w' 3028* 3029 CALL DGER( M, N, -TAU, V, INCV, WORK, 1, C, LDC ) 3030 END IF 3031 ELSE 3032* 3033* Form C * H 3034* 3035 IF( TAU.NE.ZERO ) THEN 3036* 3037* w := C * v 3038* 3039 CALL DGEMV( 'No transpose', M, N, ONE, C, LDC, V, INCV, 3040 $ ZERO, WORK, 1 ) 3041* 3042* C := C - w * v' 3043* 3044 CALL DGER( M, N, -TAU, WORK, 1, V, INCV, C, LDC ) 3045 END IF 3046 END IF 3047 RETURN 3048* 3049* End of DLARF 3050* 3051 END 3052 SUBROUTINE DSTERF( N, D, E, INFO ) 3053* 3054* -- LAPACK routine (version 1.1) -- 3055* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 3056* Courant Institute, Argonne National Lab, and Rice University 3057* March 31, 1993 3058* 3059* .. Scalar Arguments .. 3060 INTEGER INFO, N 3061* .. 3062* .. Array Arguments .. 3063 DOUBLE PRECISION D( * ), E( * ) 3064* .. 3065* 3066* Purpose 3067* ======= 3068* 3069* DSTERF computes all eigenvalues of a symmetric tridiagonal matrix 3070* using the Pal-Walker-Kahan variant of the QL or QR algorithm. 3071* 3072* Arguments 3073* ========= 3074* 3075* N (input) INTEGER 3076* The order of the matrix. N >= 0. 3077* 3078* D (input/output) DOUBLE PRECISION array, dimension (N) 3079* On entry, the n diagonal elements of the tridiagonal matrix. 3080* On exit, if INFO = 0, the eigenvalues in ascending order. 3081* 3082* E (input/output) DOUBLE PRECISION array, dimension (N-1) 3083* On entry, the (n-1) subdiagonal elements of the tridiagonal 3084* matrix. 3085* On exit, E has been destroyed. 3086* 3087* INFO (output) INTEGER 3088* = 0: successful exit 3089* < 0: if INFO = -i, the i-th argument had an illegal value 3090* > 0: the algorithm failed to find all of the eigenvalues in 3091* a total of 30*N iterations; if INFO = i, then i 3092* elements of E have not converged to zero. 3093* 3094* ===================================================================== 3095* 3096* .. Parameters .. 3097 DOUBLE PRECISION ZERO, ONE, TWO 3098 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 3099 INTEGER MAXIT 3100 PARAMETER ( MAXIT = 30 ) 3101* .. 3102* .. Local Scalars .. 3103 INTEGER I, II, J, JTOT, K, L, L1, LEND, LENDM1, LENDP1, 3104 $ LM1, M, MM1, NM1, NMAXIT 3105 DOUBLE PRECISION ALPHA, BB, C, EPS, GAMMA, OLDC, OLDGAM, P, R, 3106 $ RT1, RT2, RTE, S, SIGMA, TST 3107* .. 3108* .. External Functions .. 3109 DOUBLE PRECISION DLAMCH, DLAPY2 3110 EXTERNAL DLAMCH, DLAPY2 3111* .. 3112* .. External Subroutines .. 3113 EXTERNAL DLAE2, XERBLA 3114* .. 3115* .. Intrinsic Functions .. 3116 INTRINSIC ABS, SIGN, SQRT 3117* .. 3118* .. Executable Statements .. 3119* 3120* Test the input parameters. 3121* 3122 INFO = 0 3123* 3124* Quick return if possible 3125* 3126 IF( N.LT.0 ) THEN 3127 INFO = -1 3128 CALL XERBLA( 'DSTERF', -INFO ) 3129 RETURN 3130 END IF 3131 IF( N.LE.1 ) 3132 $ RETURN 3133* 3134* Determine the unit roundoff for this environment. 3135* 3136 EPS = DLAMCH( 'E' ) 3137* 3138* Compute the eigenvalues of the tridiagonal matrix. 3139* 3140 DO 10 I = 1, N - 1 3141 E( I ) = E( I )**2 3142 10 CONTINUE 3143* 3144 NMAXIT = N*MAXIT 3145 SIGMA = ZERO 3146 JTOT = 0 3147* 3148* Determine where the matrix splits and choose QL or QR iteration 3149* for each block, according to whether top or bottom diagonal 3150* element is smaller. 3151* 3152 L1 = 1 3153 NM1 = N - 1 3154* 3155 20 CONTINUE 3156 IF( L1.GT.N ) 3157 $ GO TO 170 3158 IF( L1.GT.1 ) 3159 $ E( L1-1 ) = ZERO 3160 IF( L1.LE.NM1 ) THEN 3161 DO 30 M = L1, NM1 3162 TST = SQRT( ABS( E( M ) ) ) 3163 IF( TST.LE.EPS*( ABS( D( M ) )+ABS( D( M+1 ) ) ) ) 3164 $ GO TO 40 3165 30 CONTINUE 3166 END IF 3167 M = N 3168* 3169 40 CONTINUE 3170 L = L1 3171 LEND = M 3172 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 3173 L = LEND 3174 LEND = L1 3175 END IF 3176 L1 = M + 1 3177* 3178 IF( LEND.GE.L ) THEN 3179* 3180* QL Iteration 3181* 3182* Look for small subdiagonal element. 3183* 3184 50 CONTINUE 3185 IF( L.NE.LEND ) THEN 3186 LENDM1 = LEND - 1 3187 DO 60 M = L, LENDM1 3188 TST = SQRT( ABS( E( M ) ) ) 3189 IF( TST.LE.EPS*( ABS( D( M ) )+ABS( D( M+1 ) ) ) ) 3190 $ GO TO 70 3191 60 CONTINUE 3192 END IF 3193* 3194 M = LEND 3195* 3196 70 CONTINUE 3197 IF( M.LT.LEND ) 3198 $ E( M ) = ZERO 3199 P = D( L ) 3200 IF( M.EQ.L ) 3201 $ GO TO 90 3202* 3203* If remaining matrix is 2 by 2, use DLAE2 to compute its 3204* eigenvalues. 3205* 3206 IF( M.EQ.L+1 ) THEN 3207 RTE = SQRT( E( L ) ) 3208 CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 ) 3209 D( L ) = RT1 3210 D( L+1 ) = RT2 3211 E( L ) = ZERO 3212 L = L + 2 3213 IF( L.LE.LEND ) 3214 $ GO TO 50 3215 GO TO 20 3216 END IF 3217* 3218 IF( JTOT.EQ.NMAXIT ) 3219 $ GO TO 150 3220 JTOT = JTOT + 1 3221* 3222* Form shift. 3223* 3224 RTE = SQRT( E( L ) ) 3225 SIGMA = ( D( L+1 )-P ) / ( TWO*RTE ) 3226 R = DLAPY2( SIGMA, ONE ) 3227 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 3228* 3229 C = ONE 3230 S = ZERO 3231 GAMMA = D( M ) - SIGMA 3232 P = GAMMA*GAMMA 3233* 3234* Inner loop 3235* 3236 MM1 = M - 1 3237 DO 80 I = MM1, L, -1 3238 BB = E( I ) 3239 R = P + BB 3240 IF( I.NE.M-1 ) 3241 $ E( I+1 ) = S*R 3242 OLDC = C 3243 C = P / R 3244 S = BB / R 3245 OLDGAM = GAMMA 3246 ALPHA = D( I ) 3247 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 3248 D( I+1 ) = OLDGAM + ( ALPHA-GAMMA ) 3249 IF( C.NE.ZERO ) THEN 3250 P = ( GAMMA*GAMMA ) / C 3251 ELSE 3252 P = OLDC*BB 3253 END IF 3254 80 CONTINUE 3255* 3256 E( L ) = S*P 3257 D( L ) = SIGMA + GAMMA 3258 GO TO 50 3259* 3260* Eigenvalue found. 3261* 3262 90 CONTINUE 3263 D( L ) = P 3264* 3265 L = L + 1 3266 IF( L.LE.LEND ) 3267 $ GO TO 50 3268 GO TO 20 3269* 3270 ELSE 3271* 3272* QR Iteration 3273* 3274* Look for small superdiagonal element. 3275* 3276 100 CONTINUE 3277 IF( L.NE.LEND ) THEN 3278 LENDP1 = LEND + 1 3279 DO 110 M = L, LENDP1, -1 3280 TST = SQRT( ABS( E( M-1 ) ) ) 3281 IF( TST.LE.EPS*( ABS( D( M ) )+ABS( D( M-1 ) ) ) ) 3282 $ GO TO 120 3283 110 CONTINUE 3284 END IF 3285* 3286 M = LEND 3287* 3288 120 CONTINUE 3289 IF( M.GT.LEND ) 3290 $ E( M-1 ) = ZERO 3291 P = D( L ) 3292 IF( M.EQ.L ) 3293 $ GO TO 140 3294* 3295* If remaining matrix is 2 by 2, use DLAE2 to compute its 3296* eigenvalues. 3297* 3298 IF( M.EQ.L-1 ) THEN 3299 RTE = SQRT( E( L-1 ) ) 3300 CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 ) 3301 D( L ) = RT1 3302 D( L-1 ) = RT2 3303 E( L-1 ) = ZERO 3304 L = L - 2 3305 IF( L.GE.LEND ) 3306 $ GO TO 100 3307 GO TO 20 3308 END IF 3309* 3310 IF( JTOT.EQ.NMAXIT ) 3311 $ GO TO 150 3312 JTOT = JTOT + 1 3313* 3314* Form shift. 3315* 3316 RTE = SQRT( E( L-1 ) ) 3317 SIGMA = ( D( L-1 )-P ) / ( TWO*RTE ) 3318 R = DLAPY2( SIGMA, ONE ) 3319 SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) 3320* 3321 C = ONE 3322 S = ZERO 3323 GAMMA = D( M ) - SIGMA 3324 P = GAMMA*GAMMA 3325* 3326* Inner loop 3327* 3328 LM1 = L - 1 3329 DO 130 I = M, LM1 3330 BB = E( I ) 3331 R = P + BB 3332 IF( I.NE.M ) 3333 $ E( I-1 ) = S*R 3334 OLDC = C 3335 C = P / R 3336 S = BB / R 3337 OLDGAM = GAMMA 3338 ALPHA = D( I+1 ) 3339 GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM 3340 D( I ) = OLDGAM + ( ALPHA-GAMMA ) 3341 IF( C.NE.ZERO ) THEN 3342 P = ( GAMMA*GAMMA ) / C 3343 ELSE 3344 P = OLDC*BB 3345 END IF 3346 130 CONTINUE 3347* 3348 E( LM1 ) = S*P 3349 D( L ) = SIGMA + GAMMA 3350 GO TO 100 3351* 3352* Eigenvalue found. 3353* 3354 140 CONTINUE 3355 D( L ) = P 3356* 3357 L = L - 1 3358 IF( L.GE.LEND ) 3359 $ GO TO 100 3360 GO TO 20 3361* 3362 END IF 3363* 3364* Set error -- no convergence to an eigenvalue after a total 3365* of N*MAXIT iterations. 3366* 3367 150 CONTINUE 3368 DO 160 I = 1, N - 1 3369 IF( E( I ).NE.ZERO ) 3370 $ INFO = INFO + 1 3371 160 CONTINUE 3372 RETURN 3373* 3374* Sort eigenvalues in increasing order. 3375* 3376 170 CONTINUE 3377 DO 190 II = 2, N 3378 I = II - 1 3379 K = I 3380 P = D( I ) 3381 DO 180 J = II, N 3382 IF( D( J ).LT.P ) THEN 3383 K = J 3384 P = D( J ) 3385 END IF 3386 180 CONTINUE 3387 IF( K.NE.I ) THEN 3388 D( K ) = D( I ) 3389 D( I ) = P 3390 END IF 3391 190 CONTINUE 3392* 3393 RETURN 3394* 3395* End of DSTERF 3396* 3397 END 3398 SUBROUTINE DLAE2( A, B, C, RT1, RT2 ) 3399* 3400* -- LAPACK auxiliary routine (version 1.1) -- 3401* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 3402* Courant Institute, Argonne National Lab, and Rice University 3403* October 31, 1992 3404* 3405* .. Scalar Arguments .. 3406 DOUBLE PRECISION A, B, C, RT1, RT2 3407* .. 3408* 3409* Purpose 3410* ======= 3411* 3412* DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix 3413* [ A B ] 3414* [ B C ]. 3415* On return, RT1 is the eigenvalue of larger absolute value, and RT2 3416* is the eigenvalue of smaller absolute value. 3417* 3418* Arguments 3419* ========= 3420* 3421* A (input) DOUBLE PRECISION 3422* The (1,1) entry of the 2-by-2 matrix. 3423* 3424* B (input) DOUBLE PRECISION 3425* The (1,2) and (2,1) entries of the 2-by-2 matrix. 3426* 3427* C (input) DOUBLE PRECISION 3428* The (2,2) entry of the 2-by-2 matrix. 3429* 3430* RT1 (output) DOUBLE PRECISION 3431* The eigenvalue of larger absolute value. 3432* 3433* RT2 (output) DOUBLE PRECISION 3434* The eigenvalue of smaller absolute value. 3435* 3436* Further Details 3437* =============== 3438* 3439* RT1 is accurate to a few ulps barring over/underflow. 3440* 3441* RT2 may be inaccurate if there is massive cancellation in the 3442* determinant A*C-B*B; higher precision or correctly rounded or 3443* correctly truncated arithmetic would be needed to compute RT2 3444* accurately in all cases. 3445* 3446* Overflow is possible only if RT1 is within a factor of 5 of overflow. 3447* Underflow is harmless if the input data is 0 or exceeds 3448* underflow_threshold / macheps. 3449* 3450* ===================================================================== 3451* 3452* .. Parameters .. 3453 DOUBLE PRECISION ONE 3454 PARAMETER ( ONE = 1.0D0 ) 3455 DOUBLE PRECISION TWO 3456 PARAMETER ( TWO = 2.0D0 ) 3457 DOUBLE PRECISION ZERO 3458 PARAMETER ( ZERO = 0.0D0 ) 3459 DOUBLE PRECISION HALF 3460 PARAMETER ( HALF = 0.5D0 ) 3461* .. 3462* .. Local Scalars .. 3463 DOUBLE PRECISION AB, ACMN, ACMX, ADF, DF, RT, SM, TB 3464* .. 3465* .. Intrinsic Functions .. 3466 INTRINSIC ABS, SQRT 3467* .. 3468* .. Executable Statements .. 3469* 3470* Compute the eigenvalues 3471* 3472 SM = A + C 3473 DF = A - C 3474 ADF = ABS( DF ) 3475 TB = B + B 3476 AB = ABS( TB ) 3477 IF( ABS( A ).GT.ABS( C ) ) THEN 3478 ACMX = A 3479 ACMN = C 3480 ELSE 3481 ACMX = C 3482 ACMN = A 3483 END IF 3484 IF( ADF.GT.AB ) THEN 3485 RT = ADF*SQRT( ONE+( AB / ADF )**2 ) 3486 ELSE IF( ADF.LT.AB ) THEN 3487 RT = AB*SQRT( ONE+( ADF / AB )**2 ) 3488 ELSE 3489* 3490* Includes case AB=ADF=0 3491* 3492 RT = AB*SQRT( TWO ) 3493 END IF 3494 IF( SM.LT.ZERO ) THEN 3495 RT1 = HALF*( SM-RT ) 3496* 3497* Order of execution important. 3498* To get fully accurate smaller eigenvalue, 3499* next line needs to be executed in higher precision. 3500* 3501 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 3502 ELSE IF( SM.GT.ZERO ) THEN 3503 RT1 = HALF*( SM+RT ) 3504* 3505* Order of execution important. 3506* To get fully accurate smaller eigenvalue, 3507* next line needs to be executed in higher precision. 3508* 3509 RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B 3510 ELSE 3511* 3512* Includes case RT1 = RT2 = 0 3513* 3514 RT1 = HALF*RT 3515 RT2 = -HALF*RT 3516 END IF 3517 RETURN 3518* 3519* End of DLAE2 3520* 3521 END 3522 SUBROUTINE DSYTRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO ) 3523* 3524* -- LAPACK routine (version 1.1) -- 3525* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 3526* Courant Institute, Argonne National Lab, and Rice University 3527* March 31, 1993 3528* 3529* .. Scalar Arguments .. 3530 CHARACTER UPLO 3531 INTEGER INFO, LDA, LWORK, N 3532* .. 3533* .. Array Arguments .. 3534 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAU( * ), 3535 $ WORK( * ) 3536* .. 3537* 3538* Purpose 3539* ======= 3540* 3541* DSYTRD reduces a real symmetric matrix A to real symmetric 3542* tridiagonal form T by an orthogonal similarity transformation: 3543* Q**T * A * Q = T. 3544* 3545* Arguments 3546* ========= 3547* 3548* UPLO (input) CHARACTER*1 3549* = 'U': Upper triangle of A is stored; 3550* = 'L': Lower triangle of A is stored. 3551* 3552* N (input) INTEGER 3553* The order of the matrix A. N >= 0. 3554* 3555* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 3556* On entry, the symmetric matrix A. If UPLO = 'U', the leading 3557* N-by-N upper triangular part of A contains the upper 3558* triangular part of the matrix A, and the strictly lower 3559* triangular part of A is not referenced. If UPLO = 'L', the 3560* leading N-by-N lower triangular part of A contains the lower 3561* triangular part of the matrix A, and the strictly upper 3562* triangular part of A is not referenced. 3563* On exit, if UPLO = 'U', the diagonal and first superdiagonal 3564* of A are overwritten by the corresponding elements of the 3565* tridiagonal matrix T, and the elements above the first 3566* superdiagonal, with the array TAU, represent the orthogonal 3567* matrix Q as a product of elementary reflectors; if UPLO 3568* = 'L', the diagonal and first subdiagonal of A are over- 3569* written by the corresponding elements of the tridiagonal 3570* matrix T, and the elements below the first subdiagonal, with 3571* the array TAU, represent the orthogonal matrix Q as a product 3572* of elementary reflectors. See Further Details. 3573* 3574* LDA (input) INTEGER 3575* The leading dimension of the array A. LDA >= max(1,N). 3576* 3577* D (output) DOUBLE PRECISION array, dimension (N) 3578* The diagonal elements of the tridiagonal matrix T: 3579* D(i) = A(i,i). 3580* 3581* E (output) DOUBLE PRECISION array, dimension (N-1) 3582* The off-diagonal elements of the tridiagonal matrix T: 3583* E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 3584* 3585* TAU (output) DOUBLE PRECISION array, dimension (N-1) 3586* The scalar factors of the elementary reflectors (see Further 3587* Details). 3588* 3589* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 3590* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 3591* 3592* LWORK (input) INTEGER 3593* The dimension of the array WORK. LWORK >= 1. 3594* For optimum performance LWORK >= N*NB, where NB is the 3595* optimal blocksize. 3596* 3597* INFO (output) INTEGER 3598* = 0: successful exit 3599* < 0: if INFO = -i, the i-th argument had an illegal value 3600* 3601* Further Details 3602* =============== 3603* 3604* If UPLO = 'U', the matrix Q is represented as a product of elementary 3605* reflectors 3606* 3607* Q = H(n-1) . . . H(2) H(1). 3608* 3609* Each H(i) has the form 3610* 3611* H(i) = I - tau * v * v' 3612* 3613* where tau is a real scalar, and v is a real vector with 3614* v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in 3615* A(1:i-1,i+1), and tau in TAU(i). 3616* 3617* If UPLO = 'L', the matrix Q is represented as a product of elementary 3618* reflectors 3619* 3620* Q = H(1) H(2) . . . H(n-1). 3621* 3622* Each H(i) has the form 3623* 3624* H(i) = I - tau * v * v' 3625* 3626* where tau is a real scalar, and v is a real vector with 3627* v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), 3628* and tau in TAU(i). 3629* 3630* The contents of A on exit are illustrated by the following examples 3631* with n = 5: 3632* 3633* if UPLO = 'U': if UPLO = 'L': 3634* 3635* ( d e v2 v3 v4 ) ( d ) 3636* ( d e v3 v4 ) ( e d ) 3637* ( d e v4 ) ( v1 e d ) 3638* ( d e ) ( v1 v2 e d ) 3639* ( d ) ( v1 v2 v3 e d ) 3640* 3641* where d and e denote diagonal and off-diagonal elements of T, and vi 3642* denotes an element of the vector defining H(i). 3643* 3644* ===================================================================== 3645* 3646* .. Parameters .. 3647 DOUBLE PRECISION ONE 3648 PARAMETER ( ONE = 1.0D0 ) 3649* .. 3650* .. Local Scalars .. 3651 LOGICAL UPPER 3652 INTEGER I, IINFO, IWS, J, KK, LDWORK, NB, NBMIN, NX 3653* .. 3654* .. External Subroutines .. 3655 EXTERNAL DLATRD, DSYR2K, DSYTD2, XERBLA 3656* .. 3657* .. Intrinsic Functions .. 3658 INTRINSIC MAX 3659* .. 3660* .. External Functions .. 3661 LOGICAL LSAME 3662 INTEGER ILAENV 3663 EXTERNAL LSAME, ILAENV 3664* .. 3665* .. Executable Statements .. 3666* 3667* Test the input parameters 3668* 3669 INFO = 0 3670 UPPER = LSAME( UPLO, 'U' ) 3671 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 3672 INFO = -1 3673 ELSE IF( N.LT.0 ) THEN 3674 INFO = -2 3675 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 3676 INFO = -4 3677 ELSE IF( LWORK.LT.1 ) THEN 3678 INFO = -9 3679 END IF 3680 IF( INFO.NE.0 ) THEN 3681 CALL XERBLA( 'DSYTRD', -INFO ) 3682 RETURN 3683 END IF 3684* 3685* Quick return if possible 3686* 3687 IF( N.EQ.0 ) THEN 3688 WORK( 1 ) = 1 3689 RETURN 3690 END IF 3691* 3692* Determine the block size. 3693* 3694 NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) 3695 NX = N 3696 IWS = 1 3697 IF( NB.GT.1 .AND. NB.LT.N ) THEN 3698* 3699* Determine when to cross over from blocked to unblocked code 3700* (last block is always handled by unblocked code). 3701* 3702 NX = MAX( NB, ILAENV( 3, 'DSYTRD', UPLO, N, -1, -1, -1 ) ) 3703 IF( NX.LT.N ) THEN 3704* 3705* Determine if workspace is large enough for blocked code. 3706* 3707 LDWORK = N 3708 IWS = LDWORK*NB 3709 IF( LWORK.LT.IWS ) THEN 3710* 3711* Not enough workspace to use optimal NB: determine the 3712* minimum value of NB, and reduce NB or force use of 3713* unblocked code by setting NX = N. 3714* 3715 NB = LWORK / LDWORK 3716 NBMIN = ILAENV( 2, 'DSYTRD', UPLO, N, -1, -1, -1 ) 3717 IF( NB.LT.NBMIN ) 3718 $ NX = N 3719 END IF 3720 ELSE 3721 NX = N 3722 END IF 3723 ELSE 3724 NB = 1 3725 END IF 3726* 3727 IF( UPPER ) THEN 3728* 3729* Reduce the upper triangle of A. 3730* Columns 1:kk are handled by the unblocked method. 3731* 3732 KK = N - ( ( N-NX+NB-1 ) / NB )*NB 3733 DO 20 I = N - NB + 1, KK + 1, -NB 3734* 3735* Reduce columns i:i+nb-1 to tridiagonal form and form the 3736* matrix W which is needed to update the unreduced part of 3737* the matrix 3738* 3739 CALL DLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK, 3740 $ LDWORK ) 3741* 3742* Update the unreduced submatrix A(1:i-1,1:i-1), using an 3743* update of the form: A := A - V*W' - W*V' 3744* 3745 CALL DSYR2K( UPLO, 'No transpose', I-1, NB, -ONE, A( 1, I ), 3746 $ LDA, WORK, LDWORK, ONE, A, LDA ) 3747* 3748* Copy superdiagonal elements back into A, and diagonal 3749* elements into D 3750* 3751 DO 10 J = I, I + NB - 1 3752 A( J-1, J ) = E( J-1 ) 3753 D( J ) = A( J, J ) 3754 10 CONTINUE 3755 20 CONTINUE 3756* 3757* Use unblocked code to reduce the last or only block 3758* 3759 CALL DSYTD2( UPLO, KK, A, LDA, D, E, TAU, IINFO ) 3760 ELSE 3761* 3762* Reduce the lower triangle of A 3763* 3764 DO 40 I = 1, N - NX, NB 3765* 3766* Reduce columns i:i+nb-1 to tridiagonal form and form the 3767* matrix W which is needed to update the unreduced part of 3768* the matrix 3769* 3770 CALL DLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ), 3771 $ TAU( I ), WORK, LDWORK ) 3772* 3773* Update the unreduced submatrix A(i+ib:n,i+ib:n), using 3774* an update of the form: A := A - V*W' - W*V' 3775* 3776 CALL DSYR2K( UPLO, 'No transpose', N-I-NB+1, NB, -ONE, 3777 $ A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE, 3778 $ A( I+NB, I+NB ), LDA ) 3779* 3780* Copy subdiagonal elements back into A, and diagonal 3781* elements into D 3782* 3783 DO 30 J = I, I + NB - 1 3784 A( J+1, J ) = E( J ) 3785 D( J ) = A( J, J ) 3786 30 CONTINUE 3787 40 CONTINUE 3788* 3789* Use unblocked code to reduce the last or only block 3790* 3791 CALL DSYTD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ), 3792 $ TAU( I ), IINFO ) 3793 END IF 3794* 3795 WORK( 1 ) = IWS 3796 RETURN 3797* 3798* End of DSYTRD 3799* 3800 END 3801 SUBROUTINE DSYTD2( UPLO, N, A, LDA, D, E, TAU, INFO ) 3802* 3803* -- LAPACK routine (version 1.1) -- 3804* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 3805* Courant Institute, Argonne National Lab, and Rice University 3806* October 31, 1992 3807* 3808* .. Scalar Arguments .. 3809 CHARACTER UPLO 3810 INTEGER INFO, LDA, N 3811* .. 3812* .. Array Arguments .. 3813 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAU( * ) 3814* .. 3815* 3816* Purpose 3817* ======= 3818* 3819* DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal 3820* form T by an orthogonal similarity transformation: Q' * A * Q = T. 3821* 3822* Arguments 3823* ========= 3824* 3825* UPLO (input) CHARACTER*1 3826* Specifies whether the upper or lower triangular part of the 3827* symmetric matrix A is stored: 3828* = 'U': Upper triangular 3829* = 'L': Lower triangular 3830* 3831* N (input) INTEGER 3832* The order of the matrix A. N >= 0. 3833* 3834* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 3835* On entry, the symmetric matrix A. If UPLO = 'U', the leading 3836* n-by-n upper triangular part of A contains the upper 3837* triangular part of the matrix A, and the strictly lower 3838* triangular part of A is not referenced. If UPLO = 'L', the 3839* leading n-by-n lower triangular part of A contains the lower 3840* triangular part of the matrix A, and the strictly upper 3841* triangular part of A is not referenced. 3842* On exit, if UPLO = 'U', the diagonal and first superdiagonal 3843* of A are overwritten by the corresponding elements of the 3844* tridiagonal matrix T, and the elements above the first 3845* superdiagonal, with the array TAU, represent the orthogonal 3846* matrix Q as a product of elementary reflectors; if UPLO 3847* = 'L', the diagonal and first subdiagonal of A are over- 3848* written by the corresponding elements of the tridiagonal 3849* matrix T, and the elements below the first subdiagonal, with 3850* the array TAU, represent the orthogonal matrix Q as a product 3851* of elementary reflectors. See Further Details. 3852* 3853* LDA (input) INTEGER 3854* The leading dimension of the array A. LDA >= max(1,N). 3855* 3856* D (output) DOUBLE PRECISION array, dimension (N) 3857* The diagonal elements of the tridiagonal matrix T: 3858* D(i) = A(i,i). 3859* 3860* E (output) DOUBLE PRECISION array, dimension (N-1) 3861* The off-diagonal elements of the tridiagonal matrix T: 3862* E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 3863* 3864* TAU (output) DOUBLE PRECISION array, dimension (N-1) 3865* The scalar factors of the elementary reflectors (see Further 3866* Details). 3867* 3868* INFO (output) INTEGER 3869* = 0: successful exit 3870* < 0: if INFO = -i, the i-th argument had an illegal value. 3871* 3872* Further Details 3873* =============== 3874* 3875* If UPLO = 'U', the matrix Q is represented as a product of elementary 3876* reflectors 3877* 3878* Q = H(n-1) . . . H(2) H(1). 3879* 3880* Each H(i) has the form 3881* 3882* H(i) = I - tau * v * v' 3883* 3884* where tau is a real scalar, and v is a real vector with 3885* v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in 3886* A(1:i-1,i+1), and tau in TAU(i). 3887* 3888* If UPLO = 'L', the matrix Q is represented as a product of elementary 3889* reflectors 3890* 3891* Q = H(1) H(2) . . . H(n-1). 3892* 3893* Each H(i) has the form 3894* 3895* H(i) = I - tau * v * v' 3896* 3897* where tau is a real scalar, and v is a real vector with 3898* v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), 3899* and tau in TAU(i). 3900* 3901* The contents of A on exit are illustrated by the following examples 3902* with n = 5: 3903* 3904* if UPLO = 'U': if UPLO = 'L': 3905* 3906* ( d e v2 v3 v4 ) ( d ) 3907* ( d e v3 v4 ) ( e d ) 3908* ( d e v4 ) ( v1 e d ) 3909* ( d e ) ( v1 v2 e d ) 3910* ( d ) ( v1 v2 v3 e d ) 3911* 3912* where d and e denote diagonal and off-diagonal elements of T, and vi 3913* denotes an element of the vector defining H(i). 3914* 3915* ===================================================================== 3916* 3917* .. Parameters .. 3918 DOUBLE PRECISION ONE, ZERO, HALF 3919 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, 3920 $ HALF = 1.0D0 / 2.0D0 ) 3921* .. 3922* .. Local Scalars .. 3923 LOGICAL UPPER 3924 INTEGER I 3925 DOUBLE PRECISION ALPHA, TAUI 3926* .. 3927* .. External Subroutines .. 3928 EXTERNAL DAXPY, DLARFG, DSYMV, DSYR2, XERBLA 3929* .. 3930* .. External Functions .. 3931 LOGICAL LSAME 3932 DOUBLE PRECISION DDOT 3933 EXTERNAL LSAME, DDOT 3934* .. 3935* .. Intrinsic Functions .. 3936 INTRINSIC MAX, MIN 3937* .. 3938* .. Executable Statements .. 3939* 3940* Test the input parameters 3941* 3942 INFO = 0 3943 UPPER = LSAME( UPLO, 'U' ) 3944 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 3945 INFO = -1 3946 ELSE IF( N.LT.0 ) THEN 3947 INFO = -2 3948 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 3949 INFO = -4 3950 END IF 3951 IF( INFO.NE.0 ) THEN 3952 CALL XERBLA( 'DSYTD2', -INFO ) 3953 RETURN 3954 END IF 3955* 3956* Quick return if possible 3957* 3958 IF( N.LE.0 ) 3959 $ RETURN 3960* 3961 IF( UPPER ) THEN 3962* 3963* Reduce the upper triangle of A 3964* 3965 DO 10 I = N - 1, 1, -1 3966* 3967* Generate elementary reflector H(i) = I - tau * v * v' 3968* to annihilate A(1:i-1,i+1) 3969* 3970 CALL DLARFG( I, A( I, I+1 ), A( 1, I+1 ), 1, TAUI ) 3971 E( I ) = A( I, I+1 ) 3972* 3973 IF( TAUI.NE.ZERO ) THEN 3974* 3975* Apply H(i) from both sides to A(1:i,1:i) 3976* 3977 A( I, I+1 ) = ONE 3978* 3979* Compute x := tau * A * v storing x in TAU(1:i) 3980* 3981 CALL DSYMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO, 3982 $ TAU, 1 ) 3983* 3984* Compute w := x - 1/2 * tau * (x'*v) * v 3985* 3986 ALPHA = -HALF*TAUI*DDOT( I, TAU, 1, A( 1, I+1 ), 1 ) 3987 CALL DAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 ) 3988* 3989* Apply the transformation as a rank-2 update: 3990* A := A - v * w' - w * v' 3991* 3992 CALL DSYR2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, 3993 $ LDA ) 3994* 3995 A( I, I+1 ) = E( I ) 3996 END IF 3997 D( I+1 ) = A( I+1, I+1 ) 3998 TAU( I ) = TAUI 3999 10 CONTINUE 4000 D( 1 ) = A( 1, 1 ) 4001 ELSE 4002* 4003* Reduce the lower triangle of A 4004* 4005 DO 20 I = 1, N - 1 4006* 4007* Generate elementary reflector H(i) = I - tau * v * v' 4008* to annihilate A(i+2:n,i) 4009* 4010 CALL DLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1, 4011 $ TAUI ) 4012 E( I ) = A( I+1, I ) 4013* 4014 IF( TAUI.NE.ZERO ) THEN 4015* 4016* Apply H(i) from both sides to A(i+1:n,i+1:n) 4017* 4018 A( I+1, I ) = ONE 4019* 4020* Compute x := tau * A * v storing y in TAU(i:n-1) 4021* 4022 CALL DSYMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA, 4023 $ A( I+1, I ), 1, ZERO, TAU( I ), 1 ) 4024* 4025* Compute w := x - 1/2 * tau * (x'*v) * v 4026* 4027 ALPHA = -HALF*TAUI*DDOT( N-I, TAU( I ), 1, A( I+1, I ), 4028 $ 1 ) 4029 CALL DAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 ) 4030* 4031* Apply the transformation as a rank-2 update: 4032* A := A - v * w' - w * v' 4033* 4034 CALL DSYR2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1, 4035 $ A( I+1, I+1 ), LDA ) 4036* 4037 A( I+1, I ) = E( I ) 4038 END IF 4039 D( I ) = A( I, I ) 4040 TAU( I ) = TAUI 4041 20 CONTINUE 4042 D( N ) = A( N, N ) 4043 END IF 4044* 4045 RETURN 4046* 4047* End of DSYTD2 4048* 4049 END 4050 SUBROUTINE XERBLA( SRNAME, INFO ) 4051* 4052* -- LAPACK auxiliary routine (version 1.1) -- 4053* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 4054* Courant Institute, Argonne National Lab, and Rice University 4055* February 29, 1992 4056* 4057* .. Scalar Arguments .. 4058 CHARACTER*6 SRNAME 4059 INTEGER INFO 4060* .. 4061* 4062* Purpose 4063* ======= 4064* 4065* XERBLA is an error handler for the LAPACK routines. 4066* It is called by an LAPACK routine if an input parameter has an 4067* invalid value. A message is printed and execution stops. 4068* 4069* Installers may consider modifying the STOP statement in order to 4070* call system-specific exception-handling facilities. 4071* 4072* Arguments 4073* ========= 4074* 4075* SRNAME (input) CHARACTER*6 4076* The name of the routine which called XERBLA. 4077* 4078* INFO (input) INTEGER 4079* The position of the invalid parameter in the parameter list 4080* of the calling routine. 4081* 4082* .. Executable Statements .. 4083* 4084 WRITE( *, FMT = 9999 )SRNAME, INFO 4085* 4086 STOP 4087* 4088 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', 4089 $ 'an illegal value' ) 4090* 4091* End of XERBLA 4092* 4093 END 4094 SUBROUTINE DLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW ) 4095* 4096* -- LAPACK auxiliary routine (version 1.1) -- 4097* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 4098* Courant Institute, Argonne National Lab, and Rice University 4099* October 31, 1992 4100* 4101* .. Scalar Arguments .. 4102 CHARACTER UPLO 4103 INTEGER LDA, LDW, N, NB 4104* .. 4105* .. Array Arguments .. 4106 DOUBLE PRECISION A( LDA, * ), E( * ), TAU( * ), W( LDW, * ) 4107* .. 4108* 4109* Purpose 4110* ======= 4111* 4112* DLATRD reduces NB rows and columns of a real symmetric matrix A to 4113* symmetric tridiagonal form by an orthogonal similarity 4114* transformation Q' * A * Q, and returns the matrices V and W which are 4115* needed to apply the transformation to the unreduced part of A. 4116* 4117* If UPLO = 'U', DLATRD reduces the last NB rows and columns of a 4118* matrix, of which the upper triangle is supplied; 4119* if UPLO = 'L', DLATRD reduces the first NB rows and columns of a 4120* matrix, of which the lower triangle is supplied. 4121* 4122* This is an auxiliary routine called by DSYTRD. 4123* 4124* Arguments 4125* ========= 4126* 4127* UPLO (input) CHARACTER 4128* Specifies whether the upper or lower triangular part of the 4129* symmetric matrix A is stored: 4130* = 'U': Upper triangular 4131* = 'L': Lower triangular 4132* 4133* N (input) INTEGER 4134* The order of the matrix A. 4135* 4136* NB (input) INTEGER 4137* The number of rows and columns to be reduced. 4138* 4139* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 4140* On entry, the symmetric matrix A. If UPLO = 'U', the leading 4141* n-by-n upper triangular part of A contains the upper 4142* triangular part of the matrix A, and the strictly lower 4143* triangular part of A is not referenced. If UPLO = 'L', the 4144* leading n-by-n lower triangular part of A contains the lower 4145* triangular part of the matrix A, and the strictly upper 4146* triangular part of A is not referenced. 4147* On exit: 4148* if UPLO = 'U', the last NB columns have been reduced to 4149* tridiagonal form, with the diagonal elements overwriting 4150* the diagonal elements of A; the elements above the diagonal 4151* with the array TAU, represent the orthogonal matrix Q as a 4152* product of elementary reflectors; 4153* if UPLO = 'L', the first NB columns have been reduced to 4154* tridiagonal form, with the diagonal elements overwriting 4155* the diagonal elements of A; the elements below the diagonal 4156* with the array TAU, represent the orthogonal matrix Q as a 4157* product of elementary reflectors. 4158* See Further Details. 4159* 4160* LDA (input) INTEGER 4161* The leading dimension of the array A. LDA >= (1,N). 4162* 4163* E (output) DOUBLE PRECISION array, dimension (N-1) 4164* If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal 4165* elements of the last NB columns of the reduced matrix; 4166* if UPLO = 'L', E(1:nb) contains the subdiagonal elements of 4167* the first NB columns of the reduced matrix. 4168* 4169* TAU (output) DOUBLE PRECISION array, dimension (N-1) 4170* The scalar factors of the elementary reflectors, stored in 4171* TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'. 4172* See Further Details. 4173* 4174* W (output) DOUBLE PRECISION array, dimension (LDW,NB) 4175* The n-by-nb matrix W required to update the unreduced part 4176* of A. 4177* 4178* LDW (input) INTEGER 4179* The leading dimension of the array W. LDW >= max(1,N). 4180* 4181* Further Details 4182* =============== 4183* 4184* If UPLO = 'U', the matrix Q is represented as a product of elementary 4185* reflectors 4186* 4187* Q = H(n) H(n-1) . . . H(n-nb+1). 4188* 4189* Each H(i) has the form 4190* 4191* H(i) = I - tau * v * v' 4192* 4193* where tau is a real scalar, and v is a real vector with 4194* v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i), 4195* and tau in TAU(i-1). 4196* 4197* If UPLO = 'L', the matrix Q is represented as a product of elementary 4198* reflectors 4199* 4200* Q = H(1) H(2) . . . H(nb). 4201* 4202* Each H(i) has the form 4203* 4204* H(i) = I - tau * v * v' 4205* 4206* where tau is a real scalar, and v is a real vector with 4207* v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i), 4208* and tau in TAU(i). 4209* 4210* The elements of the vectors v together form the n-by-nb matrix V 4211* which is needed, with W, to apply the transformation to the unreduced 4212* part of the matrix, using a symmetric rank-2k update of the form: 4213* A := A - V*W' - W*V'. 4214* 4215* The contents of A on exit are illustrated by the following examples 4216* with n = 5 and nb = 2: 4217* 4218* if UPLO = 'U': if UPLO = 'L': 4219* 4220* ( a a a v4 v5 ) ( d ) 4221* ( a a v4 v5 ) ( 1 d ) 4222* ( a 1 v5 ) ( v1 1 a ) 4223* ( d 1 ) ( v1 v2 a a ) 4224* ( d ) ( v1 v2 a a a ) 4225* 4226* where d denotes a diagonal element of the reduced matrix, a denotes 4227* an element of the original matrix that is unchanged, and vi denotes 4228* an element of the vector defining H(i). 4229* 4230* ===================================================================== 4231* 4232* .. Parameters .. 4233 DOUBLE PRECISION ZERO, ONE, HALF 4234 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 ) 4235* .. 4236* .. Local Scalars .. 4237 INTEGER I, IW 4238 DOUBLE PRECISION ALPHA 4239* .. 4240* .. External Subroutines .. 4241 EXTERNAL DAXPY, DGEMV, DLARFG, DSCAL, DSYMV 4242* .. 4243* .. External Functions .. 4244 LOGICAL LSAME 4245 DOUBLE PRECISION DDOT 4246 EXTERNAL LSAME, DDOT 4247* .. 4248* .. Intrinsic Functions .. 4249 INTRINSIC MIN 4250* .. 4251* .. Executable Statements .. 4252* 4253* Quick return if possible 4254* 4255 IF( N.LE.0 ) 4256 $ RETURN 4257* 4258 IF( LSAME( UPLO, 'U' ) ) THEN 4259* 4260* Reduce last NB columns of upper triangle 4261* 4262 DO 10 I = N, N - NB + 1, -1 4263 IW = I - N + NB 4264 IF( I.LT.N ) THEN 4265* 4266* Update A(1:i,i) 4267* 4268 CALL DGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ), 4269 $ LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 ) 4270 CALL DGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ), 4271 $ LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 ) 4272 END IF 4273 IF( I.GT.1 ) THEN 4274* 4275* Generate elementary reflector H(i) to annihilate 4276* A(1:i-2,i) 4277* 4278 CALL DLARFG( I-1, A( I-1, I ), A( 1, I ), 1, TAU( I-1 ) ) 4279 E( I-1 ) = A( I-1, I ) 4280 A( I-1, I ) = ONE 4281* 4282* Compute W(1:i-1,i) 4283* 4284 CALL DSYMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1, 4285 $ ZERO, W( 1, IW ), 1 ) 4286 IF( I.LT.N ) THEN 4287 CALL DGEMV( 'Transpose', I-1, N-I, ONE, W( 1, IW+1 ), 4288 $ LDW, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 ) 4289 CALL DGEMV( 'No transpose', I-1, N-I, -ONE, 4290 $ A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE, 4291 $ W( 1, IW ), 1 ) 4292 CALL DGEMV( 'Transpose', I-1, N-I, ONE, A( 1, I+1 ), 4293 $ LDA, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 ) 4294 CALL DGEMV( 'No transpose', I-1, N-I, -ONE, 4295 $ W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE, 4296 $ W( 1, IW ), 1 ) 4297 END IF 4298 CALL DSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 ) 4299 ALPHA = -HALF*TAU( I-1 )*DDOT( I-1, W( 1, IW ), 1, 4300 $ A( 1, I ), 1 ) 4301 CALL DAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 ) 4302 END IF 4303* 4304 10 CONTINUE 4305 ELSE 4306* 4307* Reduce first NB columns of lower triangle 4308* 4309 DO 20 I = 1, NB 4310* 4311* Update A(i:n,i) 4312* 4313 CALL DGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ), 4314 $ LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 ) 4315 CALL DGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ), 4316 $ LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 ) 4317 IF( I.LT.N ) THEN 4318* 4319* Generate elementary reflector H(i) to annihilate 4320* A(i+2:n,i) 4321* 4322 CALL DLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1, 4323 $ TAU( I ) ) 4324 E( I ) = A( I+1, I ) 4325 A( I+1, I ) = ONE 4326* 4327* Compute W(i+1:n,i) 4328* 4329 CALL DSYMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA, 4330 $ A( I+1, I ), 1, ZERO, W( I+1, I ), 1 ) 4331 CALL DGEMV( 'Transpose', N-I, I-1, ONE, W( I+1, 1 ), LDW, 4332 $ A( I+1, I ), 1, ZERO, W( 1, I ), 1 ) 4333 CALL DGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ), 4334 $ LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 ) 4335 CALL DGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA, 4336 $ A( I+1, I ), 1, ZERO, W( 1, I ), 1 ) 4337 CALL DGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ), 4338 $ LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 ) 4339 CALL DSCAL( N-I, TAU( I ), W( I+1, I ), 1 ) 4340 ALPHA = -HALF*TAU( I )*DDOT( N-I, W( I+1, I ), 1, 4341 $ A( I+1, I ), 1 ) 4342 CALL DAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 ) 4343 END IF 4344* 4345 20 CONTINUE 4346 END IF 4347* 4348 RETURN 4349* 4350* End of DLATRD 4351* 4352 END 4353 SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) 4354* 4355* -- LAPACK auxiliary routine (version 1.1) -- 4356* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 4357* Courant Institute, Argonne National Lab, and Rice University 4358* February 29, 1992 4359* 4360* .. Scalar Arguments .. 4361 INTEGER INCX, N 4362 DOUBLE PRECISION ALPHA, TAU 4363* .. 4364* .. Array Arguments .. 4365 DOUBLE PRECISION X( * ) 4366* .. 4367* 4368* Purpose 4369* ======= 4370* 4371* DLARFG generates a real elementary reflector H of order n, such 4372* that 4373* 4374* H * ( alpha ) = ( beta ), H' * H = I. 4375* ( x ) ( 0 ) 4376* 4377* where alpha and beta are scalars, and x is an (n-1)-element real 4378* vector. H is represented in the form 4379* 4380* H = I - tau * ( 1 ) * ( 1 v' ) , 4381* ( v ) 4382* 4383* where tau is a real scalar and v is a real (n-1)-element 4384* vector. 4385* 4386* If the elements of x are all zero, then tau = 0 and H is taken to be 4387* the unit matrix. 4388* 4389* Otherwise 1 <= tau <= 2. 4390* 4391* Arguments 4392* ========= 4393* 4394* N (input) INTEGER 4395* The order of the elementary reflector. 4396* 4397* ALPHA (input/output) DOUBLE PRECISION 4398* On entry, the value alpha. 4399* On exit, it is overwritten with the value beta. 4400* 4401* X (input/output) DOUBLE PRECISION array, dimension 4402* (1+(N-2)*abs(INCX)) 4403* On entry, the vector x. 4404* On exit, it is overwritten with the vector v. 4405* 4406* INCX (input) INTEGER 4407* The increment between elements of X. INCX <> 0. 4408* 4409* TAU (output) DOUBLE PRECISION 4410* The value tau. 4411* 4412* ===================================================================== 4413* 4414* .. Parameters .. 4415 DOUBLE PRECISION ONE, ZERO 4416 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 4417* .. 4418* .. Local Scalars .. 4419 INTEGER J, KNT 4420 DOUBLE PRECISION BETA, RSAFMN, SAFMIN, XNORM 4421* .. 4422* .. External Functions .. 4423 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 4424 EXTERNAL DLAMCH, DLAPY2, DNRM2 4425* .. 4426* .. Intrinsic Functions .. 4427 INTRINSIC ABS, SIGN 4428* .. 4429* .. External Subroutines .. 4430 EXTERNAL DSCAL 4431* .. 4432* .. Executable Statements .. 4433* 4434 IF( N.LE.1 ) THEN 4435 TAU = ZERO 4436 RETURN 4437 END IF 4438* 4439 XNORM = DNRM2( N-1, X, INCX ) 4440* 4441 IF( XNORM.EQ.ZERO ) THEN 4442* 4443* H = I 4444* 4445 TAU = ZERO 4446 ELSE 4447* 4448* general case 4449* 4450 BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) 4451 SAFMIN = DLAMCH( 'S' ) 4452 IF( ABS( BETA ).LT.SAFMIN ) THEN 4453* 4454* XNORM, BETA may be inaccurate; scale X and recompute them 4455* 4456 RSAFMN = ONE / SAFMIN 4457 KNT = 0 4458 10 CONTINUE 4459 KNT = KNT + 1 4460 CALL DSCAL( N-1, RSAFMN, X, INCX ) 4461 BETA = BETA*RSAFMN 4462 ALPHA = ALPHA*RSAFMN 4463 IF( ABS( BETA ).LT.SAFMIN ) 4464 $ GO TO 10 4465* 4466* New BETA is at most 1, at least SAFMIN 4467* 4468 XNORM = DNRM2( N-1, X, INCX ) 4469 BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) 4470 TAU = ( BETA-ALPHA ) / BETA 4471 CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) 4472* 4473* If ALPHA is subnormal, it may lose relative accuracy 4474* 4475 ALPHA = BETA 4476 DO 20 J = 1, KNT 4477 ALPHA = ALPHA*SAFMIN 4478 20 CONTINUE 4479 ELSE 4480 TAU = ( BETA-ALPHA ) / BETA 4481 CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) 4482 ALPHA = BETA 4483 END IF 4484 END IF 4485* 4486 RETURN 4487* 4488* End of DLARFG 4489* 4490 END 4491 DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) 4492* 4493* -- LAPACK auxiliary routine (version 1.1) -- 4494* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 4495* Courant Institute, Argonne National Lab, and Rice University 4496* October 31, 1992 4497* 4498* .. Scalar Arguments .. 4499 CHARACTER CMACH 4500* .. 4501* 4502* Purpose 4503* ======= 4504* 4505* DLAMCH determines double precision machine parameters. 4506* 4507* Arguments 4508* ========= 4509* 4510* CMACH (input) CHARACTER*1 4511* Specifies the value to be returned by DLAMCH: 4512* = 'E' or 'e', DLAMCH := eps 4513* = 'S' or 's , DLAMCH := sfmin 4514* = 'B' or 'b', DLAMCH := base 4515* = 'P' or 'p', DLAMCH := eps*base 4516* = 'N' or 'n', DLAMCH := t 4517* = 'R' or 'r', DLAMCH := rnd 4518* = 'M' or 'm', DLAMCH := emin 4519* = 'U' or 'u', DLAMCH := rmin 4520* = 'L' or 'l', DLAMCH := emax 4521* = 'O' or 'o', DLAMCH := rmax 4522* 4523* where 4524* 4525* eps = relative machine precision 4526* sfmin = safe minimum, such that 1/sfmin does not overflow 4527* base = base of the machine 4528* prec = eps*base 4529* t = number of (base) digits in the mantissa 4530* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise 4531* emin = minimum exponent before (gradual) underflow 4532* rmin = underflow threshold - base**(emin-1) 4533* emax = largest exponent before overflow 4534* rmax = overflow threshold - (base**emax)*(1-eps) 4535* 4536* ===================================================================== 4537* 4538* .. Parameters .. 4539 DOUBLE PRECISION ONE, ZERO 4540 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 4541* .. 4542* .. Local Scalars .. 4543 LOGICAL FIRST, LRND 4544 INTEGER BETA, IMAX, IMIN, IT 4545 DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, 4546 $ RND, SFMIN, SMALL, T 4547* .. 4548* .. External Functions .. 4549 LOGICAL LSAME 4550 EXTERNAL LSAME 4551* .. 4552* .. External Subroutines .. 4553 EXTERNAL DLAMC2 4554* .. 4555* .. Save statement .. 4556 SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, 4557 $ EMAX, RMAX, PREC 4558* .. 4559* .. Data statements .. 4560 DATA FIRST / .TRUE. / 4561* .. 4562* .. Executable Statements .. 4563* 4564 IF( FIRST ) THEN 4565 FIRST = .FALSE. 4566 CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) 4567 BASE = BETA 4568 T = IT 4569 IF( LRND ) THEN 4570 RND = ONE 4571 EPS = ( BASE**( 1-IT ) ) / 2 4572 ELSE 4573 RND = ZERO 4574 EPS = BASE**( 1-IT ) 4575 END IF 4576 PREC = EPS*BASE 4577 EMIN = IMIN 4578 EMAX = IMAX 4579 SFMIN = RMIN 4580 SMALL = ONE / RMAX 4581 IF( SMALL.GE.SFMIN ) THEN 4582* 4583* Use SMALL plus a bit, to avoid the possibility of rounding 4584* causing overflow when computing 1/sfmin. 4585* 4586 SFMIN = SMALL*( ONE+EPS ) 4587 END IF 4588 END IF 4589* 4590 IF( LSAME( CMACH, 'E' ) ) THEN 4591 RMACH = EPS 4592 ELSE IF( LSAME( CMACH, 'S' ) ) THEN 4593 RMACH = SFMIN 4594 ELSE IF( LSAME( CMACH, 'B' ) ) THEN 4595 RMACH = BASE 4596 ELSE IF( LSAME( CMACH, 'P' ) ) THEN 4597 RMACH = PREC 4598 ELSE IF( LSAME( CMACH, 'N' ) ) THEN 4599 RMACH = T 4600 ELSE IF( LSAME( CMACH, 'R' ) ) THEN 4601 RMACH = RND 4602 ELSE IF( LSAME( CMACH, 'M' ) ) THEN 4603 RMACH = EMIN 4604 ELSE IF( LSAME( CMACH, 'U' ) ) THEN 4605 RMACH = RMIN 4606 ELSE IF( LSAME( CMACH, 'L' ) ) THEN 4607 RMACH = EMAX 4608 ELSE IF( LSAME( CMACH, 'O' ) ) THEN 4609 RMACH = RMAX 4610 END IF 4611* 4612 DLAMCH = RMACH 4613 RETURN 4614* 4615* End of DLAMCH 4616* 4617 END 4618* 4619************************************************************************ 4620* 4621 SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) 4622* 4623* -- LAPACK auxiliary routine (version 1.1) -- 4624* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 4625* Courant Institute, Argonne National Lab, and Rice University 4626* October 31, 1992 4627* 4628* .. Scalar Arguments .. 4629 LOGICAL IEEE1, RND 4630 INTEGER BETA, T 4631* .. 4632* 4633* Purpose 4634* ======= 4635* 4636* DLAMC1 determines the machine parameters given by BETA, T, RND, and 4637* IEEE1. 4638* 4639* Arguments 4640* ========= 4641* 4642* BETA (output) INTEGER 4643* The base of the machine. 4644* 4645* T (output) INTEGER 4646* The number of ( BETA ) digits in the mantissa. 4647* 4648* RND (output) LOGICAL 4649* Specifies whether proper rounding ( RND = .TRUE. ) or 4650* chopping ( RND = .FALSE. ) occurs in addition. This may not 4651* be a reliable guide to the way in which the machine performs 4652* its arithmetic. 4653* 4654* IEEE1 (output) LOGICAL 4655* Specifies whether rounding appears to be done in the IEEE 4656* 'round to nearest' style. 4657* 4658* Further Details 4659* =============== 4660* 4661* The routine is based on the routine ENVRON by Malcolm and 4662* incorporates suggestions by Gentleman and Marovich. See 4663* 4664* Malcolm M. A. (1972) Algorithms to reveal properties of 4665* floating-point arithmetic. Comms. of the ACM, 15, 949-951. 4666* 4667* Gentleman W. M. and Marovich S. B. (1974) More on algorithms 4668* that reveal properties of floating point arithmetic units. 4669* Comms. of the ACM, 17, 276-277. 4670* 4671* ===================================================================== 4672* 4673* .. Local Scalars .. 4674 LOGICAL FIRST, LIEEE1, LRND 4675 INTEGER LBETA, LT 4676 DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 4677* .. 4678* .. External Functions .. 4679 DOUBLE PRECISION DLAMC3 4680 EXTERNAL DLAMC3 4681* .. 4682* .. Save statement .. 4683 SAVE FIRST, LIEEE1, LBETA, LRND, LT 4684* .. 4685* .. Data statements .. 4686 DATA FIRST / .TRUE. / 4687* .. 4688* .. Executable Statements .. 4689* 4690 IF( FIRST ) THEN 4691 FIRST = .FALSE. 4692 ONE = 1 4693* 4694* LBETA, LIEEE1, LT and LRND are the local values of BETA, 4695* IEEE1, T and RND. 4696* 4697* Throughout this routine we use the function DLAMC3 to ensure 4698* that relevant values are stored and not held in registers, or 4699* are not affected by optimizers. 4700* 4701* Compute a = 2.0**m with the smallest positive integer m such 4702* that 4703* 4704* fl( a + 1.0 ) = a. 4705* 4706 A = 1 4707 C = 1 4708* 4709*+ WHILE( C.EQ.ONE )LOOP 4710 10 CONTINUE 4711 IF( C.EQ.ONE ) THEN 4712 A = 2*A 4713 C = DLAMC3( A, ONE ) 4714 C = DLAMC3( C, -A ) 4715 GO TO 10 4716 END IF 4717*+ END WHILE 4718* 4719* Now compute b = 2.0**m with the smallest positive integer m 4720* such that 4721* 4722* fl( a + b ) .gt. a. 4723* 4724 B = 1 4725 C = DLAMC3( A, B ) 4726* 4727*+ WHILE( C.EQ.A )LOOP 4728 20 CONTINUE 4729 IF( C.EQ.A ) THEN 4730 B = 2*B 4731 C = DLAMC3( A, B ) 4732 GO TO 20 4733 END IF 4734*+ END WHILE 4735* 4736* Now compute the base. a and c are neighbouring floating point 4737* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so 4738* their difference is beta. Adding 0.25 to c is to ensure that it 4739* is truncated to beta and not ( beta - 1 ). 4740* 4741 QTR = ONE / 4 4742 SAVEC = C 4743 C = DLAMC3( C, -A ) 4744 LBETA = C + QTR 4745* 4746* Now determine whether rounding or chopping occurs, by adding a 4747* bit less than beta/2 and a bit more than beta/2 to a. 4748* 4749 B = LBETA 4750 F = DLAMC3( B / 2, -B / 100 ) 4751 C = DLAMC3( F, A ) 4752 IF( C.EQ.A ) THEN 4753 LRND = .TRUE. 4754 ELSE 4755 LRND = .FALSE. 4756 END IF 4757 F = DLAMC3( B / 2, B / 100 ) 4758 C = DLAMC3( F, A ) 4759 IF( ( LRND ) .AND. ( C.EQ.A ) ) 4760 $ LRND = .FALSE. 4761* 4762* Try and decide whether rounding is done in the IEEE 'round to 4763* nearest' style. B/2 is half a unit in the last place of the two 4764* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit 4765* zero, and SAVEC is odd. Thus adding B/2 to A should not change 4766* A, but adding B/2 to SAVEC should change SAVEC. 4767* 4768 T1 = DLAMC3( B / 2, A ) 4769 T2 = DLAMC3( B / 2, SAVEC ) 4770 LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND 4771* 4772* Now find the mantissa, t. It should be the integer part of 4773* log to the base beta of a, however it is safer to determine t 4774* by powering. So we find t as the smallest positive integer for 4775* which 4776* 4777* fl( beta**t + 1.0 ) = 1.0. 4778* 4779 LT = 0 4780 A = 1 4781 C = 1 4782* 4783*+ WHILE( C.EQ.ONE )LOOP 4784 30 CONTINUE 4785 IF( C.EQ.ONE ) THEN 4786 LT = LT + 1 4787 A = A*LBETA 4788 C = DLAMC3( A, ONE ) 4789 C = DLAMC3( C, -A ) 4790 GO TO 30 4791 END IF 4792*+ END WHILE 4793* 4794 END IF 4795* 4796 BETA = LBETA 4797 T = LT 4798 RND = LRND 4799 IEEE1 = LIEEE1 4800 RETURN 4801* 4802* End of DLAMC1 4803* 4804 END 4805* 4806************************************************************************ 4807* 4808 SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) 4809* 4810* -- LAPACK auxiliary routine (version 1.1) -- 4811* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 4812* Courant Institute, Argonne National Lab, and Rice University 4813* October 31, 1992 4814* 4815* .. Scalar Arguments .. 4816 LOGICAL RND 4817 INTEGER BETA, EMAX, EMIN, T 4818 DOUBLE PRECISION EPS, RMAX, RMIN 4819* .. 4820* 4821* Purpose 4822* ======= 4823* 4824* DLAMC2 determines the machine parameters specified in its argument 4825* list. 4826* 4827* Arguments 4828* ========= 4829* 4830* BETA (output) INTEGER 4831* The base of the machine. 4832* 4833* T (output) INTEGER 4834* The number of ( BETA ) digits in the mantissa. 4835* 4836* RND (output) LOGICAL 4837* Specifies whether proper rounding ( RND = .TRUE. ) or 4838* chopping ( RND = .FALSE. ) occurs in addition. This may not 4839* be a reliable guide to the way in which the machine performs 4840* its arithmetic. 4841* 4842* EPS (output) DOUBLE PRECISION 4843* The smallest positive number such that 4844* 4845* fl( 1.0 - EPS ) .LT. 1.0, 4846* 4847* where fl denotes the computed value. 4848* 4849* EMIN (output) INTEGER 4850* The minimum exponent before (gradual) underflow occurs. 4851* 4852* RMIN (output) DOUBLE PRECISION 4853* The smallest normalized number for the machine, given by 4854* BASE**( EMIN - 1 ), where BASE is the floating point value 4855* of BETA. 4856* 4857* EMAX (output) INTEGER 4858* The maximum exponent before overflow occurs. 4859* 4860* RMAX (output) DOUBLE PRECISION 4861* The largest positive number for the machine, given by 4862* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point 4863* value of BETA. 4864* 4865* Further Details 4866* =============== 4867* 4868* The computation of EPS is based on a routine PARANOIA by 4869* W. Kahan of the University of California at Berkeley. 4870* 4871* ===================================================================== 4872* 4873* .. Local Scalars .. 4874 LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND 4875 INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, 4876 $ NGNMIN, NGPMIN 4877 DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, 4878 $ SIXTH, SMALL, THIRD, TWO, ZERO 4879* .. 4880* .. External Functions .. 4881 DOUBLE PRECISION DLAMC3 4882 EXTERNAL DLAMC3 4883* .. 4884* .. External Subroutines .. 4885 EXTERNAL DLAMC1, DLAMC4, DLAMC5 4886* .. 4887* .. Intrinsic Functions .. 4888 INTRINSIC ABS, MAX, MIN 4889* .. 4890* .. Save statement .. 4891 SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, 4892 $ LRMIN, LT 4893* .. 4894* .. Data statements .. 4895 DATA FIRST / .TRUE. / , IWARN / .FALSE. / 4896* .. 4897* .. Executable Statements .. 4898* 4899 IF( FIRST ) THEN 4900 FIRST = .FALSE. 4901 ZERO = 0 4902 ONE = 1 4903 TWO = 2 4904* 4905* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of 4906* BETA, T, RND, EPS, EMIN and RMIN. 4907* 4908* Throughout this routine we use the function DLAMC3 to ensure 4909* that relevant values are stored and not held in registers, or 4910* are not affected by optimizers. 4911* 4912* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. 4913* 4914 CALL DLAMC1( LBETA, LT, LRND, LIEEE1 ) 4915* 4916* Start to find EPS. 4917* 4918 B = LBETA 4919 A = B**( -LT ) 4920 LEPS = A 4921* 4922* Try some tricks to see whether or not this is the correct EPS. 4923* 4924 B = TWO / 3 4925 HALF = ONE / 2 4926 SIXTH = DLAMC3( B, -HALF ) 4927 THIRD = DLAMC3( SIXTH, SIXTH ) 4928 B = DLAMC3( THIRD, -HALF ) 4929 B = DLAMC3( B, SIXTH ) 4930 B = ABS( B ) 4931 IF( B.LT.LEPS ) 4932 $ B = LEPS 4933* 4934 LEPS = 1 4935* 4936*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP 4937 10 CONTINUE 4938 IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN 4939 LEPS = B 4940 C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) 4941 C = DLAMC3( HALF, -C ) 4942 B = DLAMC3( HALF, C ) 4943 C = DLAMC3( HALF, -B ) 4944 B = DLAMC3( HALF, C ) 4945 GO TO 10 4946 END IF 4947*+ END WHILE 4948* 4949 IF( A.LT.LEPS ) 4950 $ LEPS = A 4951* 4952* Computation of EPS complete. 4953* 4954* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). 4955* Keep dividing A by BETA until (gradual) underflow occurs. This 4956* is detected when we cannot recover the previous A. 4957* 4958 RBASE = ONE / LBETA 4959 SMALL = ONE 4960 DO 20 I = 1, 3 4961 SMALL = DLAMC3( SMALL*RBASE, ZERO ) 4962 20 CONTINUE 4963 A = DLAMC3( ONE, SMALL ) 4964 CALL DLAMC4( NGPMIN, ONE, LBETA ) 4965 CALL DLAMC4( NGNMIN, -ONE, LBETA ) 4966 CALL DLAMC4( GPMIN, A, LBETA ) 4967 CALL DLAMC4( GNMIN, -A, LBETA ) 4968 IEEE = .FALSE. 4969* 4970 IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN 4971 IF( NGPMIN.EQ.GPMIN ) THEN 4972 LEMIN = NGPMIN 4973* ( Non twos-complement machines, no gradual underflow; 4974* e.g., VAX ) 4975 ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN 4976 LEMIN = NGPMIN - 1 + LT 4977 IEEE = .TRUE. 4978* ( Non twos-complement machines, with gradual underflow; 4979* e.g., IEEE standard followers ) 4980 ELSE 4981 LEMIN = MIN( NGPMIN, GPMIN ) 4982* ( A guess; no known machine ) 4983 IWARN = .TRUE. 4984 END IF 4985* 4986 ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN 4987 IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN 4988 LEMIN = MAX( NGPMIN, NGNMIN ) 4989* ( Twos-complement machines, no gradual underflow; 4990* e.g., CYBER 205 ) 4991 ELSE 4992 LEMIN = MIN( NGPMIN, NGNMIN ) 4993* ( A guess; no known machine ) 4994 IWARN = .TRUE. 4995 END IF 4996* 4997 ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. 4998 $ ( GPMIN.EQ.GNMIN ) ) THEN 4999 IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN 5000 LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT 5001* ( Twos-complement machines with gradual underflow; 5002* no known machine ) 5003 ELSE 5004 LEMIN = MIN( NGPMIN, NGNMIN ) 5005* ( A guess; no known machine ) 5006 IWARN = .TRUE. 5007 END IF 5008* 5009 ELSE 5010 LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) 5011* ( A guess; no known machine ) 5012 IWARN = .TRUE. 5013 END IF 5014*** 5015* Comment out this if block if EMIN is ok 5016 IF( IWARN ) THEN 5017 FIRST = .TRUE. 5018 WRITE( 6, FMT = 9999 )LEMIN 5019 END IF 5020*** 5021* 5022* Assume IEEE arithmetic if we found denormalised numbers above, 5023* or if arithmetic seems to round in the IEEE style, determined 5024* in routine DLAMC1. A true IEEE machine should have both things 5025* true; however, faulty machines may have one or the other. 5026* 5027 IEEE = IEEE .OR. LIEEE1 5028* 5029* Compute RMIN by successive division by BETA. We could compute 5030* RMIN as BASE**( EMIN - 1 ), but some machines underflow during 5031* this computation. 5032* 5033 LRMIN = 1 5034 DO 30 I = 1, 1 - LEMIN 5035 LRMIN = DLAMC3( LRMIN*RBASE, ZERO ) 5036 30 CONTINUE 5037* 5038* Finally, call DLAMC5 to compute EMAX and RMAX. 5039* 5040 CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) 5041 END IF 5042* 5043 BETA = LBETA 5044 T = LT 5045 RND = LRND 5046 EPS = LEPS 5047 EMIN = LEMIN 5048 RMIN = LRMIN 5049 EMAX = LEMAX 5050 RMAX = LRMAX 5051* 5052 RETURN 5053* 5054 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', 5055 $ ' EMIN = ', I8, / 5056 $ ' If, after inspection, the value EMIN looks', 5057 $ ' acceptable please comment out ', 5058 $ / ' the IF block as marked within the code of routine', 5059 $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / ) 5060* 5061* End of DLAMC2 5062* 5063 END 5064* 5065************************************************************************ 5066* 5067 DOUBLE PRECISION FUNCTION DLAMC3( A, B ) 5068* 5069* -- LAPACK auxiliary routine (version 1.1) -- 5070* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5071* Courant Institute, Argonne National Lab, and Rice University 5072* October 31, 1992 5073* 5074* .. Scalar Arguments .. 5075 DOUBLE PRECISION A, B 5076* .. 5077* 5078* Purpose 5079* ======= 5080* 5081* DLAMC3 is intended to force A and B to be stored prior to doing 5082* the addition of A and B , for use in situations where optimizers 5083* might hold one of these in a register. 5084* 5085* Arguments 5086* ========= 5087* 5088* A, B (input) DOUBLE PRECISION 5089* The values A and B. 5090* 5091* ===================================================================== 5092* 5093* .. Executable Statements .. 5094* 5095 DLAMC3 = A + B 5096* 5097 RETURN 5098* 5099* End of DLAMC3 5100* 5101 END 5102* 5103************************************************************************ 5104* 5105 SUBROUTINE DLAMC4( EMIN, START, BASE ) 5106* 5107* -- LAPACK auxiliary routine (version 1.1) -- 5108* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5109* Courant Institute, Argonne National Lab, and Rice University 5110* October 31, 1992 5111* 5112* .. Scalar Arguments .. 5113 INTEGER BASE, EMIN 5114 DOUBLE PRECISION START 5115* .. 5116* 5117* Purpose 5118* ======= 5119* 5120* DLAMC4 is a service routine for DLAMC2. 5121* 5122* Arguments 5123* ========= 5124* 5125* EMIN (output) EMIN 5126* The minimum exponent before (gradual) underflow, computed by 5127* setting A = START and dividing by BASE until the previous A 5128* can not be recovered. 5129* 5130* START (input) DOUBLE PRECISION 5131* The starting point for determining EMIN. 5132* 5133* BASE (input) INTEGER 5134* The base of the machine. 5135* 5136* ===================================================================== 5137* 5138* .. Local Scalars .. 5139 INTEGER I 5140 DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO 5141* .. 5142* .. External Functions .. 5143 DOUBLE PRECISION DLAMC3 5144 EXTERNAL DLAMC3 5145* .. 5146* .. Executable Statements .. 5147* 5148 A = START 5149 ONE = 1 5150 RBASE = ONE / BASE 5151 ZERO = 0 5152 EMIN = 1 5153 B1 = DLAMC3( A*RBASE, ZERO ) 5154 C1 = A 5155 C2 = A 5156 D1 = A 5157 D2 = A 5158*+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. 5159* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP 5160 10 CONTINUE 5161 IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. 5162 $ ( D2.EQ.A ) ) THEN 5163 EMIN = EMIN - 1 5164 A = B1 5165 B1 = DLAMC3( A / BASE, ZERO ) 5166 C1 = DLAMC3( B1*BASE, ZERO ) 5167 D1 = ZERO 5168 DO 20 I = 1, BASE 5169 D1 = D1 + B1 5170 20 CONTINUE 5171 B2 = DLAMC3( A*RBASE, ZERO ) 5172 C2 = DLAMC3( B2 / RBASE, ZERO ) 5173 D2 = ZERO 5174 DO 30 I = 1, BASE 5175 D2 = D2 + B2 5176 30 CONTINUE 5177 GO TO 10 5178 END IF 5179*+ END WHILE 5180* 5181 RETURN 5182* 5183* End of DLAMC4 5184* 5185 END 5186* 5187************************************************************************ 5188* 5189 SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) 5190* 5191* -- LAPACK auxiliary routine (version 1.1) -- 5192* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5193* Courant Institute, Argonne National Lab, and Rice University 5194* October 31, 1992 5195* 5196* .. Scalar Arguments .. 5197 LOGICAL IEEE 5198 INTEGER BETA, EMAX, EMIN, P 5199 DOUBLE PRECISION RMAX 5200* .. 5201* 5202* Purpose 5203* ======= 5204* 5205* DLAMC5 attempts to compute RMAX, the largest machine floating-point 5206* number, without overflow. It assumes that EMAX + abs(EMIN) sum 5207* approximately to a power of 2. It will fail on machines where this 5208* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, 5209* EMAX = 28718). It will also fail if the value supplied for EMIN is 5210* too large (i.e. too close to zero), probably with overflow. 5211* 5212* Arguments 5213* ========= 5214* 5215* BETA (input) INTEGER 5216* The base of floating-point arithmetic. 5217* 5218* P (input) INTEGER 5219* The number of base BETA digits in the mantissa of a 5220* floating-point value. 5221* 5222* EMIN (input) INTEGER 5223* The minimum exponent before (gradual) underflow. 5224* 5225* IEEE (input) LOGICAL 5226* A logical flag specifying whether or not the arithmetic 5227* system is thought to comply with the IEEE standard. 5228* 5229* EMAX (output) INTEGER 5230* The largest exponent before overflow 5231* 5232* RMAX (output) DOUBLE PRECISION 5233* The largest machine floating-point number. 5234* 5235* ===================================================================== 5236* 5237* .. Parameters .. 5238 DOUBLE PRECISION ZERO, ONE 5239 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 5240* .. 5241* .. Local Scalars .. 5242 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP 5243 DOUBLE PRECISION OLDY, RECBAS, Y, Z 5244* .. 5245* .. External Functions .. 5246 DOUBLE PRECISION DLAMC3 5247 EXTERNAL DLAMC3 5248* .. 5249* .. Intrinsic Functions .. 5250 INTRINSIC MOD 5251* .. 5252* .. Executable Statements .. 5253* 5254* First compute LEXP and UEXP, two powers of 2 that bound 5255* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum 5256* approximately to the bound that is closest to abs(EMIN). 5257* (EMAX is the exponent of the required number RMAX). 5258* 5259 LEXP = 1 5260 EXBITS = 1 5261 10 CONTINUE 5262 TRY = LEXP*2 5263 IF( TRY.LE.( -EMIN ) ) THEN 5264 LEXP = TRY 5265 EXBITS = EXBITS + 1 5266 GO TO 10 5267 END IF 5268 IF( LEXP.EQ.-EMIN ) THEN 5269 UEXP = LEXP 5270 ELSE 5271 UEXP = TRY 5272 EXBITS = EXBITS + 1 5273 END IF 5274* 5275* Now -LEXP is less than or equal to EMIN, and -UEXP is greater 5276* than or equal to EMIN. EXBITS is the number of bits needed to 5277* store the exponent. 5278* 5279 IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN 5280 EXPSUM = 2*LEXP 5281 ELSE 5282 EXPSUM = 2*UEXP 5283 END IF 5284* 5285* EXPSUM is the exponent range, approximately equal to 5286* EMAX - EMIN + 1 . 5287* 5288 EMAX = EXPSUM + EMIN - 1 5289 NBITS = 1 + EXBITS + P 5290* 5291* NBITS is the total number of bits needed to store a 5292* floating-point number. 5293* 5294 IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN 5295* 5296* Either there are an odd number of bits used to store a 5297* floating-point number, which is unlikely, or some bits are 5298* not used in the representation of numbers, which is possible, 5299* (e.g. Cray machines) or the mantissa has an implicit bit, 5300* (e.g. IEEE machines, Dec Vax machines), which is perhaps the 5301* most likely. We have to assume the last alternative. 5302* If this is true, then we need to reduce EMAX by one because 5303* there must be some way of representing zero in an implicit-bit 5304* system. On machines like Cray, we are reducing EMAX by one 5305* unnecessarily. 5306* 5307 EMAX = EMAX - 1 5308 END IF 5309* 5310 IF( IEEE ) THEN 5311* 5312* Assume we are on an IEEE machine which reserves one exponent 5313* for infinity and NaN. 5314* 5315 EMAX = EMAX - 1 5316 END IF 5317* 5318* Now create RMAX, the largest machine number, which should 5319* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . 5320* 5321* First compute 1.0 - BETA**(-P), being careful that the 5322* result is less than 1.0 . 5323* 5324 RECBAS = ONE / BETA 5325 Z = BETA - ONE 5326 Y = ZERO 5327 DO 20 I = 1, P 5328 Z = Z*RECBAS 5329 IF( Y.LT.ONE ) 5330 $ OLDY = Y 5331 Y = DLAMC3( Y, Z ) 5332 20 CONTINUE 5333 IF( Y.GE.ONE ) 5334 $ Y = OLDY 5335* 5336* Now multiply by BETA**EMAX to get RMAX. 5337* 5338 DO 30 I = 1, EMAX 5339 Y = DLAMC3( Y*BETA, ZERO ) 5340 30 CONTINUE 5341* 5342 RMAX = Y 5343 RETURN 5344* 5345* End of DLAMC5 5346* 5347 END 5348 DOUBLE PRECISION FUNCTION DLAPY2( X, Y ) 5349* 5350* -- LAPACK auxiliary routine (version 1.1) -- 5351* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5352* Courant Institute, Argonne National Lab, and Rice University 5353* October 31, 1992 5354* 5355* .. Scalar Arguments .. 5356 DOUBLE PRECISION X, Y 5357* .. 5358* 5359* Purpose 5360* ======= 5361* 5362* DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary 5363* overflow. 5364* 5365* Arguments 5366* ========= 5367* 5368* X (input) DOUBLE PRECISION 5369* Y (input) DOUBLE PRECISION 5370* X and Y specify the values x and y. 5371* 5372* ===================================================================== 5373* 5374* .. Parameters .. 5375 DOUBLE PRECISION ZERO 5376 PARAMETER ( ZERO = 0.0D0 ) 5377 DOUBLE PRECISION ONE 5378 PARAMETER ( ONE = 1.0D0 ) 5379* .. 5380* .. Local Scalars .. 5381 DOUBLE PRECISION W, XABS, YABS, Z 5382* .. 5383* .. Intrinsic Functions .. 5384 INTRINSIC ABS, MAX, MIN, SQRT 5385* .. 5386* .. Executable Statements .. 5387* 5388 XABS = ABS( X ) 5389 YABS = ABS( Y ) 5390 W = MAX( XABS, YABS ) 5391 Z = MIN( XABS, YABS ) 5392 IF( Z.EQ.ZERO ) THEN 5393 DLAPY2 = W 5394 ELSE 5395 DLAPY2 = W*SQRT( ONE+( Z / W )**2 ) 5396 END IF 5397 RETURN 5398* 5399* End of DLAPY2 5400* 5401 END 5402 INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, 5403 $ N4 ) 5404* 5405* -- LAPACK auxiliary routine (preliminary version) -- 5406* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5407* Courant Institute, Argonne National Lab, and Rice University 5408* February 20, 1992 5409* 5410* .. Scalar Arguments .. 5411 CHARACTER*( * ) NAME, OPTS 5412 INTEGER ISPEC, N1, N2, N3, N4 5413* .. 5414* 5415* Purpose 5416* ======= 5417* 5418* ILAENV is called from the LAPACK routines to choose problem-dependent 5419* parameters for the local environment. See ISPEC for a description of 5420* the parameters. 5421* 5422* This version provides a set of parameters which should give good, 5423* but not optimal, performance on many of the currently available 5424* computers. Users are encouraged to modify this subroutine to set 5425* the tuning parameters for their particular machine using the option 5426* and problem size information in the arguments. 5427* 5428* This routine will not function correctly if it is converted to all 5429* lower case. Converting it to all upper case is allowed. 5430* 5431* Arguments 5432* ========= 5433* 5434* ISPEC (input) INTEGER 5435* Specifies the parameter to be returned as the value of 5436* ILAENV. 5437* = 1: the optimal blocksize; if this value is 1, an unblocked 5438* algorithm will give the best performance. 5439* = 2: the minimum block size for which the block routine 5440* should be used; if the usable block size is less than 5441* this value, an unblocked routine should be used. 5442* = 3: the crossover point (in a block routine, for N less 5443* than this value, an unblocked routine should be used) 5444* = 4: the number of shifts, used in the nonsymmetric 5445* eigenvalue routines 5446* = 5: the minimum column dimension for blocking to be used; 5447* rectangular blocks must have dimension at least k by m, 5448* where k is given by ILAENV(2,...) and m by ILAENV(5,...) 5449* = 6: the crossover point for the SVD (when reducing an m by n 5450* matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds 5451* this value, a QR factorization is used first to reduce 5452* the matrix to a triangular form.) 5453* = 7: the number of processors 5454* = 8: the crossover point for the multishift QR and QZ methods 5455* for nonsymmetric eigenvalue problems. 5456* 5457* NAME (input) CHARACTER*(*) 5458* The name of the calling subroutine, in either upper case or 5459* lower case. 5460* 5461* OPTS (input) CHARACTER*(*) 5462* The character options to the subroutine NAME, concatenated 5463* into a single character string. For example, UPLO = 'U', 5464* TRANS = 'T', and DIAG = 'N' for a triangular routine would 5465* be specified as OPTS = 'UTN'. 5466* 5467* N1 (input) INTEGER 5468* N2 (input) INTEGER 5469* N3 (input) INTEGER 5470* N4 (input) INTEGER 5471* Problem dimensions for the subroutine NAME; these may not all 5472* be required. 5473* 5474* (ILAENV) (output) INTEGER 5475* >= 0: the value of the parameter specified by ISPEC 5476* < 0: if ILAENV = -k, the k-th argument had an illegal value. 5477* 5478* Further Details 5479* =============== 5480* 5481* The following conventions have been used when calling ILAENV from the 5482* LAPACK routines: 5483* 1) OPTS is a concatenation of all of the character options to 5484* subroutine NAME, in the same order that they appear in the 5485* argument list for NAME, even if they are not used in determining 5486* the value of the parameter specified by ISPEC. 5487* 2) The problem dimensions N1, N2, N3, N4 are specified in the order 5488* that they appear in the argument list for NAME. N1 is used 5489* first, N2 second, and so on, and unused problem dimensions are 5490* passed a value of -1. 5491* 3) The parameter value returned by ILAENV is checked for validity in 5492* the calling subroutine. For example, ILAENV is used to retrieve 5493* the optimal blocksize for STRTRI as follows: 5494* 5495* NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) 5496* IF( NB.LE.1 ) NB = MAX( 1, N ) 5497* 5498* ===================================================================== 5499* 5500* .. Local Scalars .. 5501 LOGICAL CNAME, SNAME 5502 CHARACTER*1 C1 5503 CHARACTER*2 C2, C4 5504 CHARACTER*3 C3 5505 CHARACTER*6 SUBNAM 5506 INTEGER I, IC, IZ, NB, NBMIN, NX 5507* .. 5508* .. Intrinsic Functions .. 5509 INTRINSIC CHAR, ICHAR, INT, MIN, REAL 5510* .. 5511* .. Executable Statements .. 5512* 5513 GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC 5514* 5515* Invalid value for ISPEC 5516* 5517 ILAENV = -1 5518 RETURN 5519* 5520 100 CONTINUE 5521* 5522* Convert NAME to upper case if the first character is lower case. 5523* 5524 ILAENV = 1 5525 SUBNAM = NAME 5526 IC = ICHAR( SUBNAM( 1:1 ) ) 5527 IZ = ICHAR( 'Z' ) 5528 IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN 5529* 5530* ASCII character set 5531* 5532 IF( IC.GE.97 .AND. IC.LE.122 ) THEN 5533 SUBNAM( 1:1 ) = CHAR( IC-32 ) 5534 DO 10 I = 2, 6 5535 IC = ICHAR( SUBNAM( I:I ) ) 5536 IF( IC.GE.97 .AND. IC.LE.122 ) 5537 $ SUBNAM( I:I ) = CHAR( IC-32 ) 5538 10 CONTINUE 5539 END IF 5540* 5541 ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN 5542* 5543* EBCDIC character set 5544* 5545 IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. 5546 $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. 5547 $ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN 5548 SUBNAM( 1:1 ) = CHAR( IC+64 ) 5549 DO 20 I = 2, 6 5550 IC = ICHAR( SUBNAM( I:I ) ) 5551 IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. 5552 $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. 5553 $ ( IC.GE.162 .AND. IC.LE.169 ) ) 5554 $ SUBNAM( I:I ) = CHAR( IC+64 ) 5555 20 CONTINUE 5556 END IF 5557* 5558 ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN 5559* 5560* Prime machines: ASCII+128 5561* 5562 IF( IC.GE.225 .AND. IC.LE.250 ) THEN 5563 SUBNAM( 1:1 ) = CHAR( IC-32 ) 5564 DO 30 I = 2, 6 5565 IC = ICHAR( SUBNAM( I:I ) ) 5566 IF( IC.GE.225 .AND. IC.LE.250 ) 5567 $ SUBNAM( I:I ) = CHAR( IC-32 ) 5568 30 CONTINUE 5569 END IF 5570 END IF 5571* 5572 C1 = SUBNAM( 1:1 ) 5573 SNAME = C1.EQ.'S' .OR. C1.EQ.'D' 5574 CNAME = C1.EQ.'C' .OR. C1.EQ.'Z' 5575 IF( .NOT.( CNAME .OR. SNAME ) ) 5576 $ RETURN 5577 C2 = SUBNAM( 2:3 ) 5578 C3 = SUBNAM( 4:6 ) 5579 C4 = C3( 2:3 ) 5580* 5581 GO TO ( 110, 200, 300 ) ISPEC 5582* 5583 110 CONTINUE 5584* 5585* ISPEC = 1: block size 5586* 5587* In these examples, separate code is provided for setting NB for 5588* real and complex. We assume that NB will take the same value in 5589* single or double precision. 5590* 5591 NB = 1 5592* 5593 IF( C2.EQ.'GE' ) THEN 5594 IF( C3.EQ.'TRF' ) THEN 5595 IF( SNAME ) THEN 5596 NB = 64 5597 ELSE 5598 NB = 64 5599 END IF 5600 ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. 5601 $ C3.EQ.'QLF' ) THEN 5602 IF( SNAME ) THEN 5603 NB = 32 5604 ELSE 5605 NB = 32 5606 END IF 5607 ELSE IF( C3.EQ.'HRD' ) THEN 5608 IF( SNAME ) THEN 5609 NB = 32 5610 ELSE 5611 NB = 32 5612 END IF 5613 ELSE IF( C3.EQ.'BRD' ) THEN 5614 IF( SNAME ) THEN 5615 NB = 32 5616 ELSE 5617 NB = 32 5618 END IF 5619 ELSE IF( C3.EQ.'TRI' ) THEN 5620 IF( SNAME ) THEN 5621 NB = 64 5622 ELSE 5623 NB = 64 5624 END IF 5625 END IF 5626 ELSE IF( C2.EQ.'PO' ) THEN 5627 IF( C3.EQ.'TRF' ) THEN 5628 IF( SNAME ) THEN 5629 NB = 64 5630 ELSE 5631 NB = 64 5632 END IF 5633 END IF 5634 ELSE IF( C2.EQ.'SY' ) THEN 5635 IF( C3.EQ.'TRF' ) THEN 5636 IF( SNAME ) THEN 5637 NB = 64 5638 ELSE 5639 NB = 64 5640 END IF 5641 ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN 5642 NB = 1 5643 ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN 5644 NB = 64 5645 END IF 5646 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN 5647 IF( C3.EQ.'TRF' ) THEN 5648 NB = 64 5649 ELSE IF( C3.EQ.'TRD' ) THEN 5650 NB = 1 5651 ELSE IF( C3.EQ.'GST' ) THEN 5652 NB = 64 5653 END IF 5654 ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN 5655 IF( C3( 1:1 ).EQ.'G' ) THEN 5656 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 5657 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 5658 $ C4.EQ.'BR' ) THEN 5659 NB = 32 5660 END IF 5661 ELSE IF( C3( 1:1 ).EQ.'M' ) THEN 5662 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 5663 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 5664 $ C4.EQ.'BR' ) THEN 5665 NB = 32 5666 END IF 5667 END IF 5668 ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN 5669 IF( C3( 1:1 ).EQ.'G' ) THEN 5670 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 5671 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 5672 $ C4.EQ.'BR' ) THEN 5673 NB = 32 5674 END IF 5675 ELSE IF( C3( 1:1 ).EQ.'M' ) THEN 5676 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 5677 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 5678 $ C4.EQ.'BR' ) THEN 5679 NB = 32 5680 END IF 5681 END IF 5682 ELSE IF( C2.EQ.'GB' ) THEN 5683 IF( C3.EQ.'TRF' ) THEN 5684 IF( SNAME ) THEN 5685 IF( N4.LE.64 ) THEN 5686 NB = 1 5687 ELSE 5688 NB = 32 5689 END IF 5690 ELSE 5691 IF( N4.LE.64 ) THEN 5692 NB = 1 5693 ELSE 5694 NB = 32 5695 END IF 5696 END IF 5697 END IF 5698 ELSE IF( C2.EQ.'PB' ) THEN 5699 IF( C3.EQ.'TRF' ) THEN 5700 IF( SNAME ) THEN 5701 IF( N2.LE.64 ) THEN 5702 NB = 1 5703 ELSE 5704 NB = 32 5705 END IF 5706 ELSE 5707 IF( N2.LE.64 ) THEN 5708 NB = 1 5709 ELSE 5710 NB = 32 5711 END IF 5712 END IF 5713 END IF 5714 ELSE IF( C2.EQ.'TR' ) THEN 5715 IF( C3.EQ.'TRI' ) THEN 5716 IF( SNAME ) THEN 5717 NB = 64 5718 ELSE 5719 NB = 64 5720 END IF 5721 END IF 5722 ELSE IF( C2.EQ.'LA' ) THEN 5723 IF( C3.EQ.'UUM' ) THEN 5724 IF( SNAME ) THEN 5725 NB = 64 5726 ELSE 5727 NB = 64 5728 END IF 5729 END IF 5730 ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN 5731 IF( C3.EQ.'EBZ' ) THEN 5732 NB = 1 5733 END IF 5734 END IF 5735 ILAENV = NB 5736 RETURN 5737* 5738 200 CONTINUE 5739* 5740* ISPEC = 2: minimum block size 5741* 5742 NBMIN = 2 5743 IF( C2.EQ.'GE' ) THEN 5744 IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. 5745 $ C3.EQ.'QLF' ) THEN 5746 IF( SNAME ) THEN 5747 NBMIN = 2 5748 ELSE 5749 NBMIN = 2 5750 END IF 5751 ELSE IF( C3.EQ.'HRD' ) THEN 5752 IF( SNAME ) THEN 5753 NBMIN = 2 5754 ELSE 5755 NBMIN = 2 5756 END IF 5757 ELSE IF( C3.EQ.'BRD' ) THEN 5758 IF( SNAME ) THEN 5759 NBMIN = 2 5760 ELSE 5761 NBMIN = 2 5762 END IF 5763 ELSE IF( C3.EQ.'TRI' ) THEN 5764 IF( SNAME ) THEN 5765 NBMIN = 2 5766 ELSE 5767 NBMIN = 2 5768 END IF 5769 END IF 5770 ELSE IF( C2.EQ.'SY' ) THEN 5771 IF( C3.EQ.'TRF' ) THEN 5772 IF( SNAME ) THEN 5773 NBMIN = 2 5774 ELSE 5775 NBMIN = 2 5776 END IF 5777 ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN 5778 NBMIN = 2 5779 END IF 5780 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN 5781 IF( C3.EQ.'TRD' ) THEN 5782 NBMIN = 2 5783 END IF 5784 ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN 5785 IF( C3( 1:1 ).EQ.'G' ) THEN 5786 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 5787 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 5788 $ C4.EQ.'BR' ) THEN 5789 NBMIN = 2 5790 END IF 5791 ELSE IF( C3( 1:1 ).EQ.'M' ) THEN 5792 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 5793 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 5794 $ C4.EQ.'BR' ) THEN 5795 NBMIN = 2 5796 END IF 5797 END IF 5798 ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN 5799 IF( C3( 1:1 ).EQ.'G' ) THEN 5800 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 5801 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 5802 $ C4.EQ.'BR' ) THEN 5803 NBMIN = 2 5804 END IF 5805 ELSE IF( C3( 1:1 ).EQ.'M' ) THEN 5806 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 5807 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 5808 $ C4.EQ.'BR' ) THEN 5809 NBMIN = 2 5810 END IF 5811 END IF 5812 END IF 5813 ILAENV = NBMIN 5814 RETURN 5815* 5816 300 CONTINUE 5817* 5818* ISPEC = 3: crossover point 5819* 5820 NX = 0 5821 IF( C2.EQ.'GE' ) THEN 5822 IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. 5823 $ C3.EQ.'QLF' ) THEN 5824 IF( SNAME ) THEN 5825 NX = 128 5826 ELSE 5827 NX = 128 5828 END IF 5829 ELSE IF( C3.EQ.'HRD' ) THEN 5830 IF( SNAME ) THEN 5831 NX = 128 5832 ELSE 5833 NX = 128 5834 END IF 5835 ELSE IF( C3.EQ.'BRD' ) THEN 5836 IF( SNAME ) THEN 5837 NX = 128 5838 ELSE 5839 NX = 128 5840 END IF 5841 END IF 5842 ELSE IF( C2.EQ.'SY' ) THEN 5843 IF( SNAME .AND. C3.EQ.'TRD' ) THEN 5844 NX = 1 5845 END IF 5846 ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN 5847 IF( C3.EQ.'TRD' ) THEN 5848 NX = 1 5849 END IF 5850 ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN 5851 IF( C3( 1:1 ).EQ.'G' ) THEN 5852 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 5853 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 5854 $ C4.EQ.'BR' ) THEN 5855 NX = 128 5856 END IF 5857 END IF 5858 ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN 5859 IF( C3( 1:1 ).EQ.'G' ) THEN 5860 IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. 5861 $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. 5862 $ C4.EQ.'BR' ) THEN 5863 NX = 128 5864 END IF 5865 END IF 5866 END IF 5867 ILAENV = NX 5868 RETURN 5869* 5870 400 CONTINUE 5871* 5872* ISPEC = 4: number of shifts (used by xHSEQR) 5873* 5874 ILAENV = 6 5875 RETURN 5876* 5877 500 CONTINUE 5878* 5879* ISPEC = 5: minimum column dimension (not used) 5880* 5881 ILAENV = 2 5882 RETURN 5883* 5884 600 CONTINUE 5885* 5886* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) 5887* 5888 ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) 5889 RETURN 5890* 5891 700 CONTINUE 5892* 5893* ISPEC = 7: number of processors (not used) 5894* 5895 ILAENV = 1 5896 RETURN 5897* 5898 800 CONTINUE 5899* 5900* ISPEC = 8: crossover point for multishift (used by xHSEQR) 5901* 5902 ILAENV = 50 5903 RETURN 5904* 5905* End of ILAENV 5906* 5907 END 5908 DOUBLE PRECISION FUNCTION DLANSY( NORM, UPLO, N, A, LDA, WORK ) 5909* 5910* -- LAPACK auxiliary routine (version 1.1) -- 5911* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5912* Courant Institute, Argonne National Lab, and Rice University 5913* October 31, 1992 5914* 5915* .. Scalar Arguments .. 5916 CHARACTER NORM, UPLO 5917 INTEGER LDA, N 5918* .. 5919* .. Array Arguments .. 5920 DOUBLE PRECISION A( LDA, * ), WORK( * ) 5921* .. 5922* 5923* Purpose 5924* ======= 5925* 5926* DLANSY returns the value of the one norm, or the Frobenius norm, or 5927* the infinity norm, or the element of largest absolute value of a 5928* real symmetric matrix A. 5929* 5930* Description 5931* =========== 5932* 5933* DLANSY returns the value 5934* 5935* DLANSY = ( max(abs(A(i,j))), NORM = 'M' or 'm' 5936* ( 5937* ( norm1(A), NORM = '1', 'O' or 'o' 5938* ( 5939* ( normI(A), NORM = 'I' or 'i' 5940* ( 5941* ( normF(A), NORM = 'F', 'f', 'E' or 'e' 5942* 5943* where norm1 denotes the one norm of a matrix (maximum column sum), 5944* normI denotes the infinity norm of a matrix (maximum row sum) and 5945* normF denotes the Frobenius norm of a matrix (square root of sum of 5946* squares). Note that max(abs(A(i,j))) is not a matrix norm. 5947* 5948* Arguments 5949* ========= 5950* 5951* NORM (input) CHARACTER*1 5952* Specifies the value to be returned in DLANSY as described 5953* above. 5954* 5955* UPLO (input) CHARACTER*1 5956* Specifies whether the upper or lower triangular part of the 5957* symmetric matrix A is to be referenced. 5958* = 'U': Upper triangular part of A is referenced 5959* = 'L': Lower triangular part of A is referenced 5960* 5961* N (input) INTEGER 5962* The order of the matrix A. N >= 0. When N = 0, DLANSY is 5963* set to zero. 5964* 5965* A (input) DOUBLE PRECISION array, dimension (LDA,N) 5966* The symmetric matrix A. If UPLO = 'U', the leading n by n 5967* upper triangular part of A contains the upper triangular part 5968* of the matrix A, and the strictly lower triangular part of A 5969* is not referenced. If UPLO = 'L', the leading n by n lower 5970* triangular part of A contains the lower triangular part of 5971* the matrix A, and the strictly upper triangular part of A is 5972* not referenced. 5973* 5974* LDA (input) INTEGER 5975* The leading dimension of the array A. LDA >= max(N,1). 5976* 5977* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK), 5978* where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, 5979* WORK is not referenced. 5980* 5981* ===================================================================== 5982* 5983* .. Parameters .. 5984 DOUBLE PRECISION ONE, ZERO 5985 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 5986* .. 5987* .. Local Scalars .. 5988 INTEGER I, J 5989 DOUBLE PRECISION ABSA, SCALE, SUM, VALUE 5990* .. 5991* .. External Subroutines .. 5992 EXTERNAL DLASSQ 5993* .. 5994* .. External Functions .. 5995 LOGICAL LSAME 5996 EXTERNAL LSAME 5997* .. 5998* .. Intrinsic Functions .. 5999 INTRINSIC ABS, MAX, SQRT 6000* .. 6001* .. Executable Statements .. 6002* 6003 IF( N.EQ.0 ) THEN 6004 VALUE = ZERO 6005 ELSE IF( LSAME( NORM, 'M' ) ) THEN 6006* 6007* Find max(abs(A(i,j))). 6008* 6009 VALUE = ZERO 6010 IF( LSAME( UPLO, 'U' ) ) THEN 6011 DO 20 J = 1, N 6012 DO 10 I = 1, J 6013 VALUE = MAX( VALUE, ABS( A( I, J ) ) ) 6014 10 CONTINUE 6015 20 CONTINUE 6016 ELSE 6017 DO 40 J = 1, N 6018 DO 30 I = J, N 6019 VALUE = MAX( VALUE, ABS( A( I, J ) ) ) 6020 30 CONTINUE 6021 40 CONTINUE 6022 END IF 6023 ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR. 6024 $ ( NORM.EQ.'1' ) ) THEN 6025* 6026* Find normI(A) ( = norm1(A), since A is symmetric). 6027* 6028 VALUE = ZERO 6029 IF( LSAME( UPLO, 'U' ) ) THEN 6030 DO 60 J = 1, N 6031 SUM = ZERO 6032 DO 50 I = 1, J - 1 6033 ABSA = ABS( A( I, J ) ) 6034 SUM = SUM + ABSA 6035 WORK( I ) = WORK( I ) + ABSA 6036 50 CONTINUE 6037 WORK( J ) = SUM + ABS( A( J, J ) ) 6038 60 CONTINUE 6039 DO 70 I = 1, N 6040 VALUE = MAX( VALUE, WORK( I ) ) 6041 70 CONTINUE 6042 ELSE 6043 DO 80 I = 1, N 6044 WORK( I ) = ZERO 6045 80 CONTINUE 6046 DO 100 J = 1, N 6047 SUM = WORK( J ) + ABS( A( J, J ) ) 6048 DO 90 I = J + 1, N 6049 ABSA = ABS( A( I, J ) ) 6050 SUM = SUM + ABSA 6051 WORK( I ) = WORK( I ) + ABSA 6052 90 CONTINUE 6053 VALUE = MAX( VALUE, SUM ) 6054 100 CONTINUE 6055 END IF 6056 ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN 6057* 6058* Find normF(A). 6059* 6060 SCALE = ZERO 6061 SUM = ONE 6062 IF( LSAME( UPLO, 'U' ) ) THEN 6063 DO 110 J = 2, N 6064 CALL DLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) 6065 110 CONTINUE 6066 ELSE 6067 DO 120 J = 1, N - 1 6068 CALL DLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) 6069 120 CONTINUE 6070 END IF 6071 SUM = 2*SUM 6072 CALL DLASSQ( N, A, LDA+1, SCALE, SUM ) 6073 VALUE = SCALE*SQRT( SUM ) 6074 END IF 6075* 6076 DLANSY = VALUE 6077 RETURN 6078* 6079* End of DLANSY 6080* 6081 END 6082 SUBROUTINE DLASSQ( N, X, INCX, SCALE, SUMSQ ) 6083* 6084* -- LAPACK auxiliary routine (version 1.1) -- 6085* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6086* Courant Institute, Argonne National Lab, and Rice University 6087* October 31, 1992 6088* 6089* .. Scalar Arguments .. 6090 INTEGER INCX, N 6091 DOUBLE PRECISION SCALE, SUMSQ 6092* .. 6093* .. Array Arguments .. 6094 DOUBLE PRECISION X( * ) 6095* .. 6096* 6097* Purpose 6098* ======= 6099* 6100* DLASSQ returns the values scl and smsq such that 6101* 6102* ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, 6103* 6104* where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is 6105* assumed to be non-negative and scl returns the value 6106* 6107* scl = max( scale, abs( x( i ) ) ). 6108* 6109* scale and sumsq must be supplied in SCALE and SUMSQ and 6110* scl and smsq are overwritten on SCALE and SUMSQ respectively. 6111* 6112* The routine makes only one pass through the vector x. 6113* 6114* Arguments 6115* ========= 6116* 6117* N (input) INTEGER 6118* The number of elements to be used from the vector X. 6119* 6120* X (input) DOUBLE PRECISION 6121* The vector for which a scaled sum of squares is computed. 6122* x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. 6123* 6124* INCX (input) INTEGER 6125* The increment between successive values of the vector X. 6126* INCX > 0. 6127* 6128* SCALE (input/output) DOUBLE PRECISION 6129* On entry, the value scale in the equation above. 6130* On exit, SCALE is overwritten with scl , the scaling factor 6131* for the sum of squares. 6132* 6133* SUMSQ (input/output) DOUBLE PRECISION 6134* On entry, the value sumsq in the equation above. 6135* On exit, SUMSQ is overwritten with smsq , the basic sum of 6136* squares from which scl has been factored out. 6137* 6138* ===================================================================== 6139* 6140* .. Parameters .. 6141 DOUBLE PRECISION ZERO 6142 PARAMETER ( ZERO = 0.0D+0 ) 6143* .. 6144* .. Local Scalars .. 6145 INTEGER IX 6146 DOUBLE PRECISION ABSXI 6147* .. 6148* .. Intrinsic Functions .. 6149 INTRINSIC ABS 6150* .. 6151* .. Executable Statements .. 6152* 6153 IF( N.GT.0 ) THEN 6154 DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX 6155 IF( X( IX ).NE.ZERO ) THEN 6156 ABSXI = ABS( X( IX ) ) 6157 IF( SCALE.LT.ABSXI ) THEN 6158 SUMSQ = 1 + SUMSQ*( SCALE / ABSXI )**2 6159 SCALE = ABSXI 6160 ELSE 6161 SUMSQ = SUMSQ + ( ABSXI / SCALE )**2 6162 END IF 6163 END IF 6164 10 CONTINUE 6165 END IF 6166 RETURN 6167* 6168* End of DLASSQ 6169* 6170 END 6171 DOUBLE PRECISION FUNCTION DLARAN( ISEED ) 6172* 6173* -- LAPACK auxiliary routine (version 1.1) -- 6174* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6175* Courant Institute, Argonne National Lab, and Rice University 6176* February 29, 1992 6177* 6178* .. Array Arguments .. 6179 INTEGER ISEED( 4 ) 6180* .. 6181* 6182* Purpose 6183* ======= 6184* 6185* DLARAN returns a random DOUBLE PRECISION number from a uniform (0,1) 6186* distribution. 6187* 6188* Arguments 6189* ========= 6190* 6191* ISEED (input/output) INTEGER array, dimension (4) 6192* On entry, the seed of the random number generator; the array 6193* elements must be between 0 and 4095, and ISEED(4) must be 6194* odd. 6195* On exit, the seed is updated. 6196* 6197* Further Details 6198* =============== 6199* 6200* This routine uses a multiplicative congruential method with modulus 6201* 2**48 and multiplier 33952834046453 (see G.S.Fishman, 6202* 'Multiplicative congruential random number generators with modulus 6203* 2**b: an exhaustive analysis for b = 32 and a partial analysis for 6204* b = 48', Math. Comp. 189, pp 331-344, 1990). 6205* 6206* 48-bit integers are stored in 4 integer array elements with 12 bits 6207* per element. Hence the routine is portable across machines with 6208* integers of 32 bits or more. 6209* 6210* ===================================================================== 6211* 6212* .. Parameters .. 6213 INTEGER M1, M2, M3, M4 6214 PARAMETER ( M1 = 494, M2 = 322, M3 = 2508, M4 = 2549 ) 6215 DOUBLE PRECISION ONE 6216 PARAMETER ( ONE = 1.0D+0 ) 6217 INTEGER IPW2 6218 DOUBLE PRECISION R 6219 PARAMETER ( IPW2 = 4096, R = ONE / IPW2 ) 6220* .. 6221* .. Local Scalars .. 6222 INTEGER IT1, IT2, IT3, IT4 6223* .. 6224* .. Intrinsic Functions .. 6225 INTRINSIC DBLE, MOD 6226* .. 6227* .. Executable Statements .. 6228* 6229* multiply the seed by the multiplier modulo 2**48 6230* 6231 IT4 = ISEED( 4 )*M4 6232 IT3 = IT4 / IPW2 6233 IT4 = IT4 - IPW2*IT3 6234 IT3 = IT3 + ISEED( 3 )*M4 + ISEED( 4 )*M3 6235 IT2 = IT3 / IPW2 6236 IT3 = IT3 - IPW2*IT2 6237 IT2 = IT2 + ISEED( 2 )*M4 + ISEED( 3 )*M3 + ISEED( 4 )*M2 6238 IT1 = IT2 / IPW2 6239 IT2 = IT2 - IPW2*IT1 6240 IT1 = IT1 + ISEED( 1 )*M4 + ISEED( 2 )*M3 + ISEED( 3 )*M2 + 6241 $ ISEED( 4 )*M1 6242 IT1 = MOD( IT1, IPW2 ) 6243* 6244* return updated seed 6245* 6246 ISEED( 1 ) = IT1 6247 ISEED( 2 ) = IT2 6248 ISEED( 3 ) = IT3 6249 ISEED( 4 ) = IT4 6250* 6251* convert 48-bit integer to a DOUBLE PRECISION number in the interval (0,1) 6252* 6253 DLARAN = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R* 6254 $ ( DBLE( IT4 ) ) ) ) ) 6255 RETURN 6256* 6257* End of DLARAN 6258* 6259 END 6260 LOGICAL FUNCTION LSAME( CA, CB ) 6261* 6262* -- LAPACK auxiliary routine (version 1.1) -- 6263* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6264* Courant Institute, Argonne National Lab, and Rice University 6265* February 29, 1992 6266* 6267* .. Scalar Arguments .. 6268 CHARACTER CA, CB 6269* .. 6270* 6271* Purpose 6272* ======= 6273* 6274* LSAME returns .TRUE. if CA is the same letter as CB regardless of 6275* case. 6276* 6277* Arguments 6278* ========= 6279* 6280* CA (input) CHARACTER*1 6281* CB (input) CHARACTER*1 6282* CA and CB specify the single characters to be compared. 6283* 6284* .. Intrinsic Functions .. 6285 INTRINSIC ICHAR 6286* .. 6287* .. Local Scalars .. 6288 INTEGER INTA, INTB, ZCODE 6289* .. 6290* .. Executable Statements .. 6291* 6292* Test if the characters are equal 6293* 6294 LSAME = CA.EQ.CB 6295 IF( LSAME ) 6296 $ RETURN 6297* 6298* Now test for equivalence if both characters are alphabetic. 6299* 6300 ZCODE = ICHAR( 'Z' ) 6301* 6302* Use 'Z' rather than 'A' so that ASCII can be detected on Prime 6303* machines, on which ICHAR returns a value with bit 8 set. 6304* ICHAR('A') on Prime machines returns 193 which is the same as 6305* ICHAR('A') on an EBCDIC machine. 6306* 6307 INTA = ICHAR( CA ) 6308 INTB = ICHAR( CB ) 6309* 6310 IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN 6311* 6312* ASCII is assumed - ZCODE is the ASCII code of either lower or 6313* upper case 'Z'. 6314* 6315 IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 6316 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 6317* 6318 ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN 6319* 6320* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or 6321* upper case 'Z'. 6322* 6323 IF( INTA.GE.129 .AND. INTA.LE.137 .OR. 6324 $ INTA.GE.145 .AND. INTA.LE.153 .OR. 6325 $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 6326 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. 6327 $ INTB.GE.145 .AND. INTB.LE.153 .OR. 6328 $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 6329* 6330 ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN 6331* 6332* ASCII is assumed, on Prime machines - ZCODE is the ASCII code 6333* plus 128 of either lower or upper case 'Z'. 6334* 6335 IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 6336 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 6337 END IF 6338 LSAME = INTA.EQ.INTB 6339* 6340* RETURN 6341* 6342* End of LSAME 6343* 6344 END 6345 LOGICAL FUNCTION LSAMEN( N, CA, CB ) 6346* 6347* -- LAPACK auxiliary routine (version 1.0b) -- 6348* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6349* Courant Institute, Argonne National Lab, and Rice University 6350* February 29, 1992 6351* 6352* .. Scalar Arguments .. 6353 CHARACTER*( * ) CA, CB 6354 INTEGER N 6355* .. 6356* 6357* Purpose 6358* ======= 6359* 6360* LSAMEN tests if the first N letters of CA are the same as the 6361* first N letters of CB, regardless of case. 6362* LSAMEN returns .TRUE. if CA and CB are equivalent except for case 6363* and .FALSE. otherwise. LSAMEN also returns .FALSE. if LEN( CA ) 6364* or LEN( CB ) is less than N. 6365* 6366* Arguments 6367* ========= 6368* 6369* N (input) INTEGER 6370* The number of characters in CA and CB to be compared. 6371* 6372* CA (input) CHARACTER*(*) 6373* CB (input) CHARACTER*(*) 6374* CA and CB specify two character strings of length at least N. 6375* Only the first N characters of each string will be accessed. 6376* 6377* .. Local Scalars .. 6378 INTEGER I 6379* .. 6380* .. External Functions .. 6381 LOGICAL LSAME 6382 EXTERNAL LSAME 6383* .. 6384* .. Intrinsic Functions .. 6385 INTRINSIC LEN 6386* .. 6387* .. Executable Statements .. 6388* 6389 LSAMEN = .FALSE. 6390 IF( LEN( CA ).LT.N .OR. LEN( CB ).LT.N ) 6391 $ GO TO 20 6392* 6393* Do for each character in the two strings. 6394* 6395 DO 10 I = 1, N 6396* 6397* Test if the characters are equal using LSAME. 6398* 6399 IF( .NOT.LSAME( CA( I: I ), CB( I: I ) ) ) 6400 $ GO TO 20 6401* 6402 10 CONTINUE 6403 LSAMEN = .TRUE. 6404* 6405 20 CONTINUE 6406 RETURN 6407* 6408* End of LSAMEN 6409* 6410 END 6411