1 SUBROUTINE CDIV(AR,AI,BR,BI,CR,CI) 2 DOUBLE PRECISION AR,AI,BR,BI,CR,CI 3C 4C COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI) 5C 6 DOUBLE PRECISION S,ARS,AIS,BRS,BIS 7 S = DABS(BR) + DABS(BI) 8 ARS = AR/S 9 AIS = AI/S 10 BRS = BR/S 11 BIS = BI/S 12 S = BRS**2 + BIS**2 13 CR = (ARS*BRS + AIS*BIS)/S 14 CI = (AIS*BRS - ARS*BIS)/S 15 RETURN 16 END 17 DOUBLE PRECISION FUNCTION EPSLON (X) 18 DOUBLE PRECISION X 19C 20C ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. 21C 22 DOUBLE PRECISION A,B,C,EPS 23C 24C THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS 25C SATISFYING THE FOLLOWING TWO ASSUMPTIONS, 26C 1. THE BASE USED IN REPRESENTING FLOATING POINT 27C NUMBERS IS NOT A POWER OF THREE. 28C 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO 29C THE ACCURACY USED IN FLOATING POINT VARIABLES 30C THAT ARE STORED IN MEMORY. 31C THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO 32C FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING 33C ASSUMPTION 2. 34C UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, 35C A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, 36C B HAS A ZERO FOR ITS LAST BIT OR DIGIT, 37C C IS NOT EXACTLY EQUAL TO ONE, 38C EPS MEASURES THE SEPARATION OF 1.0 FROM 39C THE NEXT LARGER FLOATING POINT NUMBER. 40C THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED 41C ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD. 42C 43C THIS VERSION DATED 4/6/83. 44C 45 A = 4.0D0/3.0D0 46 10 B = A - 1.0D0 47 C = B + B + B 48 EPS = DABS(C-1.0D0) 49 IF (EPS .EQ. 0.0D0) GO TO 10 50 EPSLON = EPS*DABS(X) 51 RETURN 52 END 53 SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR) 54C 55 INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR 56 DOUBLE PRECISION H(NM,N),WR(N),WI(N) 57 DOUBLE PRECISION P,Q,R,S,T,W,X,Y,ZZ,NORM,TST1,TST2 58 LOGICAL NOTLAS 59* 60* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 61* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 62* OPST IS USED TO ACCUMULATE SMALL CONTRIBUTIONS TO OPS 63* TO AVOID ROUNDOFF ERROR 64* .. COMMON BLOCKS .. 65 COMMON /LATIME/ OPS, ITCNT 66* .. 67* .. SCALARS IN COMMON .. 68 DOUBLE PRECISION OPS, ITCNT, OPST 69* .. 70C 71C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, 72C NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. 73C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). 74C 75C THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL 76C UPPER HESSENBERG MATRIX BY THE QR METHOD. 77C 78C ON INPUT 79C 80C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 81C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 82C DIMENSION STATEMENT. 83C 84C N IS THE ORDER OF THE MATRIX. 85C 86C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING 87C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, 88C SET LOW=1, IGH=N. 89C 90C H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT 91C THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG 92C FORM BY ELMHES OR ORTHES, IF PERFORMED, IS STORED 93C IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. 94C 95C ON OUTPUT 96C 97C H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED 98C BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND 99C BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED. 100C 101C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, 102C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES 103C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS 104C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE 105C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN 106C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT 107C FOR INDICES IERR+1,...,N. 108C 109C IERR IS SET TO 110C ZERO FOR NORMAL RETURN, 111C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED 112C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. 113C 114C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 115C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 116C 117C THIS VERSION DATED AUGUST 1983. 118C MODIFIED ON 11/1/89; ADJUSTING INDICES OF LOOPS 119C 200, 210, 230, AND 240 TO INCREASE PERFORMANCE. JACK DONGARRA 120C 121C ------------------------------------------------------------------ 122C 123* 124 EXTERNAL DLAMCH 125 DOUBLE PRECISION DLAMCH, UNFL,OVFL,ULP,SMLNUM,SMALL 126 IF (N.LE.0) RETURN 127* 128* 129* INITIALIZE 130 ITCNT = 0 131 OPST = 0 132 IERR = 0 133 K = 1 134C .......... STORE ROOTS ISOLATED BY BALANC 135C AND COMPUTE MATRIX NORM .......... 136 DO 50 I = 1, N 137 K = I 138 IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50 139 WR(I) = H(I,I) 140 WI(I) = 0.0D0 141 50 CONTINUE 142* 143* INCREMENT OPCOUNT FOR COMPUTING MATRIX NORM 144 OPS = OPS + (IGH-LOW+1)*(IGH-LOW+2)/2 145* 146* COMPUTE THE 1-NORM OF MATRIX H 147* 148 NORM = 0.0D0 149 DO 5 J = LOW, IGH 150 S = 0.0D0 151 DO 4 I = LOW, MIN(IGH,J+1) 152 S = S + DABS(H(I,J)) 153 4 CONTINUE 154 NORM = MAX(NORM, S) 155 5 CONTINUE 156* 157 UNFL = DLAMCH( 'SAFE MINIMUM' ) 158 OVFL = DLAMCH( 'OVERFLOW' ) 159 ULP = DLAMCH( 'EPSILON' )*DLAMCH( 'BASE' ) 160 SMLNUM = MAX( UNFL*( N / ULP ), N / ( ULP*OVFL ) ) 161 SMALL = MAX( SMLNUM, ULP*NORM ) 162C 163 EN = IGH 164 T = 0.0D0 165 ITN = 30*N 166C .......... SEARCH FOR NEXT EIGENVALUES .......... 167 60 IF (EN .LT. LOW) GO TO 1001 168 ITS = 0 169 NA = EN - 1 170 ENM2 = NA - 1 171C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT 172C FOR L=EN STEP -1 UNTIL LOW DO -- .......... 173* REPLACE SPLITTING CRITERION WITH NEW ONE AS IN LAPACK 174* 175 70 DO 80 LL = LOW, EN 176 L = EN + LOW - LL 177 IF (L .EQ. LOW) GO TO 100 178 S = DABS(H(L-1,L-1)) + DABS(H(L,L)) 179 IF (S .EQ. 0.0D0) S = NORM 180 IF (DABS(H(L,L-1)) .LE. MAX(ULP*S,SMALL)) GO TO 100 181 80 CONTINUE 182C .......... FORM SHIFT .......... 183 100 CONTINUE 184* 185* INCREMENT OP COUNT FOR CONVERGENCE TEST 186 OPS = OPS + 2*(EN-L+1) 187 X = H(EN,EN) 188 IF (L .EQ. EN) GO TO 270 189 Y = H(NA,NA) 190 W = H(EN,NA) * H(NA,EN) 191 IF (L .EQ. NA) GO TO 280 192 IF (ITN .EQ. 0) GO TO 1000 193 IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 194C .......... FORM EXCEPTIONAL SHIFT .......... 195* 196* INCREMENT OP COUNT FOR FORMING EXCEPTIONAL SHIFT 197 OPS = OPS + (EN-LOW+6) 198 T = T + X 199C 200 DO 120 I = LOW, EN 201 120 H(I,I) = H(I,I) - X 202C 203 S = DABS(H(EN,NA)) + DABS(H(NA,ENM2)) 204 X = 0.75D0 * S 205 Y = X 206 W = -0.4375D0 * S * S 207 130 ITS = ITS + 1 208 ITN = ITN - 1 209* 210* UPDATE ITERATION NUMBER 211 ITCNT = 30*N - ITN 212C .......... LOOK FOR TWO CONSECUTIVE SMALL 213C SUB-DIAGONAL ELEMENTS. 214C FOR M=EN-2 STEP -1 UNTIL L DO -- .......... 215* REPLACE SPLITTING CRITERION WITH NEW ONE AS IN LAPACK 216 DO 140 MM = L, ENM2 217 M = ENM2 + L - MM 218 ZZ = H(M,M) 219 R = X - ZZ 220 S = Y - ZZ 221 P = (R * S - W) / H(M+1,M) + H(M,M+1) 222 Q = H(M+1,M+1) - ZZ - R - S 223 R = H(M+2,M+1) 224 S = DABS(P) + DABS(Q) + DABS(R) 225 P = P / S 226 Q = Q / S 227 R = R / S 228 IF (M .EQ. L) GO TO 150 229 TST1 = DABS(P)*(DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1))) 230 TST2 = DABS(H(M,M-1))*(DABS(Q) + DABS(R)) 231 IF ( TST2 .LE. MAX(ULP*TST1,SMALL) ) GO TO 150 232 140 CONTINUE 233C 234 150 CONTINUE 235* 236* INCREMENT OPCOUNT FOR LOOP 140 237 OPST = OPST + 20*(ENM2-M+1) 238 MP2 = M + 2 239C 240 DO 160 I = MP2, EN 241 H(I,I-2) = 0.0D0 242 IF (I .EQ. MP2) GO TO 160 243 H(I,I-3) = 0.0D0 244 160 CONTINUE 245C .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND 246C COLUMNS M TO EN .......... 247* 248* INCREMENT OPCOUNT FOR LOOP 260 249 OPST = OPST + 18*(NA-M+1) 250 DO 260 K = M, NA 251 NOTLAS = K .NE. NA 252 IF (K .EQ. M) GO TO 170 253 P = H(K,K-1) 254 Q = H(K+1,K-1) 255 R = 0.0D0 256 IF (NOTLAS) R = H(K+2,K-1) 257 X = DABS(P) + DABS(Q) + DABS(R) 258 IF (X .EQ. 0.0D0) GO TO 260 259 P = P / X 260 Q = Q / X 261 R = R / X 262 170 S = DSIGN(DSQRT(P*P+Q*Q+R*R),P) 263 IF (K .EQ. M) GO TO 180 264 H(K,K-1) = -S * X 265 GO TO 190 266 180 IF (L .NE. M) H(K,K-1) = -H(K,K-1) 267 190 P = P + S 268 X = P / S 269 Y = Q / S 270 ZZ = R / S 271 Q = Q / P 272 R = R / P 273 IF (NOTLAS) GO TO 225 274C .......... ROW MODIFICATION .......... 275* 276* INCREMENT OPCOUNT 277 OPS = OPS + 6*(EN-K+1) 278 DO 200 J = K, EN 279 P = H(K,J) + Q * H(K+1,J) 280 H(K,J) = H(K,J) - P * X 281 H(K+1,J) = H(K+1,J) - P * Y 282 200 CONTINUE 283C 284 J = MIN0(EN,K+3) 285C .......... COLUMN MODIFICATION .......... 286* 287* INCREMENT OPCOUNT 288 OPS = OPS + 6*(J-L+1) 289 DO 210 I = L, J 290 P = X * H(I,K) + Y * H(I,K+1) 291 H(I,K) = H(I,K) - P 292 H(I,K+1) = H(I,K+1) - P * Q 293 210 CONTINUE 294 GO TO 255 295 225 CONTINUE 296C .......... ROW MODIFICATION .......... 297* 298* INCREMENT OPCOUNT 299 OPS = OPS + 10*(EN-K+1) 300 DO 230 J = K, EN 301 P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J) 302 H(K,J) = H(K,J) - P * X 303 H(K+1,J) = H(K+1,J) - P * Y 304 H(K+2,J) = H(K+2,J) - P * ZZ 305 230 CONTINUE 306C 307 J = MIN0(EN,K+3) 308C .......... COLUMN MODIFICATION .......... 309* 310* INCREMENT OPCOUNT 311 OPS = OPS + 10*(J-L+1) 312 DO 240 I = L, J 313 P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2) 314 H(I,K) = H(I,K) - P 315 H(I,K+1) = H(I,K+1) - P * Q 316 H(I,K+2) = H(I,K+2) - P * R 317 240 CONTINUE 318 255 CONTINUE 319C 320 260 CONTINUE 321C 322 GO TO 70 323C .......... ONE ROOT FOUND .......... 324 270 WR(EN) = X + T 325 WI(EN) = 0.0D0 326 EN = NA 327 GO TO 60 328C .......... TWO ROOTS FOUND .......... 329 280 P = (Y - X) / 2.0D0 330 Q = P * P + W 331 ZZ = DSQRT(DABS(Q)) 332 X = X + T 333* 334* INCREMENT OP COUNT FOR FINDING TWO ROOTS. 335 OPST = OPST + 8 336 IF (Q .LT. 0.0D0) GO TO 320 337C .......... REAL PAIR .......... 338 ZZ = P + DSIGN(ZZ,P) 339 WR(NA) = X + ZZ 340 WR(EN) = WR(NA) 341 IF (ZZ .NE. 0.0D0) WR(EN) = X - W / ZZ 342 WI(NA) = 0.0D0 343 WI(EN) = 0.0D0 344 GO TO 330 345C .......... COMPLEX PAIR .......... 346 320 WR(NA) = X + P 347 WR(EN) = X + P 348 WI(NA) = ZZ 349 WI(EN) = -ZZ 350 330 EN = ENM2 351 GO TO 60 352C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT 353C CONVERGED AFTER 30*N ITERATIONS .......... 354 1000 IERR = EN 355 1001 CONTINUE 356* 357* COMPUTE FINAL OP COUNT 358 OPS = OPS + OPST 359 RETURN 360 END 361 SUBROUTINE HQR2(NM,N,LOW,IGH,H,WR,WI,Z,IERR) 362C 363 INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN, 364 X IGH,ITN,ITS,LOW,MP2,ENM2,IERR 365 DOUBLE PRECISION H(NM,N),WR(N),WI(N),Z(NM,N) 366 DOUBLE PRECISION P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,TST1,TST2 367 LOGICAL NOTLAS 368* 369* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 370* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 371* OPST IS USED TO ACCUMULATE SMALL CONTRIBUTIONS TO OPS 372* TO AVOID ROUNDOFF ERROR 373* .. COMMON BLOCKS .. 374 COMMON /LATIME/ OPS, ITCNT 375* .. 376* .. SCALARS IN COMMON .. 377 DOUBLE PRECISION OPS, ITCNT, OPST 378* .. 379C 380C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2, 381C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. 382C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). 383C 384C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS 385C OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD. THE 386C EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND 387C IF ELMHES AND ELTRAN OR ORTHES AND ORTRAN HAVE 388C BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM 389C AND TO ACCUMULATE THE SIMILARITY TRANSFORMATIONS. 390C 391C ON INPUT 392C 393C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 394C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 395C DIMENSION STATEMENT. 396C 397C N IS THE ORDER OF THE MATRIX. 398C 399C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING 400C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, 401C SET LOW=1, IGH=N. 402C 403C H CONTAINS THE UPPER HESSENBERG MATRIX. 404C 405C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY ELTRAN 406C AFTER THE REDUCTION BY ELMHES, OR BY ORTRAN AFTER THE 407C REDUCTION BY ORTHES, IF PERFORMED. IF THE EIGENVECTORS 408C OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE 409C IDENTITY MATRIX. 410C 411C ON OUTPUT 412C 413C H HAS BEEN DESTROYED. 414C 415C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, 416C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES 417C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS 418C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE 419C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN 420C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT 421C FOR INDICES IERR+1,...,N. 422C 423C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. 424C IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z 425C CONTAINS ITS EIGENVECTOR. IF THE I-TH EIGENVALUE IS COMPLEX 426C WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH 427C COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS 428C EIGENVECTOR. THE EIGENVECTORS ARE UNNORMALIZED. IF AN 429C ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND. 430C 431C IERR IS SET TO 432C ZERO FOR NORMAL RETURN, 433C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED 434C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. 435C 436C CALLS CDIV FOR COMPLEX DIVISION. 437C 438C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 439C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 440C 441C THIS VERSION DATED AUGUST 1983. 442C 443C ------------------------------------------------------------------ 444* 445 EXTERNAL DLAMCH 446 DOUBLE PRECISION DLAMCH, UNFL,OVFL,ULP,SMLNUM,SMALL 447 IF (N.LE.0) RETURN 448* 449* INITIALIZE 450* 451 ITCNT = 0 452 OPST = 0 453C 454 IERR = 0 455 K = 1 456C .......... STORE ROOTS ISOLATED BY BALANC 457C AND COMPUTE MATRIX NORM .......... 458 DO 50 I = 1, N 459 IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50 460 WR(I) = H(I,I) 461 WI(I) = 0.0D0 462 50 CONTINUE 463 464* INCREMENT OPCOUNT FOR COMPUTING MATRIX NORM 465 OPS = OPS + (IGH-LOW+1)*(IGH-LOW+2)/2 466* 467* COMPUTE THE 1-NORM OF MATRIX H 468* 469 NORM = 0.0D0 470 DO 5 J = LOW, IGH 471 S = 0.0D0 472 DO 4 I = LOW, MIN(IGH,J+1) 473 S = S + DABS(H(I,J)) 474 4 CONTINUE 475 NORM = MAX(NORM, S) 476 5 CONTINUE 477C 478 UNFL = DLAMCH( 'SAFE MINIMUM' ) 479 OVFL = DLAMCH( 'OVERFLOW' ) 480 ULP = DLAMCH( 'EPSILON' )*DLAMCH( 'BASE' ) 481 SMLNUM = MAX( UNFL*( N / ULP ), N / ( ULP*OVFL ) ) 482 SMALL = MAX( SMLNUM, MIN( ( NORM*SMLNUM )*NORM, ULP*NORM ) ) 483C 484 EN = IGH 485 T = 0.0D0 486 ITN = 30*N 487C .......... SEARCH FOR NEXT EIGENVALUES .......... 488 60 IF (EN .LT. LOW) GO TO 340 489 ITS = 0 490 NA = EN - 1 491 ENM2 = NA - 1 492C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT 493C FOR L=EN STEP -1 UNTIL LOW DO -- .......... 494* REPLACE SPLITTING CRITERION WITH NEW ONE AS IN LAPACK 495* 496 70 DO 80 LL = LOW, EN 497 L = EN + LOW - LL 498 IF (L .EQ. LOW) GO TO 100 499 S = DABS(H(L-1,L-1)) + DABS(H(L,L)) 500 IF (S .EQ. 0.0D0) S = NORM 501 IF ( ABS(H(L,L-1)) .LE. MAX(ULP*S,SMALL) ) GO TO 100 502 80 CONTINUE 503C .......... FORM SHIFT .......... 504 100 CONTINUE 505* 506* INCREMENT OP COUNT FOR CONVERGENCE TEST 507 OPS = OPS + 2*(EN-L+1) 508 X = H(EN,EN) 509 IF (L .EQ. EN) GO TO 270 510 Y = H(NA,NA) 511 W = H(EN,NA) * H(NA,EN) 512 IF (L .EQ. NA) GO TO 280 513 IF (ITN .EQ. 0) GO TO 1000 514 IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 515C .......... FORM EXCEPTIONAL SHIFT .......... 516* 517* INCREMENT OP COUNT 518 OPS = OPS + (EN-LOW+6) 519 T = T + X 520C 521 DO 120 I = LOW, EN 522 120 H(I,I) = H(I,I) - X 523C 524 S = DABS(H(EN,NA)) + DABS(H(NA,ENM2)) 525 X = 0.75D0 * S 526 Y = X 527 W = -0.4375D0 * S * S 528 130 ITS = ITS + 1 529 ITN = ITN - 1 530* 531* UPDATE ITERATION NUMBER 532 ITCNT = 30*N - ITN 533C .......... LOOK FOR TWO CONSECUTIVE SMALL 534C SUB-DIAGONAL ELEMENTS. 535C FOR M=EN-2 STEP -1 UNTIL L DO -- .......... 536 DO 140 MM = L, ENM2 537 M = ENM2 + L - MM 538 ZZ = H(M,M) 539 R = X - ZZ 540 S = Y - ZZ 541 P = (R * S - W) / H(M+1,M) + H(M,M+1) 542 Q = H(M+1,M+1) - ZZ - R - S 543 R = H(M+2,M+1) 544 S = DABS(P) + DABS(Q) + DABS(R) 545 P = P / S 546 Q = Q / S 547 R = R / S 548 IF (M .EQ. L) GO TO 150 549 TST1 = DABS(P)*(DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1))) 550 TST2 = DABS(H(M,M-1))*(DABS(Q) + DABS(R)) 551 IF ( TST2 .LE. MAX(ULP*TST1,SMALL) ) GO TO 150 552 140 CONTINUE 553C 554 150 CONTINUE 555* 556* INCREMENT OPCOUNT FOR LOOP 140 557 OPST = OPST + 20*(ENM2-M+1) 558 MP2 = M + 2 559C 560 DO 160 I = MP2, EN 561 H(I,I-2) = 0.0D0 562 IF (I .EQ. MP2) GO TO 160 563 H(I,I-3) = 0.0D0 564 160 CONTINUE 565C .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND 566C COLUMNS M TO EN .......... 567* 568* INCREMENT OPCOUNT FOR LOOP 260 569 OPST = OPST + 18*(NA-M+1) 570 DO 260 K = M, NA 571 NOTLAS = K .NE. NA 572 IF (K .EQ. M) GO TO 170 573 P = H(K,K-1) 574 Q = H(K+1,K-1) 575 R = 0.0D0 576 IF (NOTLAS) R = H(K+2,K-1) 577 X = DABS(P) + DABS(Q) + DABS(R) 578 IF (X .EQ. 0.0D0) GO TO 260 579 P = P / X 580 Q = Q / X 581 R = R / X 582 170 S = DSIGN(DSQRT(P*P+Q*Q+R*R),P) 583 IF (K .EQ. M) GO TO 180 584 H(K,K-1) = -S * X 585 GO TO 190 586 180 IF (L .NE. M) H(K,K-1) = -H(K,K-1) 587 190 P = P + S 588 X = P / S 589 Y = Q / S 590 ZZ = R / S 591 Q = Q / P 592 R = R / P 593 IF (NOTLAS) GO TO 225 594C .......... ROW MODIFICATION .......... 595* 596* INCREMENT OP COUNT FOR LOOP 200 597 OPS = OPS + 6*(N-K+1) 598 DO 200 J = K, N 599 P = H(K,J) + Q * H(K+1,J) 600 H(K,J) = H(K,J) - P * X 601 H(K+1,J) = H(K+1,J) - P * Y 602 200 CONTINUE 603C 604 J = MIN0(EN,K+3) 605C .......... COLUMN MODIFICATION .......... 606* 607* INCREMENT OPCOUNT FOR LOOP 210 608 OPS = OPS + 6*J 609 DO 210 I = 1, J 610 P = X * H(I,K) + Y * H(I,K+1) 611 H(I,K) = H(I,K) - P 612 H(I,K+1) = H(I,K+1) - P * Q 613 210 CONTINUE 614C .......... ACCUMULATE TRANSFORMATIONS .......... 615* 616* INCREMENT OPCOUNT FOR LOOP 220 617 OPS = OPS + 6*(IGH-LOW + 1) 618 DO 220 I = LOW, IGH 619 P = X * Z(I,K) + Y * Z(I,K+1) 620 Z(I,K) = Z(I,K) - P 621 Z(I,K+1) = Z(I,K+1) - P * Q 622 220 CONTINUE 623 GO TO 255 624 225 CONTINUE 625C .......... ROW MODIFICATION .......... 626* 627* INCREMENT OPCOUNT FOR LOOP 230 628 OPS = OPS + 10*(N-K+1) 629 DO 230 J = K, N 630 P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J) 631 H(K,J) = H(K,J) - P * X 632 H(K+1,J) = H(K+1,J) - P * Y 633 H(K+2,J) = H(K+2,J) - P * ZZ 634 230 CONTINUE 635C 636 J = MIN0(EN,K+3) 637C .......... COLUMN MODIFICATION .......... 638* 639* INCREMENT OPCOUNT FOR LOOP 240 640 OPS = OPS + 10*J 641 DO 240 I = 1, J 642 P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2) 643 H(I,K) = H(I,K) - P 644 H(I,K+1) = H(I,K+1) - P * Q 645 H(I,K+2) = H(I,K+2) - P * R 646 240 CONTINUE 647C .......... ACCUMULATE TRANSFORMATIONS .......... 648* 649* INCREMENT OPCOUNT FOR LOOP 250 650 OPS = OPS + 10*(IGH-LOW+1) 651 DO 250 I = LOW, IGH 652 P = X * Z(I,K) + Y * Z(I,K+1) + ZZ * Z(I,K+2) 653 Z(I,K) = Z(I,K) - P 654 Z(I,K+1) = Z(I,K+1) - P * Q 655 Z(I,K+2) = Z(I,K+2) - P * R 656 250 CONTINUE 657 255 CONTINUE 658C 659 260 CONTINUE 660C 661 GO TO 70 662C .......... ONE ROOT FOUND .......... 663 270 H(EN,EN) = X + T 664 WR(EN) = H(EN,EN) 665 WI(EN) = 0.0D0 666 EN = NA 667 GO TO 60 668C .......... TWO ROOTS FOUND .......... 669 280 P = (Y - X) / 2.0D0 670 Q = P * P + W 671 ZZ = DSQRT(DABS(Q)) 672 H(EN,EN) = X + T 673 X = H(EN,EN) 674 H(NA,NA) = Y + T 675 IF (Q .LT. 0.0D0) GO TO 320 676C .......... REAL PAIR .......... 677 ZZ = P + DSIGN(ZZ,P) 678 WR(NA) = X + ZZ 679 WR(EN) = WR(NA) 680 IF (ZZ .NE. 0.0D0) WR(EN) = X - W / ZZ 681 WI(NA) = 0.0D0 682 WI(EN) = 0.0D0 683 X = H(EN,NA) 684 S = DABS(X) + DABS(ZZ) 685 P = X / S 686 Q = ZZ / S 687 R = DSQRT(P*P+Q*Q) 688 P = P / R 689 Q = Q / R 690* 691* INCREMENT OP COUNT FOR FINDING TWO ROOTS. 692 OPST = OPST + 18 693* 694* INCREMENT OP COUNT FOR MODIFICATION AND ACCUMULATION 695* IN LOOP 290, 300, 310 696 OPS = OPS + 6*(N-NA+1) + 6*EN + 6*(IGH-LOW+1) 697C .......... ROW MODIFICATION .......... 698 DO 290 J = NA, N 699 ZZ = H(NA,J) 700 H(NA,J) = Q * ZZ + P * H(EN,J) 701 H(EN,J) = Q * H(EN,J) - P * ZZ 702 290 CONTINUE 703C .......... COLUMN MODIFICATION .......... 704 DO 300 I = 1, EN 705 ZZ = H(I,NA) 706 H(I,NA) = Q * ZZ + P * H(I,EN) 707 H(I,EN) = Q * H(I,EN) - P * ZZ 708 300 CONTINUE 709C .......... ACCUMULATE TRANSFORMATIONS .......... 710 DO 310 I = LOW, IGH 711 ZZ = Z(I,NA) 712 Z(I,NA) = Q * ZZ + P * Z(I,EN) 713 Z(I,EN) = Q * Z(I,EN) - P * ZZ 714 310 CONTINUE 715C 716 GO TO 330 717C .......... COMPLEX PAIR .......... 718 320 WR(NA) = X + P 719 WR(EN) = X + P 720 WI(NA) = ZZ 721 WI(EN) = -ZZ 722* 723* INCREMENT OP COUNT FOR FINDING COMPLEX PAIR. 724 OPST = OPST + 9 725 330 EN = ENM2 726 GO TO 60 727C .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND 728C VECTORS OF UPPER TRIANGULAR FORM .......... 729 340 IF (NORM .EQ. 0.0D0) GO TO 1001 730C .......... FOR EN=N STEP -1 UNTIL 1 DO -- .......... 731 DO 800 NN = 1, N 732 EN = N + 1 - NN 733 P = WR(EN) 734 Q = WI(EN) 735 NA = EN - 1 736 IF (Q) 710, 600, 800 737C .......... REAL VECTOR .......... 738 600 M = EN 739 H(EN,EN) = 1.0D0 740 IF (NA .EQ. 0) GO TO 800 741C .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... 742 DO 700 II = 1, NA 743 I = EN - II 744 W = H(I,I) - P 745 R = 0.0D0 746C 747* 748* INCREMENT OP COUNT FOR LOOP 610 749 OPST = OPST + 2*(EN - M+1) 750 DO 610 J = M, EN 751 610 R = R + H(I,J) * H(J,EN) 752C 753 IF (WI(I) .GE. 0.0D0) GO TO 630 754 ZZ = W 755 S = R 756 GO TO 700 757 630 M = I 758 IF (WI(I) .NE. 0.0D0) GO TO 640 759 T = W 760 IF (T .NE. 0.0D0) GO TO 635 761 TST1 = NORM 762 T = TST1 763 632 T = 0.01D0 * T 764 TST2 = NORM + T 765 IF (TST2 .GT. TST1) GO TO 632 766 635 H(I,EN) = -R / T 767 GO TO 680 768C .......... SOLVE REAL EQUATIONS .......... 769 640 X = H(I,I+1) 770 Y = H(I+1,I) 771 Q = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) 772 T = (X * S - ZZ * R) / Q 773* 774* INCREMENT OP COUNT FOR SOLVING REAL EQUATION. 775 OPST = OPST + 13 776 H(I,EN) = T 777 IF (DABS(X) .LE. DABS(ZZ)) GO TO 650 778 H(I+1,EN) = (-R - W * T) / X 779 GO TO 680 780 650 H(I+1,EN) = (-S - Y * T) / ZZ 781C 782C .......... OVERFLOW CONTROL .......... 783 680 T = DABS(H(I,EN)) 784 IF (T .EQ. 0.0D0) GO TO 700 785 TST1 = T 786 TST2 = TST1 + 1.0D0/TST1 787 IF (TST2 .GT. TST1) GO TO 700 788* 789* INCREMENT OP COUNT. 790 OPST = OPST + (EN-I+1) 791 DO 690 J = I, EN 792 H(J,EN) = H(J,EN)/T 793 690 CONTINUE 794C 795 700 CONTINUE 796C .......... END REAL VECTOR .......... 797 GO TO 800 798C .......... COMPLEX VECTOR .......... 799 710 M = NA 800C .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT 801C EIGENVECTOR MATRIX IS TRIANGULAR .......... 802 IF (DABS(H(EN,NA)) .LE. DABS(H(NA,EN))) GO TO 720 803 H(NA,NA) = Q / H(EN,NA) 804 H(NA,EN) = -(H(EN,EN) - P) / H(EN,NA) 805* 806* INCREMENT OP COUNT. 807 OPST = OPST + 3 808 GO TO 730 809 720 CALL CDIV(0.0D0,-H(NA,EN),H(NA,NA)-P,Q,H(NA,NA),H(NA,EN)) 810* 811* INCREMENT OP COUNT IF (ABS(H(EN,NA)) .LE. ABS(H(NA,EN))) 812 OPST = OPST + 16 813 730 H(EN,NA) = 0.0D0 814 H(EN,EN) = 1.0D0 815 ENM2 = NA - 1 816 IF (ENM2 .EQ. 0) GO TO 800 817C .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- .......... 818 DO 795 II = 1, ENM2 819 I = NA - II 820 W = H(I,I) - P 821 RA = 0.0D0 822 SA = 0.0D0 823C 824* 825* INCREMENT OP COUNT FOR LOOP 760 826 OPST = OPST + 4*(EN-M+1) 827 DO 760 J = M, EN 828 RA = RA + H(I,J) * H(J,NA) 829 SA = SA + H(I,J) * H(J,EN) 830 760 CONTINUE 831C 832 IF (WI(I) .GE. 0.0D0) GO TO 770 833 ZZ = W 834 R = RA 835 S = SA 836 GO TO 795 837 770 M = I 838 IF (WI(I) .NE. 0.0D0) GO TO 780 839 CALL CDIV(-RA,-SA,W,Q,H(I,NA),H(I,EN)) 840* 841* INCREMENT OP COUNT FOR CDIV 842 OPST = OPST + 16 843 GO TO 790 844C .......... SOLVE COMPLEX EQUATIONS .......... 845 780 X = H(I,I+1) 846 Y = H(I+1,I) 847 VR = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) - Q * Q 848 VI = (WR(I) - P) * 2.0D0 * Q 849* 850* INCREMENT OPCOUNT (AVERAGE) FOR SOLVING COMPLEX EQUATIONS 851 OPST = OPST + 42 852 IF (VR .NE. 0.0D0 .OR. VI .NE. 0.0D0) GO TO 784 853 TST1 = NORM * (DABS(W) + DABS(Q) + DABS(X) 854 X + DABS(Y) + DABS(ZZ)) 855 VR = TST1 856 783 VR = 0.01D0 * VR 857 TST2 = TST1 + VR 858 IF (TST2 .GT. TST1) GO TO 783 859 784 CALL CDIV(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA,VR,VI, 860 X H(I,NA),H(I,EN)) 861 IF (DABS(X) .LE. DABS(ZZ) + DABS(Q)) GO TO 785 862 H(I+1,NA) = (-RA - W * H(I,NA) + Q * H(I,EN)) / X 863 H(I+1,EN) = (-SA - W * H(I,EN) - Q * H(I,NA)) / X 864 GO TO 790 865 785 CALL CDIV(-R-Y*H(I,NA),-S-Y*H(I,EN),ZZ,Q, 866 X H(I+1,NA),H(I+1,EN)) 867C 868C .......... OVERFLOW CONTROL .......... 869 790 T = DMAX1(DABS(H(I,NA)), DABS(H(I,EN))) 870 IF (T .EQ. 0.0D0) GO TO 795 871 TST1 = T 872 TST2 = TST1 + 1.0D0/TST1 873 IF (TST2 .GT. TST1) GO TO 795 874* 875* INCREMENT OP COUNT. 876 OPST = OPST + 2*(EN-I+1) 877 DO 792 J = I, EN 878 H(J,NA) = H(J,NA)/T 879 H(J,EN) = H(J,EN)/T 880 792 CONTINUE 881C 882 795 CONTINUE 883C .......... END COMPLEX VECTOR .......... 884 800 CONTINUE 885C .......... END BACK SUBSTITUTION. 886C VECTORS OF ISOLATED ROOTS .......... 887 DO 840 I = 1, N 888 IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840 889C 890 DO 820 J = I, N 891 820 Z(I,J) = H(I,J) 892C 893 840 CONTINUE 894C .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE 895C VECTORS OF ORIGINAL FULL MATRIX. 896C FOR J=N STEP -1 UNTIL LOW DO -- .......... 897 DO 880 JJ = LOW, N 898 J = N + LOW - JJ 899 M = MIN0(J,IGH) 900C 901* 902* INCREMENT OP COUNT. 903 OPS = OPS + 2*(IGH-LOW+1)*(M-LOW+1) 904 DO 880 I = LOW, IGH 905 ZZ = 0.0D0 906C 907 DO 860 K = LOW, M 908 860 ZZ = ZZ + Z(I,K) * H(K,J) 909C 910 Z(I,J) = ZZ 911 880 CONTINUE 912C 913 GO TO 1001 914C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT 915C CONVERGED AFTER 30*N ITERATIONS .......... 916 1000 IERR = EN 917 1001 CONTINUE 918* 919* COMPUTE FINAL OP COUNT 920 OPS = OPS + OPST 921 RETURN 922 END 923 SUBROUTINE IMTQL1(N,D,E,IERR) 924* 925* EISPACK ROUTINE 926* MODIFIED FOR COMPARISON WITH LAPACK ROUTINES. 927* 928* CONVERGENCE TEST WAS MODIFIED TO BE THE SAME AS IN DSTEQR. 929* 930C 931 INTEGER I,J,L,M,N,II,MML,IERR 932 DOUBLE PRECISION D(N),E(N) 933 DOUBLE PRECISION B,C,F,G,P,R,S,TST1,TST2,PYTHAG 934 DOUBLE PRECISION EPS, TST 935 DOUBLE PRECISION DLAMCH 936* 937* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 938* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 939* OPST IS USED TO ACCUMULATE CONTRIBUTIONS TO OPS FROM 940* FUNCTION PYTHAG. IT IS PASSED TO AND FROM PYTHAG 941* THROUGH COMMON BLOCK PYTHOP. 942* .. COMMON BLOCKS .. 943 COMMON / LATIME / OPS, ITCNT 944 COMMON / PYTHOP / OPST 945* 946* .. SCALARS IN COMMON .. 947 DOUBLE PRECISION ITCNT, OPS, OPST 948* .. 949C 950C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL1, 951C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, 952C AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. 953C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). 954C 955C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC 956C TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. 957C 958C ON INPUT 959C 960C N IS THE ORDER OF THE MATRIX. 961C 962C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. 963C 964C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX 965C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. 966C 967C ON OUTPUT 968C 969C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN 970C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND 971C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE 972C THE SMALLEST EIGENVALUES. 973C 974C E HAS BEEN DESTROYED. 975C 976C IERR IS SET TO 977C ZERO FOR NORMAL RETURN, 978C J IF THE J-TH EIGENVALUE HAS NOT BEEN 979C DETERMINED AFTER 40 ITERATIONS. 980C 981C CALLS PYTHAG FOR SQRT(A*A + B*B) . 982C 983C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 984C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 985C 986C THIS VERSION DATED AUGUST 1983. 987C 988C ------------------------------------------------------------------ 989C 990 IERR = 0 991 IF (N .EQ. 1) GO TO 1001 992* 993* INITIALIZE ITERATION COUNT AND OPST 994 ITCNT = 0 995 OPST = 0 996* 997* DETERMINE THE UNIT ROUNDOFF FOR THIS ENVIRONMENT. 998* 999 EPS = DLAMCH( 'EPSILON' ) 1000C 1001 DO 100 I = 2, N 1002 100 E(I-1) = E(I) 1003C 1004 E(N) = 0.0D0 1005C 1006 DO 290 L = 1, N 1007 J = 0 1008C .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... 1009 105 DO 110 M = L, N 1010 IF (M .EQ. N) GO TO 120 1011 TST = ABS( E(M) ) 1012 IF( TST .LE. EPS * ( ABS(D(M)) + ABS(D(M+1)) ) ) GO TO 120 1013* TST1 = ABS(D(M)) + ABS(D(M+1)) 1014* TST2 = TST1 + ABS(E(M)) 1015* IF (TST2 .EQ. TST1) GO TO 120 1016 110 CONTINUE 1017C 1018 120 P = D(L) 1019* 1020* INCREMENT OPCOUNT FOR FINDING SMALL SUBDIAGONAL ELEMENT. 1021 OPS = OPS + 2*( MIN(M,N-1)-L+1 ) 1022 IF (M .EQ. L) GO TO 215 1023 IF (J .EQ. 40) GO TO 1000 1024 J = J + 1 1025C .......... FORM SHIFT .......... 1026 G = (D(L+1) - P) / (2.0D0 * E(L)) 1027 R = PYTHAG(G,1.0D0) 1028 G = D(M) - P + E(L) / (G + DSIGN(R,G)) 1029* 1030* INCREMENT OPCOUNT FOR FORMING SHIFT. 1031 OPS = OPS + 7 1032 S = 1.0D0 1033 C = 1.0D0 1034 P = 0.0D0 1035 MML = M - L 1036C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... 1037 DO 200 II = 1, MML 1038 I = M - II 1039 F = S * E(I) 1040 B = C * E(I) 1041 R = PYTHAG(F,G) 1042 E(I+1) = R 1043 IF (R .EQ. 0.0D0) GO TO 210 1044 S = F / R 1045 C = G / R 1046 G = D(I+1) - P 1047 R = (D(I) - G) * S + 2.0D0 * C * B 1048 P = S * R 1049 D(I+1) = G + P 1050 G = C * R - B 1051 200 CONTINUE 1052C 1053 D(L) = D(L) - P 1054 E(L) = G 1055 E(M) = 0.0D0 1056* 1057* INCREMENT OPCOUNT FOR INNER LOOP. 1058 OPS = OPS + MML*14 + 1 1059* 1060* INCREMENT ITERATION COUNTER 1061 ITCNT = ITCNT + 1 1062 GO TO 105 1063C .......... RECOVER FROM UNDERFLOW .......... 1064 210 D(I+1) = D(I+1) - P 1065 E(M) = 0.0D0 1066* 1067* INCREMENT OPCOUNT FOR INNER LOOP, WHEN UNDERFLOW OCCURS. 1068 OPS = OPS + 2+(II-1)*14 + 1 1069 GO TO 105 1070C .......... ORDER EIGENVALUES .......... 1071 215 IF (L .EQ. 1) GO TO 250 1072C .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... 1073 DO 230 II = 2, L 1074 I = L + 2 - II 1075 IF (P .GE. D(I-1)) GO TO 270 1076 D(I) = D(I-1) 1077 230 CONTINUE 1078C 1079 250 I = 1 1080 270 D(I) = P 1081 290 CONTINUE 1082C 1083 GO TO 1001 1084C .......... SET ERROR -- NO CONVERGENCE TO AN 1085C EIGENVALUE AFTER 40 ITERATIONS .......... 1086 1000 IERR = L 1087 1001 CONTINUE 1088* 1089* COMPUTE FINAL OP COUNT 1090 OPS = OPS + OPST 1091 RETURN 1092 END 1093 SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR) 1094* 1095* EISPACK ROUTINE. MODIFIED FOR COMPARISON WITH LAPACK. 1096* 1097* CONVERGENCE TEST WAS MODIFIED TO BE THE SAME AS IN DSTEQR. 1098* 1099C 1100 INTEGER I,J,K,L,M,N,II,NM,MML,IERR 1101 DOUBLE PRECISION D(N),E(N),Z(NM,N) 1102 DOUBLE PRECISION B,C,F,G,P,R,S,TST1,TST2,PYTHAG 1103 DOUBLE PRECISION EPS, TST 1104 DOUBLE PRECISION DLAMCH 1105* 1106* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 1107* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 1108* OPST IS USED TO ACCUMULATE CONTRIBUTIONS TO OPS FROM 1109* FUNCTION PYTHAG. IT IS PASSED TO AND FROM PYTHAG 1110* THROUGH COMMON BLOCK PYTHOP. 1111* .. COMMON BLOCKS .. 1112 COMMON / LATIME / OPS, ITCNT 1113 COMMON / PYTHOP / OPST 1114* .. 1115* .. SCALARS IN COMMON .. 1116 DOUBLE PRECISION ITCNT, OPS, OPST 1117* .. 1118C 1119C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2, 1120C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, 1121C AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. 1122C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). 1123C 1124C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS 1125C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. 1126C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO 1127C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS 1128C FULL MATRIX TO TRIDIAGONAL FORM. 1129C 1130C ON INPUT 1131C 1132C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 1133C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 1134C DIMENSION STATEMENT. 1135C 1136C N IS THE ORDER OF THE MATRIX. 1137C 1138C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. 1139C 1140C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX 1141C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. 1142C 1143C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE 1144C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS 1145C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN 1146C THE IDENTITY MATRIX. 1147C 1148C ON OUTPUT 1149C 1150C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN 1151C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT 1152C UNORDERED FOR INDICES 1,2,...,IERR-1. 1153C 1154C E HAS BEEN DESTROYED. 1155C 1156C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC 1157C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, 1158C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED 1159C EIGENVALUES. 1160C 1161C IERR IS SET TO 1162C ZERO FOR NORMAL RETURN, 1163C J IF THE J-TH EIGENVALUE HAS NOT BEEN 1164C DETERMINED AFTER 40 ITERATIONS. 1165C 1166C CALLS PYTHAG FOR SQRT(A*A + B*B) . 1167C 1168C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 1169C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 1170C 1171C THIS VERSION DATED AUGUST 1983. 1172C 1173C ------------------------------------------------------------------ 1174C 1175 IERR = 0 1176 IF (N .EQ. 1) GO TO 1001 1177* 1178* INITIALIZE ITERATION COUNT AND OPST 1179 ITCNT = 0 1180 OPST = 0 1181* 1182* DETERMINE UNIT ROUNDOFF FOR THIS MACHINE. 1183 EPS = DLAMCH( 'EPSILON' ) 1184C 1185 DO 100 I = 2, N 1186 100 E(I-1) = E(I) 1187C 1188 E(N) = 0.0D0 1189C 1190 DO 240 L = 1, N 1191 J = 0 1192C .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... 1193 105 DO 110 M = L, N 1194 IF (M .EQ. N) GO TO 120 1195* TST1 = ABS(D(M)) + ABS(D(M+1)) 1196* TST2 = TST1 + ABS(E(M)) 1197* IF (TST2 .EQ. TST1) GO TO 120 1198 TST = ABS( E(M) ) 1199 IF( TST .LE. EPS * ( ABS(D(M)) + ABS(D(M+1)) ) ) GO TO 120 1200 110 CONTINUE 1201C 1202 120 P = D(L) 1203* 1204* INCREMENT OPCOUNT FOR FINDING SMALL SUBDIAGONAL ELEMENT. 1205 OPS = OPS + 2*( MIN(M,N)-L+1 ) 1206 IF (M .EQ. L) GO TO 240 1207 IF (J .EQ. 40) GO TO 1000 1208 J = J + 1 1209C .......... FORM SHIFT .......... 1210 G = (D(L+1) - P) / (2.0D0 * E(L)) 1211 R = PYTHAG(G,1.0D0) 1212 G = D(M) - P + E(L) / (G + DSIGN(R,G)) 1213* 1214* INCREMENT OPCOUNT FOR FORMING SHIFT. 1215 OPS = OPS + 7 1216 S = 1.0D0 1217 C = 1.0D0 1218 P = 0.0D0 1219 MML = M - L 1220C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... 1221 DO 200 II = 1, MML 1222 I = M - II 1223 F = S * E(I) 1224 B = C * E(I) 1225 R = PYTHAG(F,G) 1226 E(I+1) = R 1227 IF (R .EQ. 0.0D0) GO TO 210 1228 S = F / R 1229 C = G / R 1230 G = D(I+1) - P 1231 R = (D(I) - G) * S + 2.0D0 * C * B 1232 P = S * R 1233 D(I+1) = G + P 1234 G = C * R - B 1235C .......... FORM VECTOR .......... 1236 DO 180 K = 1, N 1237 F = Z(K,I+1) 1238 Z(K,I+1) = S * Z(K,I) + C * F 1239 Z(K,I) = C * Z(K,I) - S * F 1240 180 CONTINUE 1241C 1242 200 CONTINUE 1243C 1244 D(L) = D(L) - P 1245 E(L) = G 1246 E(M) = 0.0D0 1247* 1248* INCREMENT OPCOUNT FOR INNER LOOP. 1249 OPS = OPS + MML*( 14+6*N ) + 1 1250* 1251* INCREMENT ITERATION COUNTER 1252 ITCNT = ITCNT + 1 1253 GO TO 105 1254C .......... RECOVER FROM UNDERFLOW .......... 1255 210 D(I+1) = D(I+1) - P 1256 E(M) = 0.0D0 1257* 1258* INCREMENT OPCOUNT FOR INNER LOOP, WHEN UNDERFLOW OCCURS. 1259 OPS = OPS + 2+(II-1)*(14+6*N) + 1 1260 GO TO 105 1261 240 CONTINUE 1262C .......... ORDER EIGENVALUES AND EIGENVECTORS .......... 1263 DO 300 II = 2, N 1264 I = II - 1 1265 K = I 1266 P = D(I) 1267C 1268 DO 260 J = II, N 1269 IF (D(J) .GE. P) GO TO 260 1270 K = J 1271 P = D(J) 1272 260 CONTINUE 1273C 1274 IF (K .EQ. I) GO TO 300 1275 D(K) = D(I) 1276 D(I) = P 1277C 1278 DO 280 J = 1, N 1279 P = Z(J,I) 1280 Z(J,I) = Z(J,K) 1281 Z(J,K) = P 1282 280 CONTINUE 1283C 1284 300 CONTINUE 1285C 1286 GO TO 1001 1287C .......... SET ERROR -- NO CONVERGENCE TO AN 1288C EIGENVALUE AFTER 40 ITERATIONS .......... 1289 1000 IERR = L 1290 1001 CONTINUE 1291* 1292* COMPUTE FINAL OP COUNT 1293 OPS = OPS + OPST 1294 RETURN 1295 END 1296 SUBROUTINE INVIT(NM,N,A,WR,WI,SELECT,MM,M,Z,IERR,RM1,RV1,RV2) 1297C 1298 INTEGER I,J,K,L,M,N,S,II,IP,MM,MP,NM,NS,N1,UK,IP1,ITS,KM1,IERR 1299 DOUBLE PRECISION A(NM,N),WR(N),WI(N),Z(NM,MM),RM1(N,N), 1300 X RV1(N),RV2(N) 1301 DOUBLE PRECISION T,W,X,Y,EPS3,NORM,NORMV,GROWTO,ILAMBD, 1302 X PYTHAG,RLAMBD,UKROOT 1303 LOGICAL SELECT(N) 1304* 1305* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 1306* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 1307* OPST IS USED TO ACCUMULATE SMALL CONTRIBUTIONS TO OPS 1308* TO AVOID ROUNDOFF ERROR 1309* .. COMMON BLOCKS .. 1310 COMMON /LATIME/ OPS, ITCNT 1311* .. 1312* .. SCALARS IN COMMON .. 1313 DOUBLE PRECISION OPS, ITCNT, OPST 1314* .. 1315C 1316C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE INVIT 1317C BY PETERS AND WILKINSON. 1318C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). 1319C 1320C THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL UPPER 1321C HESSENBERG MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, 1322C USING INVERSE ITERATION. 1323C 1324C ON INPUT 1325C 1326C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 1327C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 1328C DIMENSION STATEMENT. 1329C 1330C N IS THE ORDER OF THE MATRIX. 1331C 1332C A CONTAINS THE HESSENBERG MATRIX. 1333C 1334C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, 1335C OF THE EIGENVALUES OF THE MATRIX. THE EIGENVALUES MUST BE 1336C STORED IN A MANNER IDENTICAL TO THAT OF SUBROUTINE HQR, 1337C WHICH RECOGNIZES POSSIBLE SPLITTING OF THE MATRIX. 1338C 1339C SELECT SPECIFIES THE EIGENVECTORS TO BE FOUND. THE 1340C EIGENVECTOR CORRESPONDING TO THE J-TH EIGENVALUE IS 1341C SPECIFIED BY SETTING SELECT(J) TO .TRUE.. 1342C 1343C MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF 1344C COLUMNS REQUIRED TO STORE THE EIGENVECTORS TO BE FOUND. 1345C NOTE THAT TWO COLUMNS ARE REQUIRED TO STORE THE 1346C EIGENVECTOR CORRESPONDING TO A COMPLEX EIGENVALUE. 1347C 1348C ON OUTPUT 1349C 1350C A AND WI ARE UNALTERED. 1351C 1352C WR MAY HAVE BEEN ALTERED SINCE CLOSE EIGENVALUES ARE PERTURBED 1353C SLIGHTLY IN SEARCHING FOR INDEPENDENT EIGENVECTORS. 1354C 1355C SELECT MAY HAVE BEEN ALTERED. IF THE ELEMENTS CORRESPONDING 1356C TO A PAIR OF CONJUGATE COMPLEX EIGENVALUES WERE EACH 1357C INITIALLY SET TO .TRUE., THE PROGRAM RESETS THE SECOND OF 1358C THE TWO ELEMENTS TO .FALSE.. 1359C 1360C M IS THE NUMBER OF COLUMNS ACTUALLY USED TO STORE 1361C THE EIGENVECTORS. 1362C 1363C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. 1364C IF THE NEXT SELECTED EIGENVALUE IS REAL, THE NEXT COLUMN 1365C OF Z CONTAINS ITS EIGENVECTOR. IF THE EIGENVALUE IS 1366C COMPLEX, THE NEXT TWO COLUMNS OF Z CONTAIN THE REAL AND 1367C IMAGINARY PARTS OF ITS EIGENVECTOR. THE EIGENVECTORS ARE 1368C NORMALIZED SO THAT THE COMPONENT OF LARGEST MAGNITUDE IS 1. 1369C ANY VECTOR WHICH FAILS THE ACCEPTANCE TEST IS SET TO ZERO. 1370C 1371C IERR IS SET TO 1372C ZERO FOR NORMAL RETURN, 1373C -(2*N+1) IF MORE THAN MM COLUMNS OF Z ARE NECESSARY 1374C TO STORE THE EIGENVECTORS CORRESPONDING TO 1375C THE SPECIFIED EIGENVALUES. 1376C -K IF THE ITERATION CORRESPONDING TO THE K-TH 1377C VALUE FAILS, 1378C -(N+K) IF BOTH ERROR SITUATIONS OCCUR. 1379C 1380C RM1, RV1, AND RV2 ARE TEMPORARY STORAGE ARRAYS. NOTE THAT RM1 1381C IS SQUARE OF DIMENSION N BY N AND, AUGMENTED BY TWO COLUMNS 1382C OF Z, IS THE TRANSPOSE OF THE CORRESPONDING ALGOL B ARRAY. 1383C 1384C THE ALGOL PROCEDURE GUESSVEC APPEARS IN INVIT IN LINE. 1385C 1386C CALLS CDIV FOR COMPLEX DIVISION. 1387C CALLS PYTHAG FOR SQRT(A*A + B*B) . 1388C 1389C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 1390C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 1391C 1392C THIS VERSION DATED AUGUST 1983. 1393C 1394C ------------------------------------------------------------------ 1395* 1396* GET ULP FROM DLAMCH FOR NEW SMALL PERTURBATION AS IN LAPACK 1397 EXTERNAL DLAMCH 1398 DOUBLE PRECISION DLAMCH, ULP 1399 IF (N.LE.0) RETURN 1400 ULP = DLAMCH( 'EPSILON' ) 1401C 1402* 1403* INITIALIZE 1404 OPST = 0 1405 IERR = 0 1406 UK = 0 1407 S = 1 1408C .......... IP = 0, REAL EIGENVALUE 1409C 1, FIRST OF CONJUGATE COMPLEX PAIR 1410C -1, SECOND OF CONJUGATE COMPLEX PAIR .......... 1411 IP = 0 1412 N1 = N - 1 1413C 1414 DO 980 K = 1, N 1415 IF (WI(K) .EQ. 0.0D0 .OR. IP .LT. 0) GO TO 100 1416 IP = 1 1417 IF (SELECT(K) .AND. SELECT(K+1)) SELECT(K+1) = .FALSE. 1418 100 IF (.NOT. SELECT(K)) GO TO 960 1419 IF (WI(K) .NE. 0.0D0) S = S + 1 1420 IF (S .GT. MM) GO TO 1000 1421 IF (UK .GE. K) GO TO 200 1422C .......... CHECK FOR POSSIBLE SPLITTING .......... 1423 DO 120 UK = K, N 1424 IF (UK .EQ. N) GO TO 140 1425 IF (A(UK+1,UK) .EQ. 0.0D0) GO TO 140 1426 120 CONTINUE 1427C .......... COMPUTE INFINITY NORM OF LEADING UK BY UK 1428C (HESSENBERG) MATRIX .......... 1429 140 NORM = 0.0D0 1430 MP = 1 1431C 1432* 1433* INCREMENT OPCOUNT FOR COMPUTING MATRIX NORM 1434 OPS = OPS + UK*(UK-1)/2 1435 DO 180 I = 1, UK 1436 X = 0.0D0 1437C 1438 DO 160 J = MP, UK 1439 160 X = X + DABS(A(I,J)) 1440C 1441 IF (X .GT. NORM) NORM = X 1442 MP = I 1443 180 CONTINUE 1444C .......... EPS3 REPLACES ZERO PIVOT IN DECOMPOSITION 1445C AND CLOSE ROOTS ARE MODIFIED BY EPS3 .......... 1446 IF (NORM .EQ. 0.0D0) NORM = 1.0D0 1447* EPS3 = EPSLON(NORM) 1448* 1449* INCREMENT OPCOUNT 1450 OPST = OPST + 3 1451 EPS3 = NORM*ULP 1452C .......... GROWTO IS THE CRITERION FOR THE GROWTH .......... 1453 UKROOT = UK 1454 UKROOT = DSQRT(UKROOT) 1455 GROWTO = 0.1D0 / UKROOT 1456 200 RLAMBD = WR(K) 1457 ILAMBD = WI(K) 1458 IF (K .EQ. 1) GO TO 280 1459 KM1 = K - 1 1460 GO TO 240 1461C .......... PERTURB EIGENVALUE IF IT IS CLOSE 1462C TO ANY PREVIOUS EIGENVALUE .......... 1463 220 RLAMBD = RLAMBD + EPS3 1464C .......... FOR I=K-1 STEP -1 UNTIL 1 DO -- .......... 1465 240 DO 260 II = 1, KM1 1466 I = K - II 1467 IF (SELECT(I) .AND. DABS(WR(I)-RLAMBD) .LT. EPS3 .AND. 1468 X DABS(WI(I)-ILAMBD) .LT. EPS3) GO TO 220 1469 260 CONTINUE 1470* 1471* INCREMENT OPCOUNT FOR LOOP 260 (ASSUME THAT ALL EIGENVALUES 1472* ARE DIFFERENT) 1473 OPST = OPST + 2*(K-1) 1474C 1475 WR(K) = RLAMBD 1476C .......... PERTURB CONJUGATE EIGENVALUE TO MATCH .......... 1477 IP1 = K + IP 1478 WR(IP1) = RLAMBD 1479C .......... FORM UPPER HESSENBERG A-RLAMBD*I (TRANSPOSED) 1480C AND INITIAL REAL VECTOR .......... 1481 280 MP = 1 1482C 1483* 1484* INCREMENT OP COUNT FOR LOOP 320 1485 OPS = OPS + UK 1486 DO 320 I = 1, UK 1487C 1488 DO 300 J = MP, UK 1489 300 RM1(J,I) = A(I,J) 1490C 1491 RM1(I,I) = RM1(I,I) - RLAMBD 1492 MP = I 1493 RV1(I) = EPS3 1494 320 CONTINUE 1495C 1496 ITS = 0 1497 IF (ILAMBD .NE. 0.0D0) GO TO 520 1498C .......... REAL EIGENVALUE. 1499C TRIANGULAR DECOMPOSITION WITH INTERCHANGES, 1500C REPLACING ZERO PIVOTS BY EPS3 .......... 1501 IF (UK .EQ. 1) GO TO 420 1502C 1503* 1504* INCREMENT OPCOUNT LU DECOMPOSITION 1505 OPS = OPS + (UK-1)*(UK+2) 1506 DO 400 I = 2, UK 1507 MP = I - 1 1508 IF (DABS(RM1(MP,I)) .LE. DABS(RM1(MP,MP))) GO TO 360 1509C 1510 DO 340 J = MP, UK 1511 Y = RM1(J,I) 1512 RM1(J,I) = RM1(J,MP) 1513 RM1(J,MP) = Y 1514 340 CONTINUE 1515C 1516 360 IF (RM1(MP,MP) .EQ. 0.0D0) RM1(MP,MP) = EPS3 1517 X = RM1(MP,I) / RM1(MP,MP) 1518 IF (X .EQ. 0.0D0) GO TO 400 1519C 1520 DO 380 J = I, UK 1521 380 RM1(J,I) = RM1(J,I) - X * RM1(J,MP) 1522C 1523 400 CONTINUE 1524C 1525 420 IF (RM1(UK,UK) .EQ. 0.0D0) RM1(UK,UK) = EPS3 1526C .......... BACK SUBSTITUTION FOR REAL VECTOR 1527C FOR I=UK STEP -1 UNTIL 1 DO -- .......... 1528 440 DO 500 II = 1, UK 1529 I = UK + 1 - II 1530 Y = RV1(I) 1531 IF (I .EQ. UK) GO TO 480 1532 IP1 = I + 1 1533C 1534 DO 460 J = IP1, UK 1535 460 Y = Y - RM1(J,I) * RV1(J) 1536C 1537 480 RV1(I) = Y / RM1(I,I) 1538 500 CONTINUE 1539* 1540* INCREMENT OP COUNT FOR BACK SUBSTITUTION LOOP 500 1541 OPS = OPS + UK*(UK+1) 1542C 1543 GO TO 740 1544C .......... COMPLEX EIGENVALUE. 1545C TRIANGULAR DECOMPOSITION WITH INTERCHANGES, 1546C REPLACING ZERO PIVOTS BY EPS3. STORE IMAGINARY 1547C PARTS IN UPPER TRIANGLE STARTING AT (1,3) .......... 1548 520 NS = N - S 1549 Z(1,S-1) = -ILAMBD 1550 Z(1,S) = 0.0D0 1551 IF (N .EQ. 2) GO TO 550 1552 RM1(1,3) = -ILAMBD 1553 Z(1,S-1) = 0.0D0 1554 IF (N .EQ. 3) GO TO 550 1555C 1556 DO 540 I = 4, N 1557 540 RM1(1,I) = 0.0D0 1558C 1559 550 DO 640 I = 2, UK 1560 MP = I - 1 1561 W = RM1(MP,I) 1562 IF (I .LT. N) T = RM1(MP,I+1) 1563 IF (I .EQ. N) T = Z(MP,S-1) 1564 X = RM1(MP,MP) * RM1(MP,MP) + T * T 1565 IF (W * W .LE. X) GO TO 580 1566 X = RM1(MP,MP) / W 1567 Y = T / W 1568 RM1(MP,MP) = W 1569 IF (I .LT. N) RM1(MP,I+1) = 0.0D0 1570 IF (I .EQ. N) Z(MP,S-1) = 0.0D0 1571C 1572* 1573* INCREMENT OPCOUNT FOR LOOP 560 1574 OPS = OPS + 4*(UK-I+1) 1575 DO 560 J = I, UK 1576 W = RM1(J,I) 1577 RM1(J,I) = RM1(J,MP) - X * W 1578 RM1(J,MP) = W 1579 IF (J .LT. N1) GO TO 555 1580 L = J - NS 1581 Z(I,L) = Z(MP,L) - Y * W 1582 Z(MP,L) = 0.0D0 1583 GO TO 560 1584 555 RM1(I,J+2) = RM1(MP,J+2) - Y * W 1585 RM1(MP,J+2) = 0.0D0 1586 560 CONTINUE 1587C 1588 RM1(I,I) = RM1(I,I) - Y * ILAMBD 1589 IF (I .LT. N1) GO TO 570 1590 L = I - NS 1591 Z(MP,L) = -ILAMBD 1592 Z(I,L) = Z(I,L) + X * ILAMBD 1593 GO TO 640 1594 570 RM1(MP,I+2) = -ILAMBD 1595 RM1(I,I+2) = RM1(I,I+2) + X * ILAMBD 1596 GO TO 640 1597 580 IF (X .NE. 0.0D0) GO TO 600 1598 RM1(MP,MP) = EPS3 1599 IF (I .LT. N) RM1(MP,I+1) = 0.0D0 1600 IF (I .EQ. N) Z(MP,S-1) = 0.0D0 1601 T = 0.0D0 1602 X = EPS3 * EPS3 1603 600 W = W / X 1604 X = RM1(MP,MP) * W 1605 Y = -T * W 1606C 1607* 1608* INCREMENT OPCOUNT FOR LOOP 620 1609 OPS = OPS + 6*(UK-I+1) 1610 DO 620 J = I, UK 1611 IF (J .LT. N1) GO TO 610 1612 L = J - NS 1613 T = Z(MP,L) 1614 Z(I,L) = -X * T - Y * RM1(J,MP) 1615 GO TO 615 1616 610 T = RM1(MP,J+2) 1617 RM1(I,J+2) = -X * T - Y * RM1(J,MP) 1618 615 RM1(J,I) = RM1(J,I) - X * RM1(J,MP) + Y * T 1619 620 CONTINUE 1620C 1621 IF (I .LT. N1) GO TO 630 1622 L = I - NS 1623 Z(I,L) = Z(I,L) - ILAMBD 1624 GO TO 640 1625 630 RM1(I,I+2) = RM1(I,I+2) - ILAMBD 1626 640 CONTINUE 1627* 1628* INCREMENT OP COUNT (AVERAGE) FOR COMPUTING 1629* THE SCALARS IN LOOP 640 1630 OPS = OPS + 10*(UK -1) 1631C 1632 IF (UK .LT. N1) GO TO 650 1633 L = UK - NS 1634 T = Z(UK,L) 1635 GO TO 655 1636 650 T = RM1(UK,UK+2) 1637 655 IF (RM1(UK,UK) .EQ. 0.0D0 .AND. T .EQ. 0.0D0) RM1(UK,UK) = EPS3 1638C .......... BACK SUBSTITUTION FOR COMPLEX VECTOR 1639C FOR I=UK STEP -1 UNTIL 1 DO -- .......... 1640 660 DO 720 II = 1, UK 1641 I = UK + 1 - II 1642 X = RV1(I) 1643 Y = 0.0D0 1644 IF (I .EQ. UK) GO TO 700 1645 IP1 = I + 1 1646C 1647 DO 680 J = IP1, UK 1648 IF (J .LT. N1) GO TO 670 1649 L = J - NS 1650 T = Z(I,L) 1651 GO TO 675 1652 670 T = RM1(I,J+2) 1653 675 X = X - RM1(J,I) * RV1(J) + T * RV2(J) 1654 Y = Y - RM1(J,I) * RV2(J) - T * RV1(J) 1655 680 CONTINUE 1656C 1657 700 IF (I .LT. N1) GO TO 710 1658 L = I - NS 1659 T = Z(I,L) 1660 GO TO 715 1661 710 T = RM1(I,I+2) 1662 715 CALL CDIV(X,Y,RM1(I,I),T,RV1(I),RV2(I)) 1663 720 CONTINUE 1664* 1665* INCREMENT OP COUNT FOR LOOP 720. 1666 OPS = OPS + 4*UK*(UK+3) 1667C .......... ACCEPTANCE TEST FOR REAL OR COMPLEX 1668C EIGENVECTOR AND NORMALIZATION .......... 1669 740 ITS = ITS + 1 1670 NORM = 0.0D0 1671 NORMV = 0.0D0 1672C 1673 DO 780 I = 1, UK 1674 IF (ILAMBD .EQ. 0.0D0) X = DABS(RV1(I)) 1675 IF (ILAMBD .NE. 0.0D0) X = PYTHAG(RV1(I),RV2(I)) 1676 IF (NORMV .GE. X) GO TO 760 1677 NORMV = X 1678 J = I 1679 760 NORM = NORM + X 1680 780 CONTINUE 1681* 1682* INCREMENT OP COUNT ACCEPTANCE TEST 1683 IF (ILAMBD .EQ. 0.0D0) OPS = OPS + UK 1684 IF (ILAMBD .NE. 0.0D0) OPS = OPS + 16*UK 1685C 1686 IF (NORM .LT. GROWTO) GO TO 840 1687C .......... ACCEPT VECTOR .......... 1688 X = RV1(J) 1689 IF (ILAMBD .EQ. 0.0D0) X = 1.0D0 / X 1690 IF (ILAMBD .NE. 0.0D0) Y = RV2(J) 1691C 1692* 1693* INCREMENT OPCOUNT FOR LOOP 820 1694 IF (ILAMBD .EQ. 0.0D0) OPS = OPS + UK 1695 IF (ILAMBD .NE. 0.0D0) OPS = OPS + 16*UK 1696 DO 820 I = 1, UK 1697 IF (ILAMBD .NE. 0.0D0) GO TO 800 1698 Z(I,S) = RV1(I) * X 1699 GO TO 820 1700 800 CALL CDIV(RV1(I),RV2(I),X,Y,Z(I,S-1),Z(I,S)) 1701 820 CONTINUE 1702C 1703 IF (UK .EQ. N) GO TO 940 1704 J = UK + 1 1705 GO TO 900 1706C .......... IN-LINE PROCEDURE FOR CHOOSING 1707C A NEW STARTING VECTOR .......... 1708 840 IF (ITS .GE. UK) GO TO 880 1709 X = UKROOT 1710 Y = EPS3 / (X + 1.0D0) 1711 RV1(1) = EPS3 1712C 1713 DO 860 I = 2, UK 1714 860 RV1(I) = Y 1715C 1716 J = UK - ITS + 1 1717 RV1(J) = RV1(J) - EPS3 * X 1718 IF (ILAMBD .EQ. 0.0D0) GO TO 440 1719 GO TO 660 1720C .......... SET ERROR -- UNACCEPTED EIGENVECTOR .......... 1721 880 J = 1 1722 IERR = -K 1723C .......... SET REMAINING VECTOR COMPONENTS TO ZERO .......... 1724 900 DO 920 I = J, N 1725 Z(I,S) = 0.0D0 1726 IF (ILAMBD .NE. 0.0D0) Z(I,S-1) = 0.0D0 1727 920 CONTINUE 1728C 1729 940 S = S + 1 1730 960 IF (IP .EQ. (-1)) IP = 0 1731 IF (IP .EQ. 1) IP = -1 1732 980 CONTINUE 1733C 1734 GO TO 1001 1735C .......... SET ERROR -- UNDERESTIMATE OF EIGENVECTOR 1736C SPACE REQUIRED .......... 1737 1000 IF (IERR .NE. 0) IERR = IERR - N 1738 IF (IERR .EQ. 0) IERR = -(2 * N + 1) 1739 1001 M = S - 1 - IABS(IP) 1740* 1741* COMPUTE FINAL OP COUNT 1742 OPS = OPS + OPST 1743 RETURN 1744 END 1745 SUBROUTINE ORTHES(NM,N,LOW,IGH,A,ORT) 1746C 1747 INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW 1748 DOUBLE PRECISION A(NM,N),ORT(IGH) 1749 DOUBLE PRECISION F,G,H,SCALE 1750* 1751* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 1752* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 1753* OPST IS USED TO ACCUMULATE SMALL CONTRIBUTIONS TO OPS 1754* TO AVOID ROUNDOFF ERROR 1755* .. COMMON BLOCKS .. 1756 COMMON /LATIME/ OPS, ITCNT 1757* .. 1758* .. SCALARS IN COMMON .. 1759 DOUBLE PRECISION OPS, ITCNT, OPST 1760* .. 1761C 1762C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTHES, 1763C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. 1764C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). 1765C 1766C GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE 1767C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS 1768C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY 1769C ORTHOGONAL SIMILARITY TRANSFORMATIONS. 1770C 1771C ON INPUT 1772C 1773C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 1774C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 1775C DIMENSION STATEMENT. 1776C 1777C N IS THE ORDER OF THE MATRIX. 1778C 1779C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING 1780C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, 1781C SET LOW=1, IGH=N. 1782C 1783C A CONTAINS THE INPUT MATRIX. 1784C 1785C ON OUTPUT 1786C 1787C A CONTAINS THE HESSENBERG MATRIX. INFORMATION ABOUT 1788C THE ORTHOGONAL TRANSFORMATIONS USED IN THE REDUCTION 1789C IS STORED IN THE REMAINING TRIANGLE UNDER THE 1790C HESSENBERG MATRIX. 1791C 1792C ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. 1793C ONLY ELEMENTS LOW THROUGH IGH ARE USED. 1794C 1795C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 1796C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 1797C 1798C THIS VERSION DATED AUGUST 1983. 1799C 1800C ------------------------------------------------------------------ 1801C 1802 IF (N.LE.0) RETURN 1803 LA = IGH - 1 1804 KP1 = LOW + 1 1805 IF (LA .LT. KP1) GO TO 200 1806C 1807* 1808* INCREMENT OP COUNR FOR COMPUTING G,H,ORT(M),.. IN LOOP 180 1809 OPS = OPS + 6*(LA - KP1 + 1) 1810 DO 180 M = KP1, LA 1811 H = 0.0D0 1812 ORT(M) = 0.0D0 1813 SCALE = 0.0D0 1814C .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) .......... 1815* 1816* INCREMENT OP COUNT FOR LOOP 90 1817 OPS = OPS + (IGH-M +1) 1818 DO 90 I = M, IGH 1819 90 SCALE = SCALE + DABS(A(I,M-1)) 1820C 1821 IF (SCALE .EQ. 0.0D0) GO TO 180 1822 MP = M + IGH 1823C .......... FOR I=IGH STEP -1 UNTIL M DO -- .......... 1824* 1825* INCREMENT OP COUNT FOR LOOP 100 1826 OPS = OPS + 3*(IGH-M+1) 1827 DO 100 II = M, IGH 1828 I = MP - II 1829 ORT(I) = A(I,M-1) / SCALE 1830 H = H + ORT(I) * ORT(I) 1831 100 CONTINUE 1832C 1833 G = -DSIGN(DSQRT(H),ORT(M)) 1834 H = H - ORT(M) * G 1835 ORT(M) = ORT(M) - G 1836C .......... FORM (I-(U*UT)/H) * A .......... 1837* 1838* INCREMENT OP COUNT FOR LOOP 130 AND 160 1839 OPS = OPS + (N-M+1+IGH)*(4*(IGH-M+1) + 1) 1840 DO 130 J = M, N 1841 F = 0.0D0 1842C .......... FOR I=IGH STEP -1 UNTIL M DO -- .......... 1843 DO 110 II = M, IGH 1844 I = MP - II 1845 F = F + ORT(I) * A(I,J) 1846 110 CONTINUE 1847C 1848 F = F / H 1849C 1850 DO 120 I = M, IGH 1851 120 A(I,J) = A(I,J) - F * ORT(I) 1852C 1853 130 CONTINUE 1854C .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) .......... 1855 DO 160 I = 1, IGH 1856 F = 0.0D0 1857C .......... FOR J=IGH STEP -1 UNTIL M DO -- .......... 1858 DO 140 JJ = M, IGH 1859 J = MP - JJ 1860 F = F + ORT(J) * A(I,J) 1861 140 CONTINUE 1862C 1863 F = F / H 1864C 1865 DO 150 J = M, IGH 1866 150 A(I,J) = A(I,J) - F * ORT(J) 1867C 1868 160 CONTINUE 1869C 1870 ORT(M) = SCALE * ORT(M) 1871 A(M,M-1) = SCALE * G 1872 180 CONTINUE 1873C 1874 200 RETURN 1875 END 1876 DOUBLE PRECISION FUNCTION PYTHAG(A,B) 1877 DOUBLE PRECISION A,B 1878C 1879C FINDS SQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW 1880C 1881* 1882* COMMON BLOCK TO RETURN OPERATION COUNT 1883* OPST IS ONLY INCREMENTED HERE 1884* .. COMMON BLOCKS .. 1885 COMMON / PYTHOP / OPST 1886* .. 1887* .. SCALARS IN COMMON 1888 DOUBLE PRECISION OPST 1889* .. 1890 DOUBLE PRECISION P,R,S,T,U 1891 P = DMAX1(DABS(A),DABS(B)) 1892 IF (P .EQ. 0.0D0) GO TO 20 1893 R = (DMIN1(DABS(A),DABS(B))/P)**2 1894* 1895* INCREMENT OPST 1896 OPST = OPST + 2 1897 10 CONTINUE 1898 T = 4.0D0 + R 1899 IF (T .EQ. 4.0D0) GO TO 20 1900 S = R/T 1901 U = 1.0D0 + 2.0D0*S 1902 P = U*P 1903 R = (S/U)**2 * R 1904* 1905* INCREMENT OPST 1906 OPST = OPST + 8 1907 GO TO 10 1908 20 PYTHAG = P 1909 RETURN 1910 END 1911 SUBROUTINE TQLRAT(N,D,E2,IERR) 1912* 1913* EISPACK ROUTINE. 1914* MODIFIED FOR COMPARISON WITH LAPACK ROUTINES. 1915* 1916* CONVERGENCE TEST WAS MODIFIED TO BE THE SAME AS IN DSTEQR. 1917* 1918C 1919 INTEGER I,J,L,M,N,II,L1,MML,IERR 1920 DOUBLE PRECISION D(N),E2(N) 1921 DOUBLE PRECISION B,C,F,G,H,P,R,S,T,EPSLON,PYTHAG 1922 DOUBLE PRECISION EPS, TST 1923 DOUBLE PRECISION DLAMCH 1924* 1925* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 1926* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 1927* OPST IS USED TO ACCUMULATE CONTRIBUTIONS TO OPS FROM 1928* FUNCTION PYTHAG. IT IS PASSED TO AND FROM PYTHAG 1929* THROUGH COMMON BLOCK PYTHOP. 1930* .. COMMON BLOCKS .. 1931 COMMON / LATIME / OPS, ITCNT 1932 COMMON / PYTHOP / OPST 1933* .. 1934* .. SCALARS IN COMMON .. 1935 DOUBLE PRECISION ITCNT, OPS, OPST 1936* .. 1937C 1938C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT, 1939C ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH. 1940C 1941C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC 1942C TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD. 1943C 1944C ON INPUT 1945C 1946C N IS THE ORDER OF THE MATRIX. 1947C 1948C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. 1949C 1950C E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE 1951C INPUT MATRIX IN ITS LAST N-1 POSITIONS. E2(1) IS ARBITRARY. 1952C 1953C ON OUTPUT 1954C 1955C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN 1956C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND 1957C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE 1958C THE SMALLEST EIGENVALUES. 1959C 1960C E2 HAS BEEN DESTROYED. 1961C 1962C IERR IS SET TO 1963C ZERO FOR NORMAL RETURN, 1964C J IF THE J-TH EIGENVALUE HAS NOT BEEN 1965C DETERMINED AFTER 30 ITERATIONS. 1966C 1967C CALLS PYTHAG FOR SQRT(A*A + B*B) . 1968C 1969C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 1970C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 1971C 1972C THIS VERSION DATED AUGUST 1983. 1973C 1974C ------------------------------------------------------------------ 1975C 1976 IERR = 0 1977 IF (N .EQ. 1) GO TO 1001 1978* 1979* INITIALIZE ITERATION COUNT AND OPST 1980 ITCNT = 0 1981 OPST = 0 1982* 1983* DETERMINE THE UNIT ROUNDOFF FOR THIS ENVIRONMENT. 1984* 1985 EPS = DLAMCH( 'EPSILON' ) 1986C 1987 DO 100 I = 2, N 1988 100 E2(I-1) = E2(I) 1989C 1990 F = 0.0D0 1991 T = 0.0D0 1992 E2(N) = 0.0D0 1993C 1994 DO 290 L = 1, N 1995 J = 0 1996 H = DABS(D(L)) + DSQRT(E2(L)) 1997 IF (T .GT. H) GO TO 105 1998 T = H 1999 B = EPSLON(T) 2000 C = B * B 2001* 2002* INCREMENT OPCOUNT FOR THIS SECTION. 2003* (FUNCTION EPSLON IS COUNTED AS 6 FLOPS. THIS IS THE MINIMUM 2004* NUMBER REQUIRED, BUT COUNTING THEM EXACTLY WOULD AFFECT 2005* THE TIMING.) 2006 OPS = OPS + 9 2007C .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT .......... 2008 105 DO 110 M = L, N 2009 IF( M .EQ. N ) GO TO 120 2010 TST = SQRT( ABS( E2(M) ) ) 2011 IF( TST .LE. EPS * ( ABS(D(M)) + ABS(D(M+1)) ) ) GO TO 120 2012* IF (E2(M) .LE. C) GO TO 120 2013C .......... E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT 2014C THROUGH THE BOTTOM OF THE LOOP .......... 2015 110 CONTINUE 2016C 2017 120 CONTINUE 2018* 2019* INCREMENT OPCOUNT FOR FINDING SMALL SUBDIAGONAL ELEMENT. 2020 OPS = OPS + 3*( MIN(M,N-1)-L+1 ) 2021 IF (M .EQ. L) GO TO 210 2022 130 IF (J .EQ. 30) GO TO 1000 2023 J = J + 1 2024C .......... FORM SHIFT .......... 2025 L1 = L + 1 2026 S = DSQRT(E2(L)) 2027 G = D(L) 2028 P = (D(L1) - G) / (2.0D0 * S) 2029 R = PYTHAG(P,1.0D0) 2030 D(L) = S / (P + DSIGN(R,P)) 2031 H = G - D(L) 2032C 2033 DO 140 I = L1, N 2034 140 D(I) = D(I) - H 2035C 2036 F = F + H 2037* 2038* INCREMENT OPCOUNT FOR FORMING SHIFT AND SUBTRACTING. 2039 OPS = OPS + 8 + (I-L1+1) 2040C .......... RATIONAL QL TRANSFORMATION .......... 2041 G = D(M) 2042 IF (G .EQ. 0.0D0) G = B 2043 H = G 2044 S = 0.0D0 2045 MML = M - L 2046C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... 2047 DO 200 II = 1, MML 2048 I = M - II 2049 P = G * H 2050 R = P + E2(I) 2051 E2(I+1) = S * R 2052 S = E2(I) / R 2053 D(I+1) = H + S * (H + D(I)) 2054 G = D(I) - E2(I) / G 2055 IF (G .EQ. 0.0D0) G = B 2056 H = G * P / R 2057 200 CONTINUE 2058C 2059 E2(L) = S * G 2060 D(L) = H 2061* 2062* INCREMENT OPCOUNT FOR INNER LOOP. 2063 OPS = OPS + MML*11 + 1 2064* 2065* INCREMENT ITERATION COUNTER 2066 ITCNT = ITCNT + 1 2067C .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST .......... 2068 IF (H .EQ. 0.0D0) GO TO 210 2069 IF (DABS(E2(L)) .LE. DABS(C/H)) GO TO 210 2070 E2(L) = H * E2(L) 2071 IF (E2(L) .NE. 0.0D0) GO TO 130 2072 210 P = D(L) + F 2073C .......... ORDER EIGENVALUES .......... 2074 IF (L .EQ. 1) GO TO 250 2075C .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... 2076 DO 230 II = 2, L 2077 I = L + 2 - II 2078 IF (P .GE. D(I-1)) GO TO 270 2079 D(I) = D(I-1) 2080 230 CONTINUE 2081C 2082 250 I = 1 2083 270 D(I) = P 2084 290 CONTINUE 2085C 2086 GO TO 1001 2087C .......... SET ERROR -- NO CONVERGENCE TO AN 2088C EIGENVALUE AFTER 30 ITERATIONS .......... 2089 1000 IERR = L 2090 1001 CONTINUE 2091* 2092* COMPUTE FINAL OP COUNT 2093 OPS = OPS + OPST 2094 RETURN 2095 END 2096 SUBROUTINE TRED1(NM,N,A,D,E,E2) 2097C 2098 INTEGER I,J,K,L,N,II,NM,JP1 2099 DOUBLE PRECISION A(NM,N),D(N),E(N),E2(N) 2100 DOUBLE PRECISION F,G,H,SCALE 2101* 2102* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT. 2103* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED. 2104* .. COMMON BLOCKS .. 2105 COMMON / LATIME / OPS, ITCNT 2106* .. 2107* .. SCALARS IN COMMON .. 2108 DOUBLE PRECISION ITCNT, OPS 2109* .. 2110C 2111C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1, 2112C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. 2113C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). 2114C 2115C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX 2116C TO A SYMMETRIC TRIDIAGONAL MATRIX USING 2117C ORTHOGONAL SIMILARITY TRANSFORMATIONS. 2118C 2119C ON INPUT 2120C 2121C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 2122C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 2123C DIMENSION STATEMENT. 2124C 2125C N IS THE ORDER OF THE MATRIX. 2126C 2127C A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE 2128C LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. 2129C 2130C ON OUTPUT 2131C 2132C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- 2133C FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER 2134C TRIANGLE. THE FULL UPPER TRIANGLE OF A IS UNALTERED. 2135C 2136C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. 2137C 2138C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL 2139C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. 2140C 2141C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. 2142C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. 2143C 2144C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 2145C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 2146C 2147C THIS VERSION DATED AUGUST 1983. 2148C 2149C ------------------------------------------------------------------ 2150C 2151* 2152 OPS = OPS + MAX( 0.0D0, (4.0D0/3.0D0)*DBLE(N)**3 + 2153 $ 12.0D0*DBLE(N)**2 + 2154 $ (11.0D0/3.0D0)*N - 22 ) 2155* 2156 DO 100 I = 1, N 2157 D(I) = A(N,I) 2158 A(N,I) = A(I,I) 2159 100 CONTINUE 2160C .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... 2161 DO 300 II = 1, N 2162 I = N + 1 - II 2163 L = I - 1 2164 H = 0.0D0 2165 SCALE = 0.0D0 2166 IF (L .LT. 1) GO TO 130 2167C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... 2168 DO 120 K = 1, L 2169 120 SCALE = SCALE + DABS(D(K)) 2170C 2171 IF (SCALE .NE. 0.0D0) GO TO 140 2172C 2173 DO 125 J = 1, L 2174 D(J) = A(L,J) 2175 A(L,J) = A(I,J) 2176 A(I,J) = 0.0D0 2177 125 CONTINUE 2178C 2179 130 E(I) = 0.0D0 2180 E2(I) = 0.0D0 2181 GO TO 300 2182C 2183 140 DO 150 K = 1, L 2184 D(K) = D(K) / SCALE 2185 H = H + D(K) * D(K) 2186 150 CONTINUE 2187C 2188 E2(I) = SCALE * SCALE * H 2189 F = D(L) 2190 G = -DSIGN(DSQRT(H),F) 2191 E(I) = SCALE * G 2192 H = H - F * G 2193 D(L) = F - G 2194 IF (L .EQ. 1) GO TO 285 2195C .......... FORM A*U .......... 2196 DO 170 J = 1, L 2197 170 E(J) = 0.0D0 2198C 2199 DO 240 J = 1, L 2200 F = D(J) 2201 G = E(J) + A(J,J) * F 2202 JP1 = J + 1 2203 IF (L .LT. JP1) GO TO 220 2204C 2205 DO 200 K = JP1, L 2206 G = G + A(K,J) * D(K) 2207 E(K) = E(K) + A(K,J) * F 2208 200 CONTINUE 2209C 2210 220 E(J) = G 2211 240 CONTINUE 2212C .......... FORM P .......... 2213 F = 0.0D0 2214C 2215 DO 245 J = 1, L 2216 E(J) = E(J) / H 2217 F = F + E(J) * D(J) 2218 245 CONTINUE 2219C 2220 H = F / (H + H) 2221C .......... FORM Q .......... 2222 DO 250 J = 1, L 2223 250 E(J) = E(J) - H * D(J) 2224C .......... FORM REDUCED A .......... 2225 DO 280 J = 1, L 2226 F = D(J) 2227 G = E(J) 2228C 2229 DO 260 K = J, L 2230 260 A(K,J) = A(K,J) - F * E(K) - G * D(K) 2231C 2232 280 CONTINUE 2233C 2234 285 DO 290 J = 1, L 2235 F = D(J) 2236 D(J) = A(L,J) 2237 A(L,J) = A(I,J) 2238 A(I,J) = F * SCALE 2239 290 CONTINUE 2240C 2241 300 CONTINUE 2242C 2243 RETURN 2244 END 2245 SUBROUTINE BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,IERR,RV4,RV5) 2246* 2247* EISPACK ROUTINE. 2248* MODIFIED FOR COMPARISON WITH LAPACK ROUTINES. 2249* 2250* CONVERGENCE TEST WAS MODIFIED TO BE THE SAME AS IN DSTEBZ. 2251* 2252C 2253 INTEGER I,J,K,L,M,N,P,Q,R,S,II,MM,M1,M2,TAG,IERR,ISTURM 2254 DOUBLE PRECISION D(N),E(N),E2(N),W(MM),RV4(N),RV5(N) 2255 DOUBLE PRECISION U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,TST1,TST2,EPSLON 2256 INTEGER IND(MM) 2257* 2258* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 2259* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 2260* .. COMMON BLOCKS .. 2261 COMMON / LATIME / OPS, ITCNT 2262* .. 2263* .. SCALARS IN COMMON .. 2264 DOUBLE PRECISION ITCNT, OPS 2265* .. 2266C 2267C THIS SUBROUTINE IS A TRANSLATION OF THE BISECTION TECHNIQUE 2268C IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON. 2269C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). 2270C 2271C THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL 2272C SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL, 2273C USING BISECTION. 2274C 2275C ON INPUT 2276C 2277C N IS THE ORDER OF THE MATRIX. 2278C 2279C EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED 2280C EIGENVALUES. IF THE INPUT EPS1 IS NON-POSITIVE, 2281C IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE, 2282C NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE 2283C PRECISION AND THE 1-NORM OF THE SUBMATRIX. 2284C 2285C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. 2286C 2287C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX 2288C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. 2289C 2290C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. 2291C E2(1) IS ARBITRARY. 2292C 2293C LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES. 2294C IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND. 2295C 2296C MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF 2297C EIGENVALUES IN THE INTERVAL. WARNING. IF MORE THAN 2298C MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL, 2299C AN ERROR RETURN IS MADE WITH NO EIGENVALUES FOUND. 2300C 2301C ON OUTPUT 2302C 2303C EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS 2304C (LAST) DEFAULT VALUE. 2305C 2306C D AND E ARE UNALTERED. 2307C 2308C ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED 2309C AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE 2310C MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. 2311C E2(1) IS ALSO SET TO ZERO. 2312C 2313C M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB). 2314C 2315C W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER. 2316C 2317C IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES 2318C ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- 2319C 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM 2320C THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.. 2321C 2322C IERR IS SET TO 2323C ZERO FOR NORMAL RETURN, 2324C 3*N+1 IF M EXCEEDS MM. 2325C 2326C RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS. 2327C 2328C THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM 2329C APPEARS IN BISECT IN-LINE. 2330C 2331C NOTE THAT SUBROUTINE TQL1 OR IMTQL1 IS GENERALLY FASTER THAN 2332C BISECT, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND. 2333C 2334C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 2335C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 2336C 2337C THIS VERSION DATED AUGUST 1983. 2338C 2339C ------------------------------------------------------------------ 2340C 2341 DOUBLE PRECISION ONE 2342 PARAMETER ( ONE = 1.0D0 ) 2343 DOUBLE PRECISION RELFAC 2344 PARAMETER ( RELFAC = 2.0D0 ) 2345 DOUBLE PRECISION ATOLI, RTOLI, SAFEMN, TMP1, TMP2, TNORM, ULP 2346 DOUBLE PRECISION DLAMCH, PIVMIN 2347 EXTERNAL DLAMCH 2348* INITIALIZE ITERATION COUNT. 2349 ITCNT = 0 2350 SAFEMN = DLAMCH( 'S' ) 2351 ULP = DLAMCH( 'E' )*DLAMCH( 'B' ) 2352 RTOLI = ULP*RELFAC 2353 IERR = 0 2354 TAG = 0 2355 T1 = LB 2356 T2 = UB 2357C .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES .......... 2358 DO 40 I = 1, N 2359 IF (I .EQ. 1) GO TO 20 2360CCC TST1 = DABS(D(I)) + DABS(D(I-1)) 2361CCC TST2 = TST1 + DABS(E(I)) 2362CCC IF (TST2 .GT. TST1) GO TO 40 2363 TMP1 = E( I )**2 2364 IF( ABS( D(I)*D(I-1) )*ULP**2+SAFEMN.LE.TMP1 ) 2365 $ GO TO 40 2366 20 E2(I) = 0.0D0 2367 40 CONTINUE 2368* INCREMENT OPCOUNT FOR DETERMINING IF MATRIX SPLITS. 2369 OPS = OPS + 5*( N-1 ) 2370C 2371C COMPUTE QUANTITIES NEEDED FOR CONVERGENCE TEST. 2372 TMP1 = D( 1 ) - ABS( E( 2 ) ) 2373 TMP2 = D( 1 ) + ABS( E( 2 ) ) 2374 PIVMIN = ONE 2375 DO 41 I = 2, N - 1 2376 TMP1 = MIN( TMP1, D( I )-ABS( E( I ) )-ABS( E( I+1 ) ) ) 2377 TMP2 = MAX( TMP2, D( I )+ABS( E( I ) )+ABS( E( I+1 ) ) ) 2378 PIVMIN = MAX( PIVMIN, E( I )**2 ) 2379 41 CONTINUE 2380 TMP1 = MIN( TMP1, D( N )-ABS( E( N ) ) ) 2381 TMP2 = MAX( TMP2, D( N )+ABS( E( N ) ) ) 2382 PIVMIN = MAX( PIVMIN, E( N )**2 ) 2383 PIVMIN = PIVMIN*SAFEMN 2384 TNORM = MAX( ABS(TMP1), ABS(TMP2) ) 2385 ATOLI = ULP*TNORM 2386* INCREMENT OPCOUNT FOR COMPUTING THESE QUANTITIES. 2387 OPS = OPS + 4*( N-1 ) 2388C 2389C .......... DETERMINE THE NUMBER OF EIGENVALUES 2390C IN THE INTERVAL .......... 2391 P = 1 2392 Q = N 2393 X1 = UB 2394 ISTURM = 1 2395 GO TO 320 2396 60 M = S 2397 X1 = LB 2398 ISTURM = 2 2399 GO TO 320 2400 80 M = M - S 2401 IF (M .GT. MM) GO TO 980 2402 Q = 0 2403 R = 0 2404C .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING 2405C INTERVAL BY THE GERSCHGORIN BOUNDS .......... 2406 100 IF (R .EQ. M) GO TO 1001 2407 TAG = TAG + 1 2408 P = Q + 1 2409 XU = D(P) 2410 X0 = D(P) 2411 U = 0.0D0 2412C 2413 DO 120 Q = P, N 2414 X1 = U 2415 U = 0.0D0 2416 V = 0.0D0 2417 IF (Q .EQ. N) GO TO 110 2418 U = DABS(E(Q+1)) 2419 V = E2(Q+1) 2420 110 XU = DMIN1(D(Q)-(X1+U),XU) 2421 X0 = DMAX1(D(Q)+(X1+U),X0) 2422 IF (V .EQ. 0.0D0) GO TO 140 2423 120 CONTINUE 2424* INCREMENT OPCOUNT FOR REFINING INTERVAL. 2425 OPS = OPS + ( N-P+1 )*2 2426C 2427 140 X1 = EPSLON(DMAX1(DABS(XU),DABS(X0))) 2428 IF (EPS1 .LE. 0.0D0) EPS1 = -X1 2429 IF (P .NE. Q) GO TO 180 2430C .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL .......... 2431 IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940 2432 M1 = P 2433 M2 = P 2434 RV5(P) = D(P) 2435 GO TO 900 2436 180 X1 = X1 * (Q - P + 1) 2437 LB = DMAX1(T1,XU-X1) 2438 UB = DMIN1(T2,X0+X1) 2439 X1 = LB 2440 ISTURM = 3 2441 GO TO 320 2442 200 M1 = S + 1 2443 X1 = UB 2444 ISTURM = 4 2445 GO TO 320 2446 220 M2 = S 2447 IF (M1 .GT. M2) GO TO 940 2448C .......... FIND ROOTS BY BISECTION .......... 2449 X0 = UB 2450 ISTURM = 5 2451C 2452 DO 240 I = M1, M2 2453 RV5(I) = UB 2454 RV4(I) = LB 2455 240 CONTINUE 2456C .......... LOOP FOR K-TH EIGENVALUE 2457C FOR K=M2 STEP -1 UNTIL M1 DO -- 2458C (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) .......... 2459 K = M2 2460 250 XU = LB 2461C .......... FOR I=K STEP -1 UNTIL M1 DO -- .......... 2462 DO 260 II = M1, K 2463 I = M1 + K - II 2464 IF (XU .GE. RV4(I)) GO TO 260 2465 XU = RV4(I) 2466 GO TO 280 2467 260 CONTINUE 2468C 2469 280 IF (X0 .GT. RV5(K)) X0 = RV5(K) 2470C .......... NEXT BISECTION STEP .......... 2471 300 X1 = (XU + X0) * 0.5D0 2472CCC IF ((X0 - XU) .LE. DABS(EPS1)) GO TO 420 2473CCC TST1 = 2.0D0 * (DABS(XU) + DABS(X0)) 2474CCC TST2 = TST1 + (X0 - XU) 2475CCC IF (TST2 .EQ. TST1) GO TO 420 2476 TMP1 = ABS( X0 - XU ) 2477 TMP2 = MAX( ABS( X0 ), ABS( XU ) ) 2478 IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) 2479 $ GO TO 420 2480C .......... IN-LINE PROCEDURE FOR STURM SEQUENCE .......... 2481 320 S = P - 1 2482 U = 1.0D0 2483C 2484 DO 340 I = P, Q 2485 IF (U .NE. 0.0D0) GO TO 325 2486 V = DABS(E(I)) / EPSLON(1.0D0) 2487 IF (E2(I) .EQ. 0.0D0) V = 0.0D0 2488 GO TO 330 2489 325 V = E2(I) / U 2490 330 U = D(I) - X1 - V 2491 IF (U .LT. 0.0D0) S = S + 1 2492 340 CONTINUE 2493* INCREMENT OPCOUNT FOR STURM SEQUENCE. 2494 OPS = OPS + ( Q-P+1 )*3 2495* INCREMENT ITERATION COUNTER. 2496 ITCNT = ITCNT + 1 2497C 2498 GO TO (60,80,200,220,360), ISTURM 2499C .......... REFINE INTERVALS .......... 2500 360 IF (S .GE. K) GO TO 400 2501 XU = X1 2502 IF (S .GE. M1) GO TO 380 2503 RV4(M1) = X1 2504 GO TO 300 2505 380 RV4(S+1) = X1 2506 IF (RV5(S) .GT. X1) RV5(S) = X1 2507 GO TO 300 2508 400 X0 = X1 2509 GO TO 300 2510C .......... K-TH EIGENVALUE FOUND .......... 2511 420 RV5(K) = X1 2512 K = K - 1 2513 IF (K .GE. M1) GO TO 250 2514C .......... ORDER EIGENVALUES TAGGED WITH THEIR 2515C SUBMATRIX ASSOCIATIONS .......... 2516 900 S = R 2517 R = R + M2 - M1 + 1 2518 J = 1 2519 K = M1 2520C 2521 DO 920 L = 1, R 2522 IF (J .GT. S) GO TO 910 2523 IF (K .GT. M2) GO TO 940 2524 IF (RV5(K) .GE. W(L)) GO TO 915 2525C 2526 DO 905 II = J, S 2527 I = L + S - II 2528 W(I+1) = W(I) 2529 IND(I+1) = IND(I) 2530 905 CONTINUE 2531C 2532 910 W(L) = RV5(K) 2533 IND(L) = TAG 2534 K = K + 1 2535 GO TO 920 2536 915 J = J + 1 2537 920 CONTINUE 2538C 2539 940 IF (Q .LT. N) GO TO 100 2540 GO TO 1001 2541C .......... SET ERROR -- UNDERESTIMATE OF NUMBER OF 2542C EIGENVALUES IN INTERVAL .......... 2543 980 IERR = 3 * N + 1 2544 1001 LB = T1 2545 UB = T2 2546 RETURN 2547 END 2548 SUBROUTINE TINVIT(NM,N,D,E,E2,M,W,IND,Z, 2549 X IERR,RV1,RV2,RV3,RV4,RV6) 2550* 2551* EISPACK ROUTINE. 2552* 2553* CONVERGENCE TEST WAS NOT MODIFIED, SINCE IT SHOULD GIVE 2554* APPROXIMATELY THE SAME LEVEL OF ACCURACY AS LAPACK ROUTINE, 2555* ALTHOUGH THE EIGENVECTORS MAY NOT BE AS CLOSE TO ORTHOGONAL. 2556* 2557C 2558 INTEGER I,J,M,N,P,Q,R,S,II,IP,JJ,NM,ITS,TAG,IERR,GROUP 2559 DOUBLE PRECISION D(N),E(N),E2(N),W(M),Z(NM,M), 2560 X RV1(N),RV2(N),RV3(N),RV4(N),RV6(N) 2561 DOUBLE PRECISION U,V,UK,XU,X0,X1,EPS2,EPS3,EPS4,NORM,ORDER,EPSLON, 2562 X PYTHAG 2563 INTEGER IND(M) 2564* 2565* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 2566* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 2567* .. COMMON BLOCKS .. 2568 COMMON / LATIME / OPS, ITCNT 2569 COMMON / PYTHOP / OPST 2570* .. 2571* .. SCALARS IN COMMON .. 2572 DOUBLE PRECISION ITCNT, OPS, OPST 2573* .. 2574C 2575C THIS SUBROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH- 2576C NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON. 2577C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). 2578C 2579C THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL 2580C SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, 2581C USING INVERSE ITERATION. 2582C 2583C ON INPUT 2584C 2585C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 2586C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 2587C DIMENSION STATEMENT. 2588C 2589C N IS THE ORDER OF THE MATRIX. 2590C 2591C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. 2592C 2593C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX 2594C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. 2595C 2596C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E, 2597C WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E. 2598C E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN 2599C THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM 2600C OF THE MAGNITUDES OF D(I) AND D(I-1). E2(1) MUST CONTAIN 2601C 0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0D0 2602C IF THE EIGENVALUES ARE IN DESCENDING ORDER. IF BISECT, 2603C TRIDIB, OR IMTQLV HAS BEEN USED TO FIND THE EIGENVALUES, 2604C THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE. 2605C 2606C M IS THE NUMBER OF SPECIFIED EIGENVALUES. 2607C 2608C W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER. 2609C 2610C IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES 2611C ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- 2612C 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM 2613C THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC. 2614C 2615C ON OUTPUT 2616C 2617C ALL INPUT ARRAYS ARE UNALTERED. 2618C 2619C Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS. 2620C ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO. 2621C 2622C IERR IS SET TO 2623C ZERO FOR NORMAL RETURN, 2624C -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH 2625C EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS. 2626C 2627C RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS. 2628C 2629C CALLS PYTHAG FOR DSQRT(A*A + B*B) . 2630C 2631C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 2632C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 2633C 2634C THIS VERSION DATED AUGUST 1983. 2635C 2636C ------------------------------------------------------------------ 2637C 2638* INITIALIZE ITERATION COUNT. 2639 ITCNT = 0 2640 IERR = 0 2641 IF (M .EQ. 0) GO TO 1001 2642 TAG = 0 2643 ORDER = 1.0D0 - E2(1) 2644 Q = 0 2645C .......... ESTABLISH AND PROCESS NEXT SUBMATRIX .......... 2646 100 P = Q + 1 2647C 2648 DO 120 Q = P, N 2649 IF (Q .EQ. N) GO TO 140 2650 IF (E2(Q+1) .EQ. 0.0D0) GO TO 140 2651 120 CONTINUE 2652C .......... FIND VECTORS BY INVERSE ITERATION .......... 2653 140 TAG = TAG + 1 2654 S = 0 2655C 2656 DO 920 R = 1, M 2657 IF (IND(R) .NE. TAG) GO TO 920 2658 ITS = 1 2659 X1 = W(R) 2660 IF (S .NE. 0) GO TO 510 2661C .......... CHECK FOR ISOLATED ROOT .......... 2662 XU = 1.0D0 2663 IF (P .NE. Q) GO TO 490 2664 RV6(P) = 1.0D0 2665 GO TO 870 2666 490 NORM = DABS(D(P)) 2667 IP = P + 1 2668C 2669 DO 500 I = IP, Q 2670 500 NORM = DMAX1(NORM, DABS(D(I))+DABS(E(I))) 2671C .......... EPS2 IS THE CRITERION FOR GROUPING, 2672C EPS3 REPLACES ZERO PIVOTS AND EQUAL 2673C ROOTS ARE MODIFIED BY EPS3, 2674C EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .......... 2675 EPS2 = 1.0D-3 * NORM 2676 EPS3 = EPSLON(NORM) 2677 UK = Q - P + 1 2678 EPS4 = UK * EPS3 2679 UK = EPS4 / DSQRT(UK) 2680* INCREMENT OPCOUNT FOR COMPUTING CRITERIA. 2681 OPS = OPS + ( Q-IP+4 ) 2682 S = P 2683 505 GROUP = 0 2684 GO TO 520 2685C .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... 2686 510 IF (DABS(X1-X0) .GE. EPS2) GO TO 505 2687 GROUP = GROUP + 1 2688 IF (ORDER * (X1 - X0) .LE. 0.0D0) X1 = X0 + ORDER * EPS3 2689C .......... ELIMINATION WITH INTERCHANGES AND 2690C INITIALIZATION OF VECTOR .......... 2691 520 V = 0.0D0 2692C 2693 DO 580 I = P, Q 2694 RV6(I) = UK 2695 IF (I .EQ. P) GO TO 560 2696 IF (DABS(E(I)) .LT. DABS(U)) GO TO 540 2697C .......... WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF 2698C E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY .......... 2699 XU = U / E(I) 2700 RV4(I) = XU 2701 RV1(I-1) = E(I) 2702 RV2(I-1) = D(I) - X1 2703 RV3(I-1) = 0.0D0 2704 IF (I .NE. Q) RV3(I-1) = E(I+1) 2705 U = V - XU * RV2(I-1) 2706 V = -XU * RV3(I-1) 2707 GO TO 580 2708 540 XU = E(I) / U 2709 RV4(I) = XU 2710 RV1(I-1) = U 2711 RV2(I-1) = V 2712 RV3(I-1) = 0.0D0 2713 560 U = D(I) - X1 - XU * V 2714 IF (I .NE. Q) V = E(I+1) 2715 580 CONTINUE 2716* INCREMENT OPCOUNT FOR ELIMINATION. 2717 OPS = OPS + ( Q-P+1 )*5 2718C 2719 IF (U .EQ. 0.0D0) U = EPS3 2720 RV1(Q) = U 2721 RV2(Q) = 0.0D0 2722 RV3(Q) = 0.0D0 2723C .......... BACK SUBSTITUTION 2724C FOR I=Q STEP -1 UNTIL P DO -- .......... 2725 600 DO 620 II = P, Q 2726 I = P + Q - II 2727 RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I) 2728 V = U 2729 U = RV6(I) 2730 620 CONTINUE 2731* INCREMENT OPCOUNT FOR BACK SUBSTITUTION. 2732 OPS = OPS + ( Q-P+1 )*5 2733C .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS 2734C MEMBERS OF GROUP .......... 2735 IF (GROUP .EQ. 0) GO TO 700 2736 J = R 2737C 2738 DO 680 JJ = 1, GROUP 2739 630 J = J - 1 2740 IF (IND(J) .NE. TAG) GO TO 630 2741 XU = 0.0D0 2742C 2743 DO 640 I = P, Q 2744 640 XU = XU + RV6(I) * Z(I,J) 2745C 2746 DO 660 I = P, Q 2747 660 RV6(I) = RV6(I) - XU * Z(I,J) 2748C 2749* INCREMENT OPCOUNT FOR ORTHOGONALIZING. 2750 OPS = OPS + ( Q-P+1 )*4 2751 680 CONTINUE 2752C 2753 700 NORM = 0.0D0 2754C 2755 DO 720 I = P, Q 2756 720 NORM = NORM + DABS(RV6(I)) 2757* INCREMENT OPCOUNT FOR COMPUTING NORM. 2758 OPS = OPS + ( Q-P+1 ) 2759C 2760 IF (NORM .GE. 1.0D0) GO TO 840 2761C .......... FORWARD SUBSTITUTION .......... 2762 IF (ITS .EQ. 5) GO TO 830 2763 IF (NORM .NE. 0.0D0) GO TO 740 2764 RV6(S) = EPS4 2765 S = S + 1 2766 IF (S .GT. Q) S = P 2767 GO TO 780 2768 740 XU = EPS4 / NORM 2769C 2770 DO 760 I = P, Q 2771 760 RV6(I) = RV6(I) * XU 2772C .......... ELIMINATION OPERATIONS ON NEXT VECTOR 2773C ITERATE .......... 2774 780 DO 820 I = IP, Q 2775 U = RV6(I) 2776C .......... IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE 2777C WAS PERFORMED EARLIER IN THE 2778C TRIANGULARIZATION PROCESS .......... 2779 IF (RV1(I-1) .NE. E(I)) GO TO 800 2780 U = RV6(I-1) 2781 RV6(I-1) = RV6(I) 2782 800 RV6(I) = U - RV4(I) * RV6(I-1) 2783 820 CONTINUE 2784* INCREMENT OPCOUNT FOR FORWARD SUBSTITUTION. 2785 OPS = OPS + ( Q-P+1 ) + ( Q-IP+1 )*2 2786C 2787 ITS = ITS + 1 2788* INCREMENT ITERATION COUNTER. 2789 ITCNT = ITCNT + 1 2790 GO TO 600 2791C .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... 2792 830 IERR = -R 2793 XU = 0.0D0 2794 GO TO 870 2795C .......... NORMALIZE SO THAT SUM OF SQUARES IS 2796C 1 AND EXPAND TO FULL ORDER .......... 2797 840 U = 0.0D0 2798C 2799 DO 860 I = P, Q 2800 860 U = PYTHAG(U,RV6(I)) 2801C 2802 XU = 1.0D0 / U 2803C 2804 870 DO 880 I = 1, N 2805 880 Z(I,R) = 0.0D0 2806C 2807 DO 900 I = P, Q 2808 900 Z(I,R) = RV6(I) * XU 2809* INCREMENT OPCOUNT FOR NORMALIZING. 2810 OPS = OPS + ( Q-P+1 ) 2811C 2812 X0 = X1 2813 920 CONTINUE 2814C 2815 IF (Q .LT. N) GO TO 100 2816* INCREMENT OPCOUNT FOR USE OF FUNCTION PYTHAG. 2817 OPS = OPS + OPST 2818 1001 RETURN 2819 END 2820 SUBROUTINE TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,M,W,IND,IERR,RV4,RV5) 2821* 2822* EISPACK ROUTINE. 2823* MODIFIED FOR COMPARISON WITH LAPACK ROUTINES. 2824* 2825* CONVERGENCE TEST WAS MODIFIED TO BE THE SAME AS IN DSTEBZ. 2826* 2827C 2828 INTEGER I,J,K,L,M,N,P,Q,R,S,II,M1,M2,M11,M22,TAG,IERR,ISTURM 2829 DOUBLE PRECISION D(N),E(N),E2(N),W(M),RV4(N),RV5(N) 2830 DOUBLE PRECISION U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,TST1,TST2,EPSLON 2831 INTEGER IND(M) 2832* 2833* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 2834* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 2835* .. COMMON BLOCKS .. 2836 COMMON / LATIME / OPS, ITCNT 2837* .. 2838* .. SCALARS IN COMMON .. 2839 DOUBLE PRECISION ITCNT, OPS 2840* .. 2841C 2842C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BISECT, 2843C NUM. MATH. 9, 386-393(1967) BY BARTH, MARTIN, AND WILKINSON. 2844C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 249-256(1971). 2845C 2846C THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL 2847C SYMMETRIC MATRIX BETWEEN SPECIFIED BOUNDARY INDICES, 2848C USING BISECTION. 2849C 2850C ON INPUT 2851C 2852C N IS THE ORDER OF THE MATRIX. 2853C 2854C EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED 2855C EIGENVALUES. IF THE INPUT EPS1 IS NON-POSITIVE, 2856C IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE, 2857C NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE 2858C PRECISION AND THE 1-NORM OF THE SUBMATRIX. 2859C 2860C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. 2861C 2862C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX 2863C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. 2864C 2865C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. 2866C E2(1) IS ARBITRARY. 2867C 2868C M11 SPECIFIES THE LOWER BOUNDARY INDEX FOR THE DESIRED 2869C EIGENVALUES. 2870C 2871C M SPECIFIES THE NUMBER OF EIGENVALUES DESIRED. THE UPPER 2872C BOUNDARY INDEX M22 IS THEN OBTAINED AS M22=M11+M-1. 2873C 2874C ON OUTPUT 2875C 2876C EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS 2877C (LAST) DEFAULT VALUE. 2878C 2879C D AND E ARE UNALTERED. 2880C 2881C ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED 2882C AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE 2883C MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. 2884C E2(1) IS ALSO SET TO ZERO. 2885C 2886C LB AND UB DEFINE AN INTERVAL CONTAINING EXACTLY THE DESIRED 2887C EIGENVALUES. 2888C 2889C W CONTAINS, IN ITS FIRST M POSITIONS, THE EIGENVALUES 2890C BETWEEN INDICES M11 AND M22 IN ASCENDING ORDER. 2891C 2892C IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES 2893C ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- 2894C 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM 2895C THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.. 2896C 2897C IERR IS SET TO 2898C ZERO FOR NORMAL RETURN, 2899C 3*N+1 IF MULTIPLE EIGENVALUES AT INDEX M11 MAKE 2900C UNIQUE SELECTION IMPOSSIBLE, 2901C 3*N+2 IF MULTIPLE EIGENVALUES AT INDEX M22 MAKE 2902C UNIQUE SELECTION IMPOSSIBLE. 2903C 2904C RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS. 2905C 2906C NOTE THAT SUBROUTINE TQL1, IMTQL1, OR TQLRAT IS GENERALLY FASTER 2907C THAN TRIDIB, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND. 2908C 2909C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 2910C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 2911C 2912C THIS VERSION DATED AUGUST 1983. 2913C 2914C ------------------------------------------------------------------ 2915C 2916 DOUBLE PRECISION ONE 2917 PARAMETER ( ONE = 1.0D0 ) 2918 DOUBLE PRECISION RELFAC 2919 PARAMETER ( RELFAC = 2.0D0 ) 2920 DOUBLE PRECISION ATOLI, RTOLI, SAFEMN, TMP1, TMP2, TNORM, ULP 2921 DOUBLE PRECISION DLAMCH, PIVMIN 2922 EXTERNAL DLAMCH 2923* INITIALIZE ITERATION COUNT. 2924 ITCNT = 0 2925 SAFEMN = DLAMCH( 'S' ) 2926 ULP = DLAMCH( 'E' )*DLAMCH( 'B' ) 2927 RTOLI = ULP*RELFAC 2928 IERR = 0 2929 TAG = 0 2930 XU = D(1) 2931 X0 = D(1) 2932 U = 0.0D0 2933C .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES AND DETERMINE AN 2934C INTERVAL CONTAINING ALL THE EIGENVALUES .......... 2935 PIVMIN = ONE 2936 DO 40 I = 1, N 2937 X1 = U 2938 U = 0.0D0 2939 IF (I .NE. N) U = DABS(E(I+1)) 2940 XU = DMIN1(D(I)-(X1+U),XU) 2941 X0 = DMAX1(D(I)+(X1+U),X0) 2942 IF (I .EQ. 1) GO TO 20 2943CCC TST1 = DABS(D(I)) + DABS(D(I-1)) 2944CCC TST2 = TST1 + DABS(E(I)) 2945CCC IF (TST2 .GT. TST1) GO TO 40 2946 TMP1 = E( I )**2 2947 IF( ABS( D(I)*D(I-1) )*ULP**2+SAFEMN.LE.TMP1 ) THEN 2948 PIVMIN = MAX( PIVMIN, TMP1 ) 2949 GO TO 40 2950 END IF 2951 20 E2(I) = 0.0D0 2952 40 CONTINUE 2953 PIVMIN = PIVMIN*SAFEMN 2954 TNORM = MAX( ABS( XU ), ABS( X0 ) ) 2955 ATOLI = ULP*TNORM 2956* INCREMENT OPCOUNT FOR DETERMINING IF MATRIX SPLITS. 2957 OPS = OPS + 9*( N-1 ) 2958C 2959 X1 = N 2960 X1 = X1 * EPSLON(DMAX1(DABS(XU),DABS(X0))) 2961 XU = XU - X1 2962 T1 = XU 2963 X0 = X0 + X1 2964 T2 = X0 2965C .......... DETERMINE AN INTERVAL CONTAINING EXACTLY 2966C THE DESIRED EIGENVALUES .......... 2967 P = 1 2968 Q = N 2969 M1 = M11 - 1 2970 IF (M1 .EQ. 0) GO TO 75 2971 ISTURM = 1 2972 50 V = X1 2973 X1 = XU + (X0 - XU) * 0.5D0 2974 IF (X1 .EQ. V) GO TO 980 2975 GO TO 320 2976 60 IF (S - M1) 65, 73, 70 2977 65 XU = X1 2978 GO TO 50 2979 70 X0 = X1 2980 GO TO 50 2981 73 XU = X1 2982 T1 = X1 2983 75 M22 = M1 + M 2984 IF (M22 .EQ. N) GO TO 90 2985 X0 = T2 2986 ISTURM = 2 2987 GO TO 50 2988 80 IF (S - M22) 65, 85, 70 2989 85 T2 = X1 2990 90 Q = 0 2991 R = 0 2992C .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING 2993C INTERVAL BY THE GERSCHGORIN BOUNDS .......... 2994 100 IF (R .EQ. M) GO TO 1001 2995 TAG = TAG + 1 2996 P = Q + 1 2997 XU = D(P) 2998 X0 = D(P) 2999 U = 0.0D0 3000C 3001 DO 120 Q = P, N 3002 X1 = U 3003 U = 0.0D0 3004 V = 0.0D0 3005 IF (Q .EQ. N) GO TO 110 3006 U = DABS(E(Q+1)) 3007 V = E2(Q+1) 3008 110 XU = DMIN1(D(Q)-(X1+U),XU) 3009 X0 = DMAX1(D(Q)+(X1+U),X0) 3010 IF (V .EQ. 0.0D0) GO TO 140 3011 120 CONTINUE 3012* INCREMENT OPCOUNT FOR REFINING INTERVAL. 3013 OPS = OPS + ( N-P+1 )*2 3014C 3015 140 X1 = EPSLON(DMAX1(DABS(XU),DABS(X0))) 3016 IF (EPS1 .LE. 0.0D0) EPS1 = -X1 3017 IF (P .NE. Q) GO TO 180 3018C .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL .......... 3019 IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940 3020 M1 = P 3021 M2 = P 3022 RV5(P) = D(P) 3023 GO TO 900 3024 180 X1 = X1 * (Q - P + 1) 3025 LB = DMAX1(T1,XU-X1) 3026 UB = DMIN1(T2,X0+X1) 3027 X1 = LB 3028 ISTURM = 3 3029 GO TO 320 3030 200 M1 = S + 1 3031 X1 = UB 3032 ISTURM = 4 3033 GO TO 320 3034 220 M2 = S 3035 IF (M1 .GT. M2) GO TO 940 3036C .......... FIND ROOTS BY BISECTION .......... 3037 X0 = UB 3038 ISTURM = 5 3039C 3040 DO 240 I = M1, M2 3041 RV5(I) = UB 3042 RV4(I) = LB 3043 240 CONTINUE 3044C .......... LOOP FOR K-TH EIGENVALUE 3045C FOR K=M2 STEP -1 UNTIL M1 DO -- 3046C (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) .......... 3047 K = M2 3048 250 XU = LB 3049C .......... FOR I=K STEP -1 UNTIL M1 DO -- .......... 3050 DO 260 II = M1, K 3051 I = M1 + K - II 3052 IF (XU .GE. RV4(I)) GO TO 260 3053 XU = RV4(I) 3054 GO TO 280 3055 260 CONTINUE 3056C 3057 280 IF (X0 .GT. RV5(K)) X0 = RV5(K) 3058C .......... NEXT BISECTION STEP .......... 3059 300 X1 = (XU + X0) * 0.5D0 3060CCC IF ((X0 - XU) .LE. DABS(EPS1)) GO TO 420 3061CCC TST1 = 2.0D0 * (DABS(XU) + DABS(X0)) 3062CCC TST2 = TST1 + (X0 - XU) 3063CCC IF (TST2 .EQ. TST1) GO TO 420 3064 TMP1 = ABS( X0 - XU ) 3065 TMP2 = MAX( ABS( X0 ), ABS( XU ) ) 3066 IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) 3067 $ GO TO 420 3068C .......... IN-LINE PROCEDURE FOR STURM SEQUENCE .......... 3069 320 S = P - 1 3070 U = 1.0D0 3071C 3072 DO 340 I = P, Q 3073 IF (U .NE. 0.0D0) GO TO 325 3074 V = DABS(E(I)) / EPSLON(1.0D0) 3075 IF (E2(I) .EQ. 0.0D0) V = 0.0D0 3076 GO TO 330 3077 325 V = E2(I) / U 3078 330 U = D(I) - X1 - V 3079 IF (U .LT. 0.0D0) S = S + 1 3080 340 CONTINUE 3081* INCREMENT OPCOUNT FOR STURM SEQUENCE. 3082 OPS = OPS + ( Q-P+1 )*3 3083* INCREMENT ITERATION COUNTER. 3084 ITCNT = ITCNT + 1 3085C 3086 GO TO (60,80,200,220,360), ISTURM 3087C .......... REFINE INTERVALS .......... 3088 360 IF (S .GE. K) GO TO 400 3089 XU = X1 3090 IF (S .GE. M1) GO TO 380 3091 RV4(M1) = X1 3092 GO TO 300 3093 380 RV4(S+1) = X1 3094 IF (RV5(S) .GT. X1) RV5(S) = X1 3095 GO TO 300 3096 400 X0 = X1 3097 GO TO 300 3098C .......... K-TH EIGENVALUE FOUND .......... 3099 420 RV5(K) = X1 3100 K = K - 1 3101 IF (K .GE. M1) GO TO 250 3102C .......... ORDER EIGENVALUES TAGGED WITH THEIR 3103C SUBMATRIX ASSOCIATIONS .......... 3104 900 S = R 3105 R = R + M2 - M1 + 1 3106 J = 1 3107 K = M1 3108C 3109 DO 920 L = 1, R 3110 IF (J .GT. S) GO TO 910 3111 IF (K .GT. M2) GO TO 940 3112 IF (RV5(K) .GE. W(L)) GO TO 915 3113C 3114 DO 905 II = J, S 3115 I = L + S - II 3116 W(I+1) = W(I) 3117 IND(I+1) = IND(I) 3118 905 CONTINUE 3119C 3120 910 W(L) = RV5(K) 3121 IND(L) = TAG 3122 K = K + 1 3123 GO TO 920 3124 915 J = J + 1 3125 920 CONTINUE 3126C 3127 940 IF (Q .LT. N) GO TO 100 3128 GO TO 1001 3129C .......... SET ERROR -- INTERVAL CANNOT BE FOUND CONTAINING 3130C EXACTLY THE DESIRED EIGENVALUES .......... 3131 980 IERR = 3 * N + ISTURM 3132 1001 LB = T1 3133 UB = T2 3134 RETURN 3135 END 3136 SUBROUTINE DSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO) 3137 INTEGER LDX,N,P,LDU,LDV,JOB,INFO 3138 DOUBLE PRECISION X(LDX,*),S(*),E(*),U(LDU,*),V(LDV,*),WORK(*) 3139* 3140* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 3141* ITCNT IS INITIALIZED TO 0, IOPS IS ONLY INCREMENTED 3142* IOPST IS USED TO ACCUMULATE SMALL CONTRIBUTIONS TO IOPS 3143* TO AVOID ROUNDOFF ERROR 3144* .. COMMON BLOCKS .. 3145 COMMON /LATIME/ IOPS, ITCNT 3146* .. 3147* .. SCALARS IN COMMON .. 3148 DOUBLE PRECISION IOPS, ITCNT, IOPST 3149* .. 3150C 3151C 3152C DSVDC IS A SUBROUTINE TO REDUCE A DOUBLE PRECISION NXP MATRIX X 3153C BY ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM. THE 3154C DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X. THE 3155C COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS, 3156C AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS. 3157C 3158C ON ENTRY 3159C 3160C X DOUBLE PRECISION(LDX,P), WHERE LDX.GE.N. 3161C X CONTAINS THE MATRIX WHOSE SINGULAR VALUE 3162C DECOMPOSITION IS TO BE COMPUTED. X IS 3163C DESTROYED BY DSVDC. 3164C 3165C LDX INTEGER. 3166C LDX IS THE LEADING DIMENSION OF THE ARRAY X. 3167C 3168C N INTEGER. 3169C N IS THE NUMBER OF ROWS OF THE MATRIX X. 3170C 3171C P INTEGER. 3172C P IS THE NUMBER OF COLUMNS OF THE MATRIX X. 3173C 3174C LDU INTEGER. 3175C LDU IS THE LEADING DIMENSION OF THE ARRAY U. 3176C (SEE BELOW). 3177C 3178C LDV INTEGER. 3179C LDV IS THE LEADING DIMENSION OF THE ARRAY V. 3180C (SEE BELOW). 3181C 3182C WORK DOUBLE PRECISION(N). 3183C WORK IS A SCRATCH ARRAY. 3184C 3185C JOB INTEGER. 3186C JOB CONTROLS THE COMPUTATION OF THE SINGULAR 3187C VECTORS. IT HAS THE DECIMAL EXPANSION AB 3188C WITH THE FOLLOWING MEANING 3189C 3190C A.EQ.0 DO NOT COMPUTE THE LEFT SINGULAR 3191C VECTORS. 3192C A.EQ.1 RETURN THE N LEFT SINGULAR VECTORS 3193C IN U. 3194C A.GE.2 RETURN THE FIRST MIN(N,P) SINGULAR 3195C VECTORS IN U. 3196C B.EQ.0 DO NOT COMPUTE THE RIGHT SINGULAR 3197C VECTORS. 3198C B.EQ.1 RETURN THE RIGHT SINGULAR VECTORS 3199C IN V. 3200C 3201C ON RETURN 3202C 3203C S DOUBLE PRECISION(MM), WHERE MM=MIN(N+1,P). 3204C THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE 3205C SINGULAR VALUES OF X ARRANGED IN DESCENDING 3206C ORDER OF MAGNITUDE. 3207C 3208C E DOUBLE PRECISION(P), 3209C E ORDINARILY CONTAINS ZEROS. HOWEVER SEE THE 3210C DISCUSSION OF INFO FOR EXCEPTIONS. 3211C 3212C U DOUBLE PRECISION(LDU,K), WHERE LDU.GE.N. IF 3213C JOBA.EQ.1 THEN K.EQ.N, IF JOBA.GE.2 3214C THEN K.EQ.MIN(N,P). 3215C U CONTAINS THE MATRIX OF LEFT SINGULAR VECTORS. 3216C U IS NOT REFERENCED IF JOBA.EQ.0. IF N.LE.P 3217C OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X 3218C IN THE SUBROUTINE CALL. 3219C 3220C V DOUBLE PRECISION(LDV,P), WHERE LDV.GE.P. 3221C V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS. 3222C V IS NOT REFERENCED IF JOB.EQ.0. IF P.LE.N, 3223C THEN V MAY BE IDENTIFIED WITH X IN THE 3224C SUBROUTINE CALL. 3225C 3226C INFO INTEGER. 3227C THE SINGULAR VALUES (AND THEIR CORRESPONDING 3228C SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M) 3229C ARE CORRECT (HERE M=MIN(N,P)). THUS IF 3230C INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR 3231C VECTORS ARE CORRECT. IN ANY EVENT, THE MATRIX 3232C B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX 3233C WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE 3234C ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U) 3235C IS THE TRANSPOSE OF U). THUS THE SINGULAR 3236C VALUES OF X AND B ARE THE SAME. 3237C 3238C LINPACK. THIS VERSION DATED 08/14/78 . 3239C CORRECTION MADE TO SHIFT 2/84. 3240C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. 3241C 3242C DSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. 3243C 3244C EXTERNAL DROT 3245C BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2,DROTG 3246C FORTRAN DABS,DMAX1,MAX0,MIN0,MOD,DSQRT 3247C 3248C INTERNAL VARIABLES 3249C 3250 INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT, 3251 * MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1 3252 DOUBLE PRECISION DDOT,T 3253 DOUBLE PRECISION B,C,CS,EL,EMM1,F,G,DNRM2,SCALE,SHIFT,SL,SM,SN, 3254 * SMM1,T1,TEST 3255* DOUBLE PRECISION ZTEST,R 3256 LOGICAL WANTU,WANTV 3257* 3258* GET EPS FROM DLAMCH FOR NEW STOPPING CRITERION 3259 EXTERNAL DLAMCH 3260 DOUBLE PRECISION DLAMCH, EPS 3261 IF (N.LE.0 .OR. P.LE.0) RETURN 3262 EPS = DLAMCH( 'EPSILON' ) 3263* 3264C 3265C 3266C SET THE MAXIMUM NUMBER OF ITERATIONS. 3267C 3268 MAXIT = 50 3269C 3270C DETERMINE WHAT IS TO BE COMPUTED. 3271C 3272 WANTU = .FALSE. 3273 WANTV = .FALSE. 3274 JOBU = MOD(JOB,100)/10 3275 NCU = N 3276 IF (JOBU .GT. 1) NCU = MIN0(N,P) 3277 IF (JOBU .NE. 0) WANTU = .TRUE. 3278 IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE. 3279C 3280C REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS 3281C IN S AND THE SUPER-DIAGONAL ELEMENTS IN E. 3282C 3283* 3284* INITIALIZE OP COUNT 3285 IOPST = 0 3286 INFO = 0 3287 NCT = MIN0(N-1,P) 3288 NRT = MAX0(0,MIN0(P-2,N)) 3289 LU = MAX0(NCT,NRT) 3290 IF (LU .LT. 1) GO TO 170 3291 DO 160 L = 1, LU 3292 LP1 = L + 1 3293 IF (L .GT. NCT) GO TO 20 3294C 3295C COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND 3296C PLACE THE L-TH DIAGONAL IN S(L). 3297C 3298* 3299* INCREMENT OP COUNT 3300 IOPS = IOPS + (2*(N-L+1)+1) 3301 S(L) = DNRM2(N-L+1,X(L,L),1) 3302 IF (S(L) .EQ. 0.0D0) GO TO 10 3303 IF (X(L,L) .NE. 0.0D0) S(L) = DSIGN(S(L),X(L,L)) 3304* 3305* INCREMENT OP COUNT 3306 IOPS = IOPS + (N-L+3) 3307 CALL DSCAL(N-L+1,1.0D0/S(L),X(L,L),1) 3308 X(L,L) = 1.0D0 + X(L,L) 3309 10 CONTINUE 3310 S(L) = -S(L) 3311 20 CONTINUE 3312 IF (P .LT. LP1) GO TO 50 3313 DO 40 J = LP1, P 3314 IF (L .GT. NCT) GO TO 30 3315 IF (S(L) .EQ. 0.0D0) GO TO 30 3316C 3317C APPLY THE TRANSFORMATION. 3318C 3319* 3320* INCREMENT OP COUNT 3321 IOPS = IOPS + (4*(N-L)+5) 3322 T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) 3323 CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1) 3324 30 CONTINUE 3325C 3326C PLACE THE L-TH ROW OF X INTO E FOR THE 3327C SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION. 3328C 3329 E(J) = X(L,J) 3330 40 CONTINUE 3331 50 CONTINUE 3332 IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70 3333C 3334C PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK 3335C MULTIPLICATION. 3336C 3337 DO 60 I = L, N 3338 U(I,L) = X(I,L) 3339 60 CONTINUE 3340 70 CONTINUE 3341 IF (L .GT. NRT) GO TO 150 3342C 3343C COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE 3344C L-TH SUPER-DIAGONAL IN E(L). 3345C 3346* 3347* INCREMENT OP COUNT 3348 IOPS = IOPS + (2*(P-L)+1) 3349 E(L) = DNRM2(P-L,E(LP1),1) 3350 IF (E(L) .EQ. 0.0D0) GO TO 80 3351 IF (E(LP1) .NE. 0.0D0) E(L) = DSIGN(E(L),E(LP1)) 3352* 3353* INCREMENT OP COUNT 3354 IOPS = IOPS + (P-L+2) 3355 CALL DSCAL(P-L,1.0D0/E(L),E(LP1),1) 3356 E(LP1) = 1.0D0 + E(LP1) 3357 80 CONTINUE 3358 E(L) = -E(L) 3359 IF (LP1 .GT. N .OR. E(L) .EQ. 0.0D0) GO TO 120 3360C 3361C APPLY THE TRANSFORMATION. 3362C 3363 DO 90 I = LP1, N 3364 WORK(I) = 0.0D0 3365 90 CONTINUE 3366* 3367* INCREMENT OP COUNT 3368 IOPS = IOPS + DBLE(4*(N-L)+1)*(P-L) 3369 DO 100 J = LP1, P 3370 CALL DAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1) 3371 100 CONTINUE 3372 DO 110 J = LP1, P 3373 CALL DAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1) 3374 110 CONTINUE 3375 120 CONTINUE 3376 IF (.NOT.WANTV) GO TO 140 3377C 3378C PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT 3379C BACK MULTIPLICATION. 3380C 3381 DO 130 I = LP1, P 3382 V(I,L) = E(I) 3383 130 CONTINUE 3384 140 CONTINUE 3385 150 CONTINUE 3386 160 CONTINUE 3387 170 CONTINUE 3388C 3389C SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M. 3390C 3391 M = MIN0(P,N+1) 3392 NCTP1 = NCT + 1 3393 NRTP1 = NRT + 1 3394 IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1) 3395 IF (N .LT. M) S(M) = 0.0D0 3396 IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M) 3397 E(M) = 0.0D0 3398C 3399C IF REQUIRED, GENERATE U. 3400C 3401 IF (.NOT.WANTU) GO TO 300 3402 IF (NCU .LT. NCTP1) GO TO 200 3403 DO 190 J = NCTP1, NCU 3404 DO 180 I = 1, N 3405 U(I,J) = 0.0D0 3406 180 CONTINUE 3407 U(J,J) = 1.0D0 3408 190 CONTINUE 3409 200 CONTINUE 3410 IF (NCT .LT. 1) GO TO 290 3411 DO 280 LL = 1, NCT 3412 L = NCT - LL + 1 3413 IF (S(L) .EQ. 0.0D0) GO TO 250 3414 LP1 = L + 1 3415 IF (NCU .LT. LP1) GO TO 220 3416* 3417* INCREMENT OP COUNT 3418 IOPS = IOPS + (DBLE(4*(N-L)+5)*(NCU-L)+(N-L+2)) 3419 DO 210 J = LP1, NCU 3420 T = -DDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L) 3421 CALL DAXPY(N-L+1,T,U(L,L),1,U(L,J),1) 3422 210 CONTINUE 3423 220 CONTINUE 3424 CALL DSCAL(N-L+1,-1.0D0,U(L,L),1) 3425 U(L,L) = 1.0D0 + U(L,L) 3426 LM1 = L - 1 3427 IF (LM1 .LT. 1) GO TO 240 3428 DO 230 I = 1, LM1 3429 U(I,L) = 0.0D0 3430 230 CONTINUE 3431 240 CONTINUE 3432 GO TO 270 3433 250 CONTINUE 3434 DO 260 I = 1, N 3435 U(I,L) = 0.0D0 3436 260 CONTINUE 3437 U(L,L) = 1.0D0 3438 270 CONTINUE 3439 280 CONTINUE 3440 290 CONTINUE 3441 300 CONTINUE 3442C 3443C IF IT IS REQUIRED, GENERATE V. 3444C 3445 IF (.NOT.WANTV) GO TO 350 3446 DO 340 LL = 1, P 3447 L = P - LL + 1 3448 LP1 = L + 1 3449 IF (L .GT. NRT) GO TO 320 3450 IF (E(L) .EQ. 0.0D0) GO TO 320 3451* 3452* INCREMENT OP COUNT 3453 IOPS = IOPS + DBLE(4*(P-L)+1)*(P-L) 3454 DO 310 J = LP1, P 3455 T = -DDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L) 3456 CALL DAXPY(P-L,T,V(LP1,L),1,V(LP1,J),1) 3457 310 CONTINUE 3458 320 CONTINUE 3459 DO 330 I = 1, P 3460 V(I,L) = 0.0D0 3461 330 CONTINUE 3462 V(L,L) = 1.0D0 3463 340 CONTINUE 3464 350 CONTINUE 3465C 3466C MAIN ITERATION LOOP FOR THE SINGULAR VALUES. 3467C 3468 MM = M 3469* 3470* INITIALIZE ITERATION COUNTER 3471 ITCNT = 0 3472 ITER = 0 3473 360 CONTINUE 3474C 3475C QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND. 3476C 3477C ...EXIT 3478 IF (M .EQ. 0) GO TO 620 3479C 3480C IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET 3481C FLAG AND RETURN. 3482C 3483* 3484* UPDATE ITERATION COUNTER 3485 ITCNT = ITER 3486 IF (ITER .LT. MAXIT) GO TO 370 3487 INFO = M 3488C ......EXIT 3489 GO TO 620 3490 370 CONTINUE 3491C 3492C THIS SECTION OF THE PROGRAM INSPECTS FOR 3493C NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS. ON 3494C COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS. 3495C 3496C KASE = 1 IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M 3497C KASE = 2 IF S(L) IS NEGLIGIBLE AND L.LT.M 3498C KASE = 3 IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND 3499C S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP). 3500C KASE = 4 IF E(M-1) IS NEGLIGIBLE (CONVERGENCE). 3501C 3502 DO 390 LL = 1, M 3503 L = M - LL 3504C ...EXIT 3505 IF (L .EQ. 0) GO TO 400 3506* 3507* INCREMENT OP COUNT 3508 IOPST = IOPST + 2 3509 TEST = DABS(S(L)) + DABS(S(L+1)) 3510* 3511* REPLACE STOPPING CRITERION WITH NEW ONE AS IN LAPACK 3512* 3513* ZTEST = TEST + DABS(E(L)) 3514* IF (ZTEST .NE. TEST) GO TO 380 3515 IF (DABS(E(L)) .GT. EPS * TEST) GOTO 380 3516* 3517 E(L) = 0.0D0 3518C ......EXIT 3519 GO TO 400 3520 380 CONTINUE 3521 390 CONTINUE 3522 400 CONTINUE 3523 IF (L .NE. M - 1) GO TO 410 3524 KASE = 4 3525 GO TO 480 3526 410 CONTINUE 3527 LP1 = L + 1 3528 MP1 = M + 1 3529 DO 430 LLS = LP1, MP1 3530 LS = M - LLS + LP1 3531C ...EXIT 3532 IF (LS .EQ. L) GO TO 440 3533 TEST = 0.0D0 3534* 3535* INCREMENT OP COUNT 3536 IOPST = IOPST + 3 3537 IF (LS .NE. M) TEST = TEST + DABS(E(LS)) 3538 IF (LS .NE. L + 1) TEST = TEST + DABS(E(LS-1)) 3539* 3540* REPLACE STOPPING CRITERION WITH NEW ONE AS IN LAPACK 3541* 3542* ZTEST = TEST + DABS(S(LS)) 3543* IF (ZTEST .NE. TEST) GO TO 420 3544 IF (DABS(S(LS)) .GT. EPS * TEST) GOTO 420 3545* 3546 S(LS) = 0.0D0 3547C ......EXIT 3548 GO TO 440 3549 420 CONTINUE 3550 430 CONTINUE 3551 440 CONTINUE 3552 IF (LS .NE. L) GO TO 450 3553 KASE = 3 3554 GO TO 470 3555 450 CONTINUE 3556 IF (LS .NE. M) GO TO 460 3557 KASE = 1 3558 GO TO 470 3559 460 CONTINUE 3560 KASE = 2 3561 L = LS 3562 470 CONTINUE 3563 480 CONTINUE 3564 L = L + 1 3565C 3566C PERFORM THE TASK INDICATED BY KASE. 3567C 3568 GO TO (490,520,540,570), KASE 3569C 3570C DEFLATE NEGLIGIBLE S(M). 3571C 3572 490 CONTINUE 3573 MM1 = M - 1 3574 F = E(M-1) 3575 E(M-1) = 0.0D0 3576* 3577* INCREMENT OP COUNT 3578 IOPS = IOPS + ((MM1-L+1)*13 - 2) 3579 IF (WANTV) IOPS = IOPS + DBLE(MM1-L+1)*6*P 3580 DO 510 KK = L, MM1 3581 K = MM1 - KK + L 3582 T1 = S(K) 3583 CALL DROTG(T1,F,CS,SN) 3584 S(K) = T1 3585 IF (K .EQ. L) GO TO 500 3586 F = -SN*E(K-1) 3587 E(K-1) = CS*E(K-1) 3588 500 CONTINUE 3589 IF (WANTV) CALL DROT(P,V(1,K),1,V(1,M),1,CS,SN) 3590 510 CONTINUE 3591 GO TO 610 3592C 3593C SPLIT AT NEGLIGIBLE S(L). 3594C 3595 520 CONTINUE 3596 F = E(L-1) 3597 E(L-1) = 0.0D0 3598* 3599* INCREMENT OP COUNT 3600 IOPS = IOPS + (M-L+1)*13 3601 IF (WANTU) IOPS = IOPS + DBLE(M-L+1)*6*N 3602 DO 530 K = L, M 3603 T1 = S(K) 3604 CALL DROTG(T1,F,CS,SN) 3605 S(K) = T1 3606 F = -SN*E(K) 3607 E(K) = CS*E(K) 3608 IF (WANTU) CALL DROT(N,U(1,K),1,U(1,L-1),1,CS,SN) 3609 530 CONTINUE 3610 GO TO 610 3611C 3612C PERFORM ONE QR STEP. 3613C 3614 540 CONTINUE 3615C 3616C CALCULATE THE SHIFT. 3617C 3618* 3619* INCREMENT OP COUNT 3620 IOPST = IOPST + 23 3621 SCALE = DMAX1(DABS(S(M)),DABS(S(M-1)),DABS(E(M-1)), 3622 * DABS(S(L)),DABS(E(L))) 3623 SM = S(M)/SCALE 3624 SMM1 = S(M-1)/SCALE 3625 EMM1 = E(M-1)/SCALE 3626 SL = S(L)/SCALE 3627 EL = E(L)/SCALE 3628 B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0 3629 C = (SM*EMM1)**2 3630 SHIFT = 0.0D0 3631 IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 550 3632 SHIFT = DSQRT(B**2+C) 3633 IF (B .LT. 0.0D0) SHIFT = -SHIFT 3634 SHIFT = C/(B + SHIFT) 3635 550 CONTINUE 3636 F = (SL + SM)*(SL - SM) + SHIFT 3637 G = SL*EL 3638C 3639C CHASE ZEROS. 3640C 3641 MM1 = M - 1 3642* 3643* INCREMENT OP COUNT 3644 IOPS = IOPS + (MM1-L+1)*38 3645 IF (WANTV) IOPS = IOPS+DBLE(MM1-L+1)*6*P 3646 IF (WANTU) IOPS = IOPS+DBLE(MAX((MIN(MM1,N-1)-L+1),0))*6*N 3647 DO 560 K = L, MM1 3648 CALL DROTG(F,G,CS,SN) 3649 IF (K .NE. L) E(K-1) = F 3650 F = CS*S(K) + SN*E(K) 3651 E(K) = CS*E(K) - SN*S(K) 3652 G = SN*S(K+1) 3653 S(K+1) = CS*S(K+1) 3654 IF (WANTV) CALL DROT(P,V(1,K),1,V(1,K+1),1,CS,SN) 3655 CALL DROTG(F,G,CS,SN) 3656 S(K) = F 3657 F = CS*E(K) + SN*S(K+1) 3658 S(K+1) = -SN*E(K) + CS*S(K+1) 3659 G = SN*E(K+1) 3660 E(K+1) = CS*E(K+1) 3661 IF (WANTU .AND. K .LT. N) 3662 * CALL DROT(N,U(1,K),1,U(1,K+1),1,CS,SN) 3663 560 CONTINUE 3664 E(M-1) = F 3665 ITER = ITER + 1 3666 GO TO 610 3667C 3668C CONVERGENCE. 3669C 3670 570 CONTINUE 3671C 3672C MAKE THE SINGULAR VALUE POSITIVE. 3673C 3674 IF (S(L) .GE. 0.0D0) GO TO 580 3675 S(L) = -S(L) 3676* 3677* INCREMENT OP COUNT 3678 IF (WANTV) IOPS = IOPS + P 3679 IF (WANTV) CALL DSCAL(P,-1.0D0,V(1,L),1) 3680 580 CONTINUE 3681C 3682C ORDER THE SINGULAR VALUE. 3683C 3684 590 IF (L .EQ. MM) GO TO 600 3685C ...EXIT 3686 IF (S(L) .GE. S(L+1)) GO TO 600 3687 T = S(L) 3688 S(L) = S(L+1) 3689 S(L+1) = T 3690 IF (WANTV .AND. L .LT. P) 3691 * CALL DSWAP(P,V(1,L),1,V(1,L+1),1) 3692 IF (WANTU .AND. L .LT. N) 3693 * CALL DSWAP(N,U(1,L),1,U(1,L+1),1) 3694 L = L + 1 3695 GO TO 590 3696 600 CONTINUE 3697 ITER = 0 3698 M = M - 1 3699 610 CONTINUE 3700 GO TO 360 3701 620 CONTINUE 3702* 3703* COMPUTE FINAL OPCOUNT 3704 IOPS = IOPS + IOPST 3705 RETURN 3706 END 3707 SUBROUTINE QZHES(NM,N,A,B,MATZ,Z) 3708C 3709 INTEGER I,J,K,L,N,LB,L1,NM,NK1,NM1,NM2 3710 DOUBLE PRECISION A(NM,N),B(NM,N),Z(NM,N) 3711 DOUBLE PRECISION R,S,T,U1,U2,V1,V2,RHO 3712 LOGICAL MATZ 3713* 3714* ---------------------- BEGIN TIMING CODE ------------------------- 3715* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 3716* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 3717* OPST IS USED TO ACCUMULATE SMALL CONTRIBUTIONS TO OPS 3718* TO AVOID ROUNDOFF ERROR 3719* .. COMMON BLOCKS .. 3720 COMMON / LATIME / OPS, ITCNT 3721* .. 3722* .. SCALARS IN COMMON .. 3723 DOUBLE PRECISION ITCNT, OPS 3724* .. 3725* ----------------------- END TIMING CODE -------------------------- 3726* 3727C 3728C THIS SUBROUTINE IS THE FIRST STEP OF THE QZ ALGORITHM 3729C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, 3730C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. 3731C 3732C THIS SUBROUTINE ACCEPTS A PAIR OF REAL GENERAL MATRICES AND 3733C REDUCES ONE OF THEM TO UPPER HESSENBERG FORM AND THE OTHER 3734C TO UPPER TRIANGULAR FORM USING ORTHOGONAL TRANSFORMATIONS. 3735C IT IS USUALLY FOLLOWED BY QZIT, QZVAL AND, POSSIBLY, QZVEC. 3736C 3737C ON INPUT 3738C 3739C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 3740C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 3741C DIMENSION STATEMENT. 3742C 3743C N IS THE ORDER OF THE MATRICES. 3744C 3745C A CONTAINS A REAL GENERAL MATRIX. 3746C 3747C B CONTAINS A REAL GENERAL MATRIX. 3748C 3749C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS 3750C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING 3751C EIGENVECTORS, AND TO .FALSE. OTHERWISE. 3752C 3753C ON OUTPUT 3754C 3755C A HAS BEEN REDUCED TO UPPER HESSENBERG FORM. THE ELEMENTS 3756C BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO. 3757C 3758C B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS 3759C BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO. 3760C 3761C Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS IF 3762C MATZ HAS BEEN SET TO .TRUE. OTHERWISE, Z IS NOT REFERENCED. 3763C 3764C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 3765C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 3766C 3767C THIS VERSION DATED AUGUST 1983. 3768C 3769C ------------------------------------------------------------------ 3770C 3771C .......... INITIALIZE Z .......... 3772 IF (.NOT. MATZ) GO TO 10 3773C 3774 DO 3 J = 1, N 3775C 3776 DO 2 I = 1, N 3777 Z(I,J) = 0.0D0 3778 2 CONTINUE 3779C 3780 Z(J,J) = 1.0D0 3781 3 CONTINUE 3782C .......... REDUCE B TO UPPER TRIANGULAR FORM .......... 3783 10 IF (N .LE. 1) GO TO 170 3784 NM1 = N - 1 3785C 3786 DO 100 L = 1, NM1 3787 L1 = L + 1 3788 S = 0.0D0 3789C 3790 DO 20 I = L1, N 3791 S = S + DABS(B(I,L)) 3792 20 CONTINUE 3793C 3794 IF (S .EQ. 0.0D0) GO TO 100 3795 S = S + DABS(B(L,L)) 3796 R = 0.0D0 3797C 3798 DO 25 I = L, N 3799 B(I,L) = B(I,L) / S 3800 R = R + B(I,L)**2 3801 25 CONTINUE 3802C 3803 R = DSIGN(DSQRT(R),B(L,L)) 3804 B(L,L) = B(L,L) + R 3805 RHO = R * B(L,L) 3806C 3807 DO 50 J = L1, N 3808 T = 0.0D0 3809C 3810 DO 30 I = L, N 3811 T = T + B(I,L) * B(I,J) 3812 30 CONTINUE 3813C 3814 T = -T / RHO 3815C 3816 DO 40 I = L, N 3817 B(I,J) = B(I,J) + T * B(I,L) 3818 40 CONTINUE 3819C 3820 50 CONTINUE 3821C 3822 DO 80 J = 1, N 3823 T = 0.0D0 3824C 3825 DO 60 I = L, N 3826 T = T + B(I,L) * A(I,J) 3827 60 CONTINUE 3828C 3829 T = -T / RHO 3830C 3831 DO 70 I = L, N 3832 A(I,J) = A(I,J) + T * B(I,L) 3833 70 CONTINUE 3834C 3835 80 CONTINUE 3836C 3837 B(L,L) = -S * R 3838C 3839 DO 90 I = L1, N 3840 B(I,L) = 0.0D0 3841 90 CONTINUE 3842C 3843 100 CONTINUE 3844* 3845* ---------------------- BEGIN TIMING CODE ------------------------- 3846 OPS = OPS + DBLE( 8*N**2 + 17*N + 24 )*DBLE( N-1 ) / 3.0D0 3847* ----------------------- END TIMING CODE -------------------------- 3848* 3849C .......... REDUCE A TO UPPER HESSENBERG FORM, WHILE 3850C KEEPING B TRIANGULAR .......... 3851 IF (N .EQ. 2) GO TO 170 3852 NM2 = N - 2 3853C 3854 DO 160 K = 1, NM2 3855 NK1 = NM1 - K 3856C .......... FOR L=N-1 STEP -1 UNTIL K+1 DO -- .......... 3857 DO 150 LB = 1, NK1 3858 L = N - LB 3859 L1 = L + 1 3860C .......... ZERO A(L+1,K) .......... 3861 S = DABS(A(L,K)) + DABS(A(L1,K)) 3862 IF (S .EQ. 0.0D0) GO TO 150 3863 U1 = A(L,K) / S 3864 U2 = A(L1,K) / S 3865 R = DSIGN(DSQRT(U1*U1+U2*U2),U1) 3866 V1 = -(U1 + R) / R 3867 V2 = -U2 / R 3868 U2 = V2 / V1 3869C 3870 DO 110 J = K, N 3871 T = A(L,J) + U2 * A(L1,J) 3872 A(L,J) = A(L,J) + T * V1 3873 A(L1,J) = A(L1,J) + T * V2 3874 110 CONTINUE 3875C 3876 A(L1,K) = 0.0D0 3877C 3878 DO 120 J = L, N 3879 T = B(L,J) + U2 * B(L1,J) 3880 B(L,J) = B(L,J) + T * V1 3881 B(L1,J) = B(L1,J) + T * V2 3882 120 CONTINUE 3883C .......... ZERO B(L+1,L) .......... 3884 S = DABS(B(L1,L1)) + DABS(B(L1,L)) 3885 IF (S .EQ. 0.0D0) GO TO 150 3886 U1 = B(L1,L1) / S 3887 U2 = B(L1,L) / S 3888 R = DSIGN(DSQRT(U1*U1+U2*U2),U1) 3889 V1 = -(U1 + R) / R 3890 V2 = -U2 / R 3891 U2 = V2 / V1 3892C 3893 DO 130 I = 1, L1 3894 T = B(I,L1) + U2 * B(I,L) 3895 B(I,L1) = B(I,L1) + T * V1 3896 B(I,L) = B(I,L) + T * V2 3897 130 CONTINUE 3898C 3899 B(L1,L) = 0.0D0 3900C 3901 DO 140 I = 1, N 3902 T = A(I,L1) + U2 * A(I,L) 3903 A(I,L1) = A(I,L1) + T * V1 3904 A(I,L) = A(I,L) + T * V2 3905 140 CONTINUE 3906C 3907 IF (.NOT. MATZ) GO TO 150 3908C 3909 DO 145 I = 1, N 3910 T = Z(I,L1) + U2 * Z(I,L) 3911 Z(I,L1) = Z(I,L1) + T * V1 3912 Z(I,L) = Z(I,L) + T * V2 3913 145 CONTINUE 3914C 3915 150 CONTINUE 3916C 3917 160 CONTINUE 3918C 3919* 3920* ---------------------- BEGIN TIMING CODE ------------------------- 3921 IF( MATZ ) THEN 3922 OPS = OPS + DBLE( 11*N + 20 )*DBLE( N-1 )*DBLE( N-2 ) 3923 ELSE 3924 OPS = OPS + DBLE( 8*N + 20 )*DBLE( N-1 )*DBLE( N-2 ) 3925 END IF 3926* ----------------------- END TIMING CODE -------------------------- 3927* 3928 170 RETURN 3929 END 3930 SUBROUTINE QZIT(NM,N,A,B,EPS1,MATZ,Z,IERR) 3931C 3932 INTEGER I,J,K,L,N,EN,K1,K2,LD,LL,L1,NA,NM,ISH,ITN,ITS,KM1,LM1, 3933 X ENM2,IERR,LOR1,ENORN 3934 DOUBLE PRECISION A(NM,N),B(NM,N),Z(NM,N) 3935 DOUBLE PRECISION R,S,T,A1,A2,A3,EP,SH,U1,U2,U3,V1,V2,V3,ANI,A11, 3936 X A12,A21,A22,A33,A34,A43,A44,BNI,B11,B12,B22,B33,B34, 3937 X B44,EPSA,EPSB,EPS1,ANORM,BNORM,EPSLON 3938 LOGICAL MATZ,NOTLAS 3939* 3940* ---------------------- BEGIN TIMING CODE ------------------------- 3941* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 3942* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 3943* OPST IS USED TO ACCUMULATE SMALL CONTRIBUTIONS TO OPS 3944* TO AVOID ROUNDOFF ERROR 3945* .. COMMON BLOCKS .. 3946 COMMON / LATIME / OPS, ITCNT 3947* .. 3948* .. SCALARS IN COMMON .. 3949 DOUBLE PRECISION ITCNT, OPS 3950* .. 3951 DOUBLE PRECISION OPST 3952* ----------------------- END TIMING CODE -------------------------- 3953* 3954C 3955C THIS SUBROUTINE IS THE SECOND STEP OF THE QZ ALGORITHM 3956C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, 3957C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART, 3958C AS MODIFIED IN TECHNICAL NOTE NASA TN D-7305(1973) BY WARD. 3959C 3960C THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM 3961C IN UPPER HESSENBERG FORM AND THE OTHER IN UPPER TRIANGULAR FORM. 3962C IT REDUCES THE HESSENBERG MATRIX TO QUASI-TRIANGULAR FORM USING 3963C ORTHOGONAL TRANSFORMATIONS WHILE MAINTAINING THE TRIANGULAR FORM 3964C OF THE OTHER MATRIX. IT IS USUALLY PRECEDED BY QZHES AND 3965C FOLLOWED BY QZVAL AND, POSSIBLY, QZVEC. 3966C 3967C ON INPUT 3968C 3969C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 3970C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 3971C DIMENSION STATEMENT. 3972C 3973C N IS THE ORDER OF THE MATRICES. 3974C 3975C A CONTAINS A REAL UPPER HESSENBERG MATRIX. 3976C 3977C B CONTAINS A REAL UPPER TRIANGULAR MATRIX. 3978C 3979C EPS1 IS A TOLERANCE USED TO DETERMINE NEGLIGIBLE ELEMENTS. 3980C EPS1 = 0.0 (OR NEGATIVE) MAY BE INPUT, IN WHICH CASE AN 3981C ELEMENT WILL BE NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF 3982C ERROR TIMES THE NORM OF ITS MATRIX. IF THE INPUT EPS1 IS 3983C POSITIVE, THEN AN ELEMENT WILL BE CONSIDERED NEGLIGIBLE 3984C IF IT IS LESS THAN EPS1 TIMES THE NORM OF ITS MATRIX. A 3985C POSITIVE VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION, 3986C BUT LESS ACCURATE RESULTS. 3987C 3988C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS 3989C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING 3990C EIGENVECTORS, AND TO .FALSE. OTHERWISE. 3991C 3992C Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE 3993C TRANSFORMATION MATRIX PRODUCED IN THE REDUCTION 3994C BY QZHES, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. 3995C IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. 3996C 3997C ON OUTPUT 3998C 3999C A HAS BEEN REDUCED TO QUASI-TRIANGULAR FORM. THE ELEMENTS 4000C BELOW THE FIRST SUBDIAGONAL ARE STILL ZERO AND NO TWO 4001C CONSECUTIVE SUBDIAGONAL ELEMENTS ARE NONZERO. 4002C 4003C B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS 4004C HAVE BEEN ALTERED. THE LOCATION B(N,1) IS USED TO STORE 4005C EPS1 TIMES THE NORM OF B FOR LATER USE BY QZVAL AND QZVEC. 4006C 4007C Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS 4008C (FOR BOTH STEPS) IF MATZ HAS BEEN SET TO .TRUE.. 4009C 4010C IERR IS SET TO 4011C ZERO FOR NORMAL RETURN, 4012C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED 4013C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. 4014C 4015C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 4016C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 4017C 4018C THIS VERSION DATED AUGUST 1983. 4019C 4020C ------------------------------------------------------------------ 4021C 4022 IERR = 0 4023C .......... COMPUTE EPSA,EPSB .......... 4024 ANORM = 0.0D0 4025 BNORM = 0.0D0 4026C 4027 DO 30 I = 1, N 4028 ANI = 0.0D0 4029 IF (I .NE. 1) ANI = DABS(A(I,I-1)) 4030 BNI = 0.0D0 4031C 4032 DO 20 J = I, N 4033 ANI = ANI + DABS(A(I,J)) 4034 BNI = BNI + DABS(B(I,J)) 4035 20 CONTINUE 4036C 4037 IF (ANI .GT. ANORM) ANORM = ANI 4038 IF (BNI .GT. BNORM) BNORM = BNI 4039 30 CONTINUE 4040* 4041* ---------------------- BEGIN TIMING CODE ------------------------- 4042 OPS = OPS + DBLE( N*( N+1 ) ) 4043 OPST = 0.0D0 4044 ITCNT = 0 4045* ----------------------- END TIMING CODE -------------------------- 4046* 4047C 4048 IF (ANORM .EQ. 0.0D0) ANORM = 1.0D0 4049 IF (BNORM .EQ. 0.0D0) BNORM = 1.0D0 4050 EP = EPS1 4051 IF (EP .GT. 0.0D0) GO TO 50 4052C .......... USE ROUNDOFF LEVEL IF EPS1 IS ZERO .......... 4053 EP = EPSLON(1.0D0) 4054 50 EPSA = EP * ANORM 4055 EPSB = EP * BNORM 4056C .......... REDUCE A TO QUASI-TRIANGULAR FORM, WHILE 4057C KEEPING B TRIANGULAR .......... 4058 LOR1 = 1 4059 ENORN = N 4060 EN = N 4061 ITN = 30*N 4062C .......... BEGIN QZ STEP .......... 4063 60 IF (EN .LE. 2) GO TO 1001 4064 IF (.NOT. MATZ) ENORN = EN 4065 ITS = 0 4066 NA = EN - 1 4067 ENM2 = NA - 1 4068 70 ISH = 2 4069* 4070* ---------------------- BEGIN TIMING CODE ------------------------- 4071 OPS = OPS + OPST 4072 OPST = 0.0D0 4073 ITCNT = ITCNT + 1 4074* ----------------------- END TIMING CODE -------------------------- 4075* 4076C .......... CHECK FOR CONVERGENCE OR REDUCIBILITY. 4077C FOR L=EN STEP -1 UNTIL 1 DO -- .......... 4078 DO 80 LL = 1, EN 4079 LM1 = EN - LL 4080 L = LM1 + 1 4081 IF (L .EQ. 1) GO TO 95 4082 IF (DABS(A(L,LM1)) .LE. EPSA) GO TO 90 4083 80 CONTINUE 4084C 4085 90 A(L,LM1) = 0.0D0 4086 IF (L .LT. NA) GO TO 95 4087C .......... 1-BY-1 OR 2-BY-2 BLOCK ISOLATED .......... 4088 EN = LM1 4089 GO TO 60 4090C .......... CHECK FOR SMALL TOP OF B .......... 4091 95 LD = L 4092 100 L1 = L + 1 4093 B11 = B(L,L) 4094 IF (DABS(B11) .GT. EPSB) GO TO 120 4095 B(L,L) = 0.0D0 4096 S = DABS(A(L,L)) + DABS(A(L1,L)) 4097 U1 = A(L,L) / S 4098 U2 = A(L1,L) / S 4099 R = DSIGN(DSQRT(U1*U1+U2*U2),U1) 4100 V1 = -(U1 + R) / R 4101 V2 = -U2 / R 4102 U2 = V2 / V1 4103C 4104 DO 110 J = L, ENORN 4105 T = A(L,J) + U2 * A(L1,J) 4106 A(L,J) = A(L,J) + T * V1 4107 A(L1,J) = A(L1,J) + T * V2 4108 T = B(L,J) + U2 * B(L1,J) 4109 B(L,J) = B(L,J) + T * V1 4110 B(L1,J) = B(L1,J) + T * V2 4111 110 CONTINUE 4112C 4113* ---------------------- BEGIN TIMING CODE ------------------------- 4114 OPST = OPST + DBLE( 12*( ENORN+1-L ) + 11 ) 4115* ----------------------- END TIMING CODE -------------------------- 4116 IF (L .NE. 1) A(L,LM1) = -A(L,LM1) 4117 LM1 = L 4118 L = L1 4119 GO TO 90 4120 120 A11 = A(L,L) / B11 4121 A21 = A(L1,L) / B11 4122 IF (ISH .EQ. 1) GO TO 140 4123C .......... ITERATION STRATEGY .......... 4124 IF (ITN .EQ. 0) GO TO 1000 4125 IF (ITS .EQ. 10) GO TO 155 4126C .......... DETERMINE TYPE OF SHIFT .......... 4127 B22 = B(L1,L1) 4128 IF (DABS(B22) .LT. EPSB) B22 = EPSB 4129 B33 = B(NA,NA) 4130 IF (DABS(B33) .LT. EPSB) B33 = EPSB 4131 B44 = B(EN,EN) 4132 IF (DABS(B44) .LT. EPSB) B44 = EPSB 4133 A33 = A(NA,NA) / B33 4134 A34 = A(NA,EN) / B44 4135 A43 = A(EN,NA) / B33 4136 A44 = A(EN,EN) / B44 4137 B34 = B(NA,EN) / B44 4138 T = 0.5D0 * (A43 * B34 - A33 - A44) 4139 R = T * T + A34 * A43 - A33 * A44 4140* ---------------------- BEGIN TIMING CODE ------------------------- 4141 OPST = OPST + DBLE( 16 ) 4142* ----------------------- END TIMING CODE -------------------------- 4143 IF (R .LT. 0.0D0) GO TO 150 4144C .......... DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A .......... 4145 ISH = 1 4146 R = DSQRT(R) 4147 SH = -T + R 4148 S = -T - R 4149 IF (DABS(S-A44) .LT. DABS(SH-A44)) SH = S 4150C .......... LOOK FOR TWO CONSECUTIVE SMALL 4151C SUB-DIAGONAL ELEMENTS OF A. 4152C FOR L=EN-2 STEP -1 UNTIL LD DO -- .......... 4153 DO 130 LL = LD, ENM2 4154 L = ENM2 + LD - LL 4155 IF (L .EQ. LD) GO TO 140 4156 LM1 = L - 1 4157 L1 = L + 1 4158 T = A(L,L) 4159 IF (DABS(B(L,L)) .GT. EPSB) T = T - SH * B(L,L) 4160* --------------------- BEGIN TIMING CODE ----------------------- 4161 IF (DABS(A(L,LM1)) .LE. DABS(T/A(L1,L)) * EPSA) THEN 4162 OPST = OPST + DBLE( 5 + 4*( LL+1-LD ) ) 4163 GO TO 100 4164 END IF 4165* ---------------------- END TIMING CODE ------------------------ 4166 130 CONTINUE 4167* ---------------------- BEGIN TIMING CODE ------------------------- 4168 OPST = OPST + DBLE( 5 + 4*( ENM2+1-LD ) ) 4169* ----------------------- END TIMING CODE -------------------------- 4170C 4171 140 A1 = A11 - SH 4172 A2 = A21 4173 IF (L .NE. LD) A(L,LM1) = -A(L,LM1) 4174 GO TO 160 4175C .......... DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A .......... 4176 150 A12 = A(L,L1) / B22 4177 A22 = A(L1,L1) / B22 4178 B12 = B(L,L1) / B22 4179 A1 = ((A33 - A11) * (A44 - A11) - A34 * A43 + A43 * B34 * A11) 4180 X / A21 + A12 - A11 * B12 4181 A2 = (A22 - A11) - A21 * B12 - (A33 - A11) - (A44 - A11) 4182 X + A43 * B34 4183 A3 = A(L1+1,L1) / B22 4184* ---------------------- BEGIN TIMING CODE ------------------------- 4185 OPST = OPST + DBLE( 25 ) 4186* ----------------------- END TIMING CODE -------------------------- 4187 GO TO 160 4188C .......... AD HOC SHIFT .......... 4189 155 A1 = 0.0D0 4190 A2 = 1.0D0 4191 A3 = 1.1605D0 4192 160 ITS = ITS + 1 4193 ITN = ITN - 1 4194 IF (.NOT. MATZ) LOR1 = LD 4195C .......... MAIN LOOP .......... 4196 DO 260 K = L, NA 4197 NOTLAS = K .NE. NA .AND. ISH .EQ. 2 4198 K1 = K + 1 4199 K2 = K + 2 4200 KM1 = MAX0(K-1,L) 4201 LL = MIN0(EN,K1+ISH) 4202 IF (NOTLAS) GO TO 190 4203C .......... ZERO A(K+1,K-1) .......... 4204 IF (K .EQ. L) GO TO 170 4205 A1 = A(K,KM1) 4206 A2 = A(K1,KM1) 4207 170 S = DABS(A1) + DABS(A2) 4208 IF (S .EQ. 0.0D0) GO TO 70 4209 U1 = A1 / S 4210 U2 = A2 / S 4211 R = DSIGN(DSQRT(U1*U1+U2*U2),U1) 4212 V1 = -(U1 + R) / R 4213 V2 = -U2 / R 4214 U2 = V2 / V1 4215C 4216 DO 180 J = KM1, ENORN 4217 T = A(K,J) + U2 * A(K1,J) 4218 A(K,J) = A(K,J) + T * V1 4219 A(K1,J) = A(K1,J) + T * V2 4220 T = B(K,J) + U2 * B(K1,J) 4221 B(K,J) = B(K,J) + T * V1 4222 B(K1,J) = B(K1,J) + T * V2 4223 180 CONTINUE 4224C 4225* --------------------- BEGIN TIMING CODE ----------------------- 4226 OPST = OPST + DBLE( 11 + 12*( ENORN+1-KM1 ) ) 4227* ---------------------- END TIMING CODE ------------------------ 4228 IF (K .NE. L) A(K1,KM1) = 0.0D0 4229 GO TO 240 4230C .......... ZERO A(K+1,K-1) AND A(K+2,K-1) .......... 4231 190 IF (K .EQ. L) GO TO 200 4232 A1 = A(K,KM1) 4233 A2 = A(K1,KM1) 4234 A3 = A(K2,KM1) 4235 200 S = DABS(A1) + DABS(A2) + DABS(A3) 4236 IF (S .EQ. 0.0D0) GO TO 260 4237 U1 = A1 / S 4238 U2 = A2 / S 4239 U3 = A3 / S 4240 R = DSIGN(DSQRT(U1*U1+U2*U2+U3*U3),U1) 4241 V1 = -(U1 + R) / R 4242 V2 = -U2 / R 4243 V3 = -U3 / R 4244 U2 = V2 / V1 4245 U3 = V3 / V1 4246C 4247 DO 210 J = KM1, ENORN 4248 T = A(K,J) + U2 * A(K1,J) + U3 * A(K2,J) 4249 A(K,J) = A(K,J) + T * V1 4250 A(K1,J) = A(K1,J) + T * V2 4251 A(K2,J) = A(K2,J) + T * V3 4252 T = B(K,J) + U2 * B(K1,J) + U3 * B(K2,J) 4253 B(K,J) = B(K,J) + T * V1 4254 B(K1,J) = B(K1,J) + T * V2 4255 B(K2,J) = B(K2,J) + T * V3 4256 210 CONTINUE 4257* --------------------- BEGIN TIMING CODE ----------------------- 4258 OPST = OPST + DBLE( 17 + 20*( ENORN+1-KM1 ) ) 4259* ---------------------- END TIMING CODE ------------------------ 4260C 4261 IF (K .EQ. L) GO TO 220 4262 A(K1,KM1) = 0.0D0 4263 A(K2,KM1) = 0.0D0 4264C .......... ZERO B(K+2,K+1) AND B(K+2,K) .......... 4265 220 S = DABS(B(K2,K2)) + DABS(B(K2,K1)) + DABS(B(K2,K)) 4266 IF (S .EQ. 0.0D0) GO TO 240 4267 U1 = B(K2,K2) / S 4268 U2 = B(K2,K1) / S 4269 U3 = B(K2,K) / S 4270 R = DSIGN(DSQRT(U1*U1+U2*U2+U3*U3),U1) 4271 V1 = -(U1 + R) / R 4272 V2 = -U2 / R 4273 V3 = -U3 / R 4274 U2 = V2 / V1 4275 U3 = V3 / V1 4276C 4277 DO 230 I = LOR1, LL 4278 T = A(I,K2) + U2 * A(I,K1) + U3 * A(I,K) 4279 A(I,K2) = A(I,K2) + T * V1 4280 A(I,K1) = A(I,K1) + T * V2 4281 A(I,K) = A(I,K) + T * V3 4282 T = B(I,K2) + U2 * B(I,K1) + U3 * B(I,K) 4283 B(I,K2) = B(I,K2) + T * V1 4284 B(I,K1) = B(I,K1) + T * V2 4285 B(I,K) = B(I,K) + T * V3 4286 230 CONTINUE 4287* --------------------- BEGIN TIMING CODE ----------------------- 4288 OPST = OPST + DBLE( 17 + 20*( LL+1-LOR1 ) ) 4289* ---------------------- END TIMING CODE ------------------------ 4290C 4291 B(K2,K) = 0.0D0 4292 B(K2,K1) = 0.0D0 4293 IF (.NOT. MATZ) GO TO 240 4294C 4295 DO 235 I = 1, N 4296 T = Z(I,K2) + U2 * Z(I,K1) + U3 * Z(I,K) 4297 Z(I,K2) = Z(I,K2) + T * V1 4298 Z(I,K1) = Z(I,K1) + T * V2 4299 Z(I,K) = Z(I,K) + T * V3 4300 235 CONTINUE 4301* --------------------- BEGIN TIMING CODE ----------------------- 4302 OPST = OPST + DBLE( 10*N ) 4303* ---------------------- END TIMING CODE ------------------------ 4304C .......... ZERO B(K+1,K) .......... 4305 240 S = DABS(B(K1,K1)) + DABS(B(K1,K)) 4306 IF (S .EQ. 0.0D0) GO TO 260 4307 U1 = B(K1,K1) / S 4308 U2 = B(K1,K) / S 4309 R = DSIGN(DSQRT(U1*U1+U2*U2),U1) 4310 V1 = -(U1 + R) / R 4311 V2 = -U2 / R 4312 U2 = V2 / V1 4313C 4314 DO 250 I = LOR1, LL 4315 T = A(I,K1) + U2 * A(I,K) 4316 A(I,K1) = A(I,K1) + T * V1 4317 A(I,K) = A(I,K) + T * V2 4318 T = B(I,K1) + U2 * B(I,K) 4319 B(I,K1) = B(I,K1) + T * V1 4320 B(I,K) = B(I,K) + T * V2 4321 250 CONTINUE 4322* --------------------- BEGIN TIMING CODE ----------------------- 4323 OPST = OPST + DBLE( 11 + 12*( LL+1-LOR1 ) ) 4324* ---------------------- END TIMING CODE ------------------------ 4325C 4326 B(K1,K) = 0.0D0 4327 IF (.NOT. MATZ) GO TO 260 4328C 4329 DO 255 I = 1, N 4330 T = Z(I,K1) + U2 * Z(I,K) 4331 Z(I,K1) = Z(I,K1) + T * V1 4332 Z(I,K) = Z(I,K) + T * V2 4333 255 CONTINUE 4334* --------------------- BEGIN TIMING CODE ----------------------- 4335 OPST = OPST + DBLE( 6*N ) 4336* ---------------------- END TIMING CODE ------------------------ 4337C 4338 260 CONTINUE 4339C .......... END QZ STEP .......... 4340 GO TO 70 4341C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT 4342C CONVERGED AFTER 30*N ITERATIONS .......... 4343 1000 IERR = EN 4344C .......... SAVE EPSB FOR USE BY QZVAL AND QZVEC .......... 4345 1001 IF (N .GT. 1) B(N,1) = EPSB 4346* 4347* ---------------------- BEGIN TIMING CODE ------------------------- 4348 OPS = OPS + OPST 4349 OPST = 0.0D0 4350* ----------------------- END TIMING CODE -------------------------- 4351* 4352 RETURN 4353 END 4354 SUBROUTINE QZVAL(NM,N,A,B,ALFR,ALFI,BETA,MATZ,Z) 4355C 4356 INTEGER I,J,N,EN,NA,NM,NN,ISW 4357 DOUBLE PRECISION A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N) 4358 DOUBLE PRECISION C,D,E,R,S,T,AN,A1,A2,BN,CQ,CZ,DI,DR,EI,TI,TR,U1, 4359 X U2,V1,V2,A1I,A11,A12,A2I,A21,A22,B11,B12,B22,SQI,SQR, 4360 X SSI,SSR,SZI,SZR,A11I,A11R,A12I,A12R,A22I,A22R,EPSB 4361 LOGICAL MATZ 4362* 4363* ---------------------- BEGIN TIMING CODE ------------------------- 4364* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 4365* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 4366* OPST IS USED TO ACCUMULATE SMALL CONTRIBUTIONS TO OPS 4367* TO AVOID ROUNDOFF ERROR 4368* .. COMMON BLOCKS .. 4369 COMMON / LATIME / OPS, ITCNT 4370* .. 4371* .. SCALARS IN COMMON .. 4372 DOUBLE PRECISION ITCNT, OPS 4373* .. 4374 DOUBLE PRECISION OPST, OPST2 4375* ----------------------- END TIMING CODE -------------------------- 4376* 4377C 4378C THIS SUBROUTINE IS THE THIRD STEP OF THE QZ ALGORITHM 4379C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, 4380C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. 4381C 4382C THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM 4383C IN QUASI-TRIANGULAR FORM AND THE OTHER IN UPPER TRIANGULAR FORM. 4384C IT REDUCES THE QUASI-TRIANGULAR MATRIX FURTHER, SO THAT ANY 4385C REMAINING 2-BY-2 BLOCKS CORRESPOND TO PAIRS OF COMPLEX 4386C EIGENVALUES, AND RETURNS QUANTITIES WHOSE RATIOS GIVE THE 4387C GENERALIZED EIGENVALUES. IT IS USUALLY PRECEDED BY QZHES 4388C AND QZIT AND MAY BE FOLLOWED BY QZVEC. 4389C 4390C ON INPUT 4391C 4392C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 4393C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 4394C DIMENSION STATEMENT. 4395C 4396C N IS THE ORDER OF THE MATRICES. 4397C 4398C A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX. 4399C 4400C B CONTAINS A REAL UPPER TRIANGULAR MATRIX. IN ADDITION, 4401C LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) 4402C COMPUTED AND SAVED IN QZIT. 4403C 4404C MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS 4405C ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING 4406C EIGENVECTORS, AND TO .FALSE. OTHERWISE. 4407C 4408C Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE 4409C TRANSFORMATION MATRIX PRODUCED IN THE REDUCTIONS BY QZHES 4410C AND QZIT, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. 4411C IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. 4412C 4413C ON OUTPUT 4414C 4415C A HAS BEEN REDUCED FURTHER TO A QUASI-TRIANGULAR MATRIX 4416C IN WHICH ALL NONZERO SUBDIAGONAL ELEMENTS CORRESPOND TO 4417C PAIRS OF COMPLEX EIGENVALUES. 4418C 4419C B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS 4420C HAVE BEEN ALTERED. B(N,1) IS UNALTERED. 4421C 4422C ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS OF THE 4423C DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX THAT WOULD BE 4424C OBTAINED IF A WERE REDUCED COMPLETELY TO TRIANGULAR FORM 4425C BY UNITARY TRANSFORMATIONS. NON-ZERO VALUES OF ALFI OCCUR 4426C IN PAIRS, THE FIRST MEMBER POSITIVE AND THE SECOND NEGATIVE. 4427C 4428C BETA CONTAINS THE DIAGONAL ELEMENTS OF THE CORRESPONDING B, 4429C NORMALIZED TO BE REAL AND NON-NEGATIVE. THE GENERALIZED 4430C EIGENVALUES ARE THEN THE RATIOS ((ALFR+I*ALFI)/BETA). 4431C 4432C Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS 4433C (FOR ALL THREE STEPS) IF MATZ HAS BEEN SET TO .TRUE. 4434C 4435C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 4436C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 4437C 4438C THIS VERSION DATED AUGUST 1983. 4439C 4440C ------------------------------------------------------------------ 4441C 4442 EPSB = B(N,1) 4443 ISW = 1 4444C .......... FIND EIGENVALUES OF QUASI-TRIANGULAR MATRICES. 4445C FOR EN=N STEP -1 UNTIL 1 DO -- .......... 4446* 4447* ---------------------- BEGIN TIMING CODE ------------------------- 4448 OPST = 0.0D0 4449 OPST2 = 0.0D0 4450* ----------------------- END TIMING CODE -------------------------- 4451* 4452 DO 510 NN = 1, N 4453* 4454* --------------------- BEGIN TIMING CODE ----------------------- 4455 OPST = OPST + OPST2 4456 OPST2 = 0.0D0 4457* ---------------------- END TIMING CODE ------------------------ 4458* 4459 EN = N + 1 - NN 4460 NA = EN - 1 4461 IF (ISW .EQ. 2) GO TO 505 4462 IF (EN .EQ. 1) GO TO 410 4463 IF (A(EN,NA) .NE. 0.0D0) GO TO 420 4464C .......... 1-BY-1 BLOCK, ONE REAL ROOT .......... 4465 410 ALFR(EN) = A(EN,EN) 4466 IF (B(EN,EN) .LT. 0.0D0) ALFR(EN) = -ALFR(EN) 4467 BETA(EN) = DABS(B(EN,EN)) 4468 ALFI(EN) = 0.0D0 4469 GO TO 510 4470C .......... 2-BY-2 BLOCK .......... 4471 420 IF (DABS(B(NA,NA)) .LE. EPSB) GO TO 455 4472 IF (DABS(B(EN,EN)) .GT. EPSB) GO TO 430 4473 A1 = A(EN,EN) 4474 A2 = A(EN,NA) 4475 BN = 0.0D0 4476 GO TO 435 4477 430 AN = DABS(A(NA,NA)) + DABS(A(NA,EN)) + DABS(A(EN,NA)) 4478 X + DABS(A(EN,EN)) 4479 BN = DABS(B(NA,NA)) + DABS(B(NA,EN)) + DABS(B(EN,EN)) 4480 A11 = A(NA,NA) / AN 4481 A12 = A(NA,EN) / AN 4482 A21 = A(EN,NA) / AN 4483 A22 = A(EN,EN) / AN 4484 B11 = B(NA,NA) / BN 4485 B12 = B(NA,EN) / BN 4486 B22 = B(EN,EN) / BN 4487 E = A11 / B11 4488 EI = A22 / B22 4489 S = A21 / (B11 * B22) 4490 T = (A22 - E * B22) / B22 4491 IF (DABS(E) .LE. DABS(EI)) GO TO 431 4492 E = EI 4493 T = (A11 - E * B11) / B11 4494 431 C = 0.5D0 * (T - S * B12) 4495 D = C * C + S * (A12 - E * B12) 4496* --------------------- BEGIN TIMING CODE ----------------------- 4497 OPST2 = OPST2 + DBLE( 28 ) 4498* ---------------------- END TIMING CODE ------------------------ 4499 IF (D .LT. 0.0D0) GO TO 480 4500C .......... TWO REAL ROOTS. 4501C ZERO BOTH A(EN,NA) AND B(EN,NA) .......... 4502 E = E + (C + DSIGN(DSQRT(D),C)) 4503 A11 = A11 - E * B11 4504 A12 = A12 - E * B12 4505 A22 = A22 - E * B22 4506* --------------------- BEGIN TIMING CODE ----------------------- 4507 OPST2 = OPST2 + DBLE( 11 ) 4508* ---------------------- END TIMING CODE ------------------------ 4509 IF (DABS(A11) + DABS(A12) .LT. 4510 X DABS(A21) + DABS(A22)) GO TO 432 4511 A1 = A12 4512 A2 = A11 4513 GO TO 435 4514 432 A1 = A22 4515 A2 = A21 4516C .......... CHOOSE AND APPLY REAL Z .......... 4517 435 S = DABS(A1) + DABS(A2) 4518 U1 = A1 / S 4519 U2 = A2 / S 4520 R = DSIGN(DSQRT(U1*U1+U2*U2),U1) 4521 V1 = -(U1 + R) / R 4522 V2 = -U2 / R 4523 U2 = V2 / V1 4524C 4525 DO 440 I = 1, EN 4526 T = A(I,EN) + U2 * A(I,NA) 4527 A(I,EN) = A(I,EN) + T * V1 4528 A(I,NA) = A(I,NA) + T * V2 4529 T = B(I,EN) + U2 * B(I,NA) 4530 B(I,EN) = B(I,EN) + T * V1 4531 B(I,NA) = B(I,NA) + T * V2 4532 440 CONTINUE 4533* --------------------- BEGIN TIMING CODE ----------------------- 4534 OPST2 = OPST2 + DBLE( 11 + 12*EN ) 4535* ---------------------- END TIMING CODE ------------------------ 4536C 4537 IF (.NOT. MATZ) GO TO 450 4538C 4539 DO 445 I = 1, N 4540 T = Z(I,EN) + U2 * Z(I,NA) 4541 Z(I,EN) = Z(I,EN) + T * V1 4542 Z(I,NA) = Z(I,NA) + T * V2 4543 445 CONTINUE 4544* --------------------- BEGIN TIMING CODE ----------------------- 4545 OPST2 = OPST2 + DBLE( 6*N ) 4546* ---------------------- END TIMING CODE ------------------------ 4547C 4548 450 IF (BN .EQ. 0.0D0) GO TO 475 4549 IF (AN .LT. DABS(E) * BN) GO TO 455 4550 A1 = B(NA,NA) 4551 A2 = B(EN,NA) 4552 GO TO 460 4553 455 A1 = A(NA,NA) 4554 A2 = A(EN,NA) 4555C .......... CHOOSE AND APPLY REAL Q .......... 4556 460 S = DABS(A1) + DABS(A2) 4557 IF (S .EQ. 0.0D0) GO TO 475 4558 U1 = A1 / S 4559 U2 = A2 / S 4560 R = DSIGN(DSQRT(U1*U1+U2*U2),U1) 4561 V1 = -(U1 + R) / R 4562 V2 = -U2 / R 4563 U2 = V2 / V1 4564C 4565 DO 470 J = NA, N 4566 T = A(NA,J) + U2 * A(EN,J) 4567 A(NA,J) = A(NA,J) + T * V1 4568 A(EN,J) = A(EN,J) + T * V2 4569 T = B(NA,J) + U2 * B(EN,J) 4570 B(NA,J) = B(NA,J) + T * V1 4571 B(EN,J) = B(EN,J) + T * V2 4572 470 CONTINUE 4573* --------------------- BEGIN TIMING CODE ----------------------- 4574 OPST2 = OPST2 + DBLE( 11 + 12*( N+1-NA ) ) 4575* ---------------------- END TIMING CODE ------------------------ 4576C 4577 475 A(EN,NA) = 0.0D0 4578 B(EN,NA) = 0.0D0 4579 ALFR(NA) = A(NA,NA) 4580 ALFR(EN) = A(EN,EN) 4581 IF (B(NA,NA) .LT. 0.0D0) ALFR(NA) = -ALFR(NA) 4582 IF (B(EN,EN) .LT. 0.0D0) ALFR(EN) = -ALFR(EN) 4583 BETA(NA) = DABS(B(NA,NA)) 4584 BETA(EN) = DABS(B(EN,EN)) 4585 ALFI(EN) = 0.0D0 4586 ALFI(NA) = 0.0D0 4587 GO TO 505 4588C .......... TWO COMPLEX ROOTS .......... 4589 480 E = E + C 4590 EI = DSQRT(-D) 4591 A11R = A11 - E * B11 4592 A11I = EI * B11 4593 A12R = A12 - E * B12 4594 A12I = EI * B12 4595 A22R = A22 - E * B22 4596 A22I = EI * B22 4597 IF (DABS(A11R) + DABS(A11I) + DABS(A12R) + DABS(A12I) .LT. 4598 X DABS(A21) + DABS(A22R) + DABS(A22I)) GO TO 482 4599 A1 = A12R 4600 A1I = A12I 4601 A2 = -A11R 4602 A2I = -A11I 4603 GO TO 485 4604 482 A1 = A22R 4605 A1I = A22I 4606 A2 = -A21 4607 A2I = 0.0D0 4608C .......... CHOOSE COMPLEX Z .......... 4609 485 CZ = DSQRT(A1*A1+A1I*A1I) 4610 IF (CZ .EQ. 0.0D0) GO TO 487 4611 SZR = (A1 * A2 + A1I * A2I) / CZ 4612 SZI = (A1 * A2I - A1I * A2) / CZ 4613 R = DSQRT(CZ*CZ+SZR*SZR+SZI*SZI) 4614 CZ = CZ / R 4615 SZR = SZR / R 4616 SZI = SZI / R 4617 GO TO 490 4618 487 SZR = 1.0D0 4619 SZI = 0.0D0 4620 490 IF (AN .LT. (DABS(E) + EI) * BN) GO TO 492 4621 A1 = CZ * B11 + SZR * B12 4622 A1I = SZI * B12 4623 A2 = SZR * B22 4624 A2I = SZI * B22 4625 GO TO 495 4626 492 A1 = CZ * A11 + SZR * A12 4627 A1I = SZI * A12 4628 A2 = CZ * A21 + SZR * A22 4629 A2I = SZI * A22 4630C .......... CHOOSE COMPLEX Q .......... 4631 495 CQ = DSQRT(A1*A1+A1I*A1I) 4632 IF (CQ .EQ. 0.0D0) GO TO 497 4633 SQR = (A1 * A2 + A1I * A2I) / CQ 4634 SQI = (A1 * A2I - A1I * A2) / CQ 4635 R = DSQRT(CQ*CQ+SQR*SQR+SQI*SQI) 4636 CQ = CQ / R 4637 SQR = SQR / R 4638 SQI = SQI / R 4639 GO TO 500 4640 497 SQR = 1.0D0 4641 SQI = 0.0D0 4642C .......... COMPUTE DIAGONAL ELEMENTS THAT WOULD RESULT 4643C IF TRANSFORMATIONS WERE APPLIED .......... 4644 500 SSR = SQR * SZR + SQI * SZI 4645 SSI = SQR * SZI - SQI * SZR 4646 I = 1 4647 TR = CQ * CZ * A11 + CQ * SZR * A12 + SQR * CZ * A21 4648 X + SSR * A22 4649 TI = CQ * SZI * A12 - SQI * CZ * A21 + SSI * A22 4650 DR = CQ * CZ * B11 + CQ * SZR * B12 + SSR * B22 4651 DI = CQ * SZI * B12 + SSI * B22 4652 GO TO 503 4653 502 I = 2 4654 TR = SSR * A11 - SQR * CZ * A12 - CQ * SZR * A21 4655 X + CQ * CZ * A22 4656 TI = -SSI * A11 - SQI * CZ * A12 + CQ * SZI * A21 4657 DR = SSR * B11 - SQR * CZ * B12 + CQ * CZ * B22 4658 DI = -SSI * B11 - SQI * CZ * B12 4659 503 T = TI * DR - TR * DI 4660 J = NA 4661 IF (T .LT. 0.0D0) J = EN 4662 R = DSQRT(DR*DR+DI*DI) 4663 BETA(J) = BN * R 4664 ALFR(J) = AN * (TR * DR + TI * DI) / R 4665 ALFI(J) = AN * T / R 4666 IF (I .EQ. 1) GO TO 502 4667* --------------------- BEGIN TIMING CODE ----------------------- 4668 OPST2 = OPST2 + DBLE( 151 ) 4669* ---------------------- END TIMING CODE ------------------------ 4670 505 ISW = 3 - ISW 4671 510 CONTINUE 4672* 4673* ---------------------- BEGIN TIMING CODE ------------------------- 4674 OPS = OPS + ( OPST + OPST2 ) 4675* ----------------------- END TIMING CODE -------------------------- 4676* 4677 B(N,1) = EPSB 4678C 4679 RETURN 4680 END 4681 SUBROUTINE QZVEC(NM,N,A,B,ALFR,ALFI,BETA,Z) 4682C 4683 INTEGER I,J,K,M,N,EN,II,JJ,NA,NM,NN,ISW,ENM2 4684 DOUBLE PRECISION A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N) 4685 DOUBLE PRECISION D,Q,R,S,T,W,X,Y,DI,DR,RA,RR,SA,TI,TR,T1,T2,W1,X1, 4686 X ZZ,Z1,ALFM,ALMI,ALMR,BETM,EPSB 4687* 4688* ---------------------- BEGIN TIMING CODE ------------------------- 4689* COMMON BLOCK TO RETURN OPERATION COUNT AND ITERATION COUNT 4690* ITCNT IS INITIALIZED TO 0, OPS IS ONLY INCREMENTED 4691* OPST IS USED TO ACCUMULATE SMALL CONTRIBUTIONS TO OPS 4692* TO AVOID ROUNDOFF ERROR 4693* .. COMMON BLOCKS .. 4694 COMMON / LATIME / OPS, ITCNT 4695* .. 4696* .. SCALARS IN COMMON .. 4697 DOUBLE PRECISION ITCNT, OPS 4698* .. 4699 INTEGER IN2BY2 4700* ----------------------- END TIMING CODE -------------------------- 4701* 4702C 4703C THIS SUBROUTINE IS THE OPTIONAL FOURTH STEP OF THE QZ ALGORITHM 4704C FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, 4705C SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. 4706C 4707C THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM IN 4708C QUASI-TRIANGULAR FORM (IN WHICH EACH 2-BY-2 BLOCK CORRESPONDS TO 4709C A PAIR OF COMPLEX EIGENVALUES) AND THE OTHER IN UPPER TRIANGULAR 4710C FORM. IT COMPUTES THE EIGENVECTORS OF THE TRIANGULAR PROBLEM AND 4711C TRANSFORMS THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM. 4712C IT IS USUALLY PRECEDED BY QZHES, QZIT, AND QZVAL. 4713C 4714C ON INPUT 4715C 4716C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 4717C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 4718C DIMENSION STATEMENT. 4719C 4720C N IS THE ORDER OF THE MATRICES. 4721C 4722C A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX. 4723C 4724C B CONTAINS A REAL UPPER TRIANGULAR MATRIX. IN ADDITION, 4725C LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) 4726C COMPUTED AND SAVED IN QZIT. 4727C 4728C ALFR, ALFI, AND BETA ARE VECTORS WITH COMPONENTS WHOSE 4729C RATIOS ((ALFR+I*ALFI)/BETA) ARE THE GENERALIZED 4730C EIGENVALUES. THEY ARE USUALLY OBTAINED FROM QZVAL. 4731C 4732C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE 4733C REDUCTIONS BY QZHES, QZIT, AND QZVAL, IF PERFORMED. 4734C IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM ARE 4735C DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX. 4736C 4737C ON OUTPUT 4738C 4739C A IS UNALTERED. ITS SUBDIAGONAL ELEMENTS PROVIDE INFORMATION 4740C ABOUT THE STORAGE OF THE COMPLEX EIGENVECTORS. 4741C 4742C B HAS BEEN DESTROYED. 4743C 4744C ALFR, ALFI, AND BETA ARE UNALTERED. 4745C 4746C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. 4747C IF ALFI(I) .EQ. 0.0, THE I-TH EIGENVALUE IS REAL AND 4748C THE I-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR. 4749C IF ALFI(I) .NE. 0.0, THE I-TH EIGENVALUE IS COMPLEX. 4750C IF ALFI(I) .GT. 0.0, THE EIGENVALUE IS THE FIRST OF 4751C A COMPLEX PAIR AND THE I-TH AND (I+1)-TH COLUMNS 4752C OF Z CONTAIN ITS EIGENVECTOR. 4753C IF ALFI(I) .LT. 0.0, THE EIGENVALUE IS THE SECOND OF 4754C A COMPLEX PAIR AND THE (I-1)-TH AND I-TH COLUMNS 4755C OF Z CONTAIN THE CONJUGATE OF ITS EIGENVECTOR. 4756C EACH EIGENVECTOR IS NORMALIZED SO THAT THE MODULUS 4757C OF ITS LARGEST COMPONENT IS 1.0 . 4758C 4759C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 4760C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 4761C 4762C THIS VERSION DATED AUGUST 1983. 4763C 4764C ------------------------------------------------------------------ 4765C 4766 EPSB = B(N,1) 4767 ISW = 1 4768C .......... FOR EN=N STEP -1 UNTIL 1 DO -- .......... 4769 DO 800 NN = 1, N 4770* --------------------- BEGIN TIMING CODE ----------------------- 4771 IN2BY2 = 0 4772* ---------------------- END TIMING CODE ------------------------ 4773 EN = N + 1 - NN 4774 NA = EN - 1 4775 IF (ISW .EQ. 2) GO TO 795 4776 IF (ALFI(EN) .NE. 0.0D0) GO TO 710 4777C .......... REAL VECTOR .......... 4778 M = EN 4779 B(EN,EN) = 1.0D0 4780 IF (NA .EQ. 0) GO TO 800 4781 ALFM = ALFR(M) 4782 BETM = BETA(M) 4783C .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... 4784 DO 700 II = 1, NA 4785 I = EN - II 4786 W = BETM * A(I,I) - ALFM * B(I,I) 4787 R = 0.0D0 4788C 4789 DO 610 J = M, EN 4790 610 R = R + (BETM * A(I,J) - ALFM * B(I,J)) * B(J,EN) 4791C 4792 IF (I .EQ. 1 .OR. ISW .EQ. 2) GO TO 630 4793 IF (BETM * A(I,I-1) .EQ. 0.0D0) GO TO 630 4794 ZZ = W 4795 S = R 4796 GO TO 690 4797 630 M = I 4798 IF (ISW .EQ. 2) GO TO 640 4799C .......... REAL 1-BY-1 BLOCK .......... 4800 T = W 4801 IF (W .EQ. 0.0D0) T = EPSB 4802 B(I,EN) = -R / T 4803 GO TO 700 4804C .......... REAL 2-BY-2 BLOCK .......... 4805 640 X = BETM * A(I,I+1) - ALFM * B(I,I+1) 4806 Y = BETM * A(I+1,I) 4807 Q = W * ZZ - X * Y 4808 T = (X * S - ZZ * R) / Q 4809 B(I,EN) = T 4810* ------------------- BEGIN TIMING CODE ---------------------- 4811 IN2BY2 = IN2BY2 + 1 4812* -------------------- END TIMING CODE ----------------------- 4813 IF (DABS(X) .LE. DABS(ZZ)) GO TO 650 4814 B(I+1,EN) = (-R - W * T) / X 4815 GO TO 690 4816 650 B(I+1,EN) = (-S - Y * T) / ZZ 4817 690 ISW = 3 - ISW 4818 700 CONTINUE 4819C .......... END REAL VECTOR .......... 4820* --------------------- BEGIN TIMING CODE ----------------------- 4821 OPS = OPS + ( 5.0D0/2.0D0 )*DBLE( ( EN+2 )*( EN-1 ) + IN2BY2 ) 4822* ---------------------- END TIMING CODE ------------------------ 4823 GO TO 800 4824C .......... COMPLEX VECTOR .......... 4825 710 M = NA 4826 ALMR = ALFR(M) 4827 ALMI = ALFI(M) 4828 BETM = BETA(M) 4829C .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT 4830C EIGENVECTOR MATRIX IS TRIANGULAR .......... 4831 Y = BETM * A(EN,NA) 4832 B(NA,NA) = -ALMI * B(EN,EN) / Y 4833 B(NA,EN) = (ALMR * B(EN,EN) - BETM * A(EN,EN)) / Y 4834 B(EN,NA) = 0.0D0 4835 B(EN,EN) = 1.0D0 4836 ENM2 = NA - 1 4837 IF (ENM2 .EQ. 0) GO TO 795 4838C .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- .......... 4839 DO 790 II = 1, ENM2 4840 I = NA - II 4841 W = BETM * A(I,I) - ALMR * B(I,I) 4842 W1 = -ALMI * B(I,I) 4843 RA = 0.0D0 4844 SA = 0.0D0 4845C 4846 DO 760 J = M, EN 4847 X = BETM * A(I,J) - ALMR * B(I,J) 4848 X1 = -ALMI * B(I,J) 4849 RA = RA + X * B(J,NA) - X1 * B(J,EN) 4850 SA = SA + X * B(J,EN) + X1 * B(J,NA) 4851 760 CONTINUE 4852C 4853 IF (I .EQ. 1 .OR. ISW .EQ. 2) GO TO 770 4854 IF (BETM * A(I,I-1) .EQ. 0.0D0) GO TO 770 4855 ZZ = W 4856 Z1 = W1 4857 R = RA 4858 S = SA 4859 ISW = 2 4860 GO TO 790 4861 770 M = I 4862 IF (ISW .EQ. 2) GO TO 780 4863C .......... COMPLEX 1-BY-1 BLOCK .......... 4864 TR = -RA 4865 TI = -SA 4866 773 DR = W 4867 DI = W1 4868C .......... COMPLEX DIVIDE (T1,T2) = (TR,TI) / (DR,DI) .......... 4869 775 IF (DABS(DI) .GT. DABS(DR)) GO TO 777 4870 RR = DI / DR 4871 D = DR + DI * RR 4872 T1 = (TR + TI * RR) / D 4873 T2 = (TI - TR * RR) / D 4874 GO TO (787,782), ISW 4875 777 RR = DR / DI 4876 D = DR * RR + DI 4877 T1 = (TR * RR + TI) / D 4878 T2 = (TI * RR - TR) / D 4879 GO TO (787,782), ISW 4880C .......... COMPLEX 2-BY-2 BLOCK .......... 4881 780 X = BETM * A(I,I+1) - ALMR * B(I,I+1) 4882 X1 = -ALMI * B(I,I+1) 4883 Y = BETM * A(I+1,I) 4884 TR = Y * RA - W * R + W1 * S 4885 TI = Y * SA - W * S - W1 * R 4886 DR = W * ZZ - W1 * Z1 - X * Y 4887 DI = W * Z1 + W1 * ZZ - X1 * Y 4888* ------------------- BEGIN TIMING CODE ---------------------- 4889 IN2BY2 = IN2BY2 + 1 4890* -------------------- END TIMING CODE ----------------------- 4891 IF (DR .EQ. 0.0D0 .AND. DI .EQ. 0.0D0) DR = EPSB 4892 GO TO 775 4893 782 B(I+1,NA) = T1 4894 B(I+1,EN) = T2 4895 ISW = 1 4896 IF (DABS(Y) .GT. DABS(W) + DABS(W1)) GO TO 785 4897 TR = -RA - X * B(I+1,NA) + X1 * B(I+1,EN) 4898 TI = -SA - X * B(I+1,EN) - X1 * B(I+1,NA) 4899 GO TO 773 4900 785 T1 = (-R - ZZ * B(I+1,NA) + Z1 * B(I+1,EN)) / Y 4901 T2 = (-S - ZZ * B(I+1,EN) - Z1 * B(I+1,NA)) / Y 4902 787 B(I,NA) = T1 4903 B(I,EN) = T2 4904 790 CONTINUE 4905* --------------------- BEGIN TIMING CODE ----------------------- 4906 OPS = OPS + DBLE( ( 6*EN-7 )*( EN-2 ) + 31*IN2BY2 ) 4907* ---------------------- END TIMING CODE ------------------------ 4908C .......... END COMPLEX VECTOR .......... 4909 795 ISW = 3 - ISW 4910 800 CONTINUE 4911C .......... END BACK SUBSTITUTION. 4912C TRANSFORM TO ORIGINAL COORDINATE SYSTEM. 4913C FOR J=N STEP -1 UNTIL 1 DO -- .......... 4914 DO 880 JJ = 1, N 4915 J = N + 1 - JJ 4916C 4917 DO 880 I = 1, N 4918 ZZ = 0.0D0 4919C 4920 DO 860 K = 1, J 4921 860 ZZ = ZZ + Z(I,K) * B(K,J) 4922C 4923 Z(I,J) = ZZ 4924 880 CONTINUE 4925* ----------------------- BEGIN TIMING CODE ------------------------ 4926 OPS = OPS + DBLE( N**2 )*DBLE( N+1 ) 4927* ------------------------ END TIMING CODE ------------------------- 4928C .......... NORMALIZE SO THAT MODULUS OF LARGEST 4929C COMPONENT OF EACH VECTOR IS 1. 4930C (ISW IS 1 INITIALLY FROM BEFORE) .......... 4931* ------------------------ BEGIN TIMING CODE ----------------------- 4932 IN2BY2 = 0 4933* ------------------------- END TIMING CODE ------------------------ 4934 DO 950 J = 1, N 4935 D = 0.0D0 4936 IF (ISW .EQ. 2) GO TO 920 4937 IF (ALFI(J) .NE. 0.0D0) GO TO 945 4938C 4939 DO 890 I = 1, N 4940 IF (DABS(Z(I,J)) .GT. D) D = DABS(Z(I,J)) 4941 890 CONTINUE 4942C 4943 DO 900 I = 1, N 4944 900 Z(I,J) = Z(I,J) / D 4945C 4946 GO TO 950 4947C 4948 920 DO 930 I = 1, N 4949 R = DABS(Z(I,J-1)) + DABS(Z(I,J)) 4950 IF (R .NE. 0.0D0) R = R * DSQRT((Z(I,J-1)/R)**2 4951 X +(Z(I,J)/R)**2) 4952 IF (R .GT. D) D = R 4953 930 CONTINUE 4954C 4955 DO 940 I = 1, N 4956 Z(I,J-1) = Z(I,J-1) / D 4957 Z(I,J) = Z(I,J) / D 4958 940 CONTINUE 4959* ---------------------- BEGIN TIMING CODE ---------------------- 4960 IN2BY2 = IN2BY2 + 1 4961* ----------------------- END TIMING CODE ----------------------- 4962C 4963 945 ISW = 3 - ISW 4964 950 CONTINUE 4965* ------------------------ BEGIN TIMING CODE ----------------------- 4966 OPS = OPS + DBLE( N*( N + 5*IN2BY2 ) ) 4967* ------------------------- END TIMING CODE ------------------------ 4968C 4969 RETURN 4970 END 4971