pythiaeRHIC/pythia/pycmqr.f

0001
0002 C*********************************************************************
0003
0004 C...PYCMQR
0005 C...Auxiliary to PYEICG.
0006 C
0007 C     THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
0008 C     ALGOL PROCEDURE  COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN
0009 C     AND WILKINSON.
0010 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971).
0011 C     THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS
0012 C     (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM.
0013 C
0014 C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX
0015 C     UPPER HESSENBERG MATRIX BY THE QR METHOD.
0016 C
0017 C     ON INPUT
0018 C
0019 C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
0020 C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
0021 C          DIMENSION STATEMENT.
0022 C
0023 C        N IS THE ORDER OF THE MATRIX.
0024 C
0025 C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
0026 C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
0027 C          SET LOW=1, IGH=N.
0028 C
0029 C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
0030 C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
0031 C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN
0032 C          INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN
0033 C          THE REDUCTION BY  CORTH, IF PERFORMED.
0034 C
0035 C     ON OUTPUT
0036 C
0037 C        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN
0038 C          DESTROYED.  THEREFORE, THEY MUST BE SAVED BEFORE
0039 C          CALLING  COMQR  IF SUBSEQUENT CALCULATION OF
0040 C          EIGENVECTORS IS TO BE PERFORMED.
0041 C
0042 C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
0043 C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR
0044 C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
0045 C          FOR INDICES IERR+1,...,N.
0046 C
0047 C        IERR IS SET TO
0048 C          ZERO       FOR NORMAL RETURN,
0049 C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
0050 C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
0051 C
0052 C     CALLS PYCDIV FOR COMPLEX DIVISION.
0053 C     CALLS PYCSRT FOR COMPLEX SQUARE ROOT.
0054 C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
0055 C
0056 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
0057 C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
0058 C
0059 C     THIS VERSION DATED AUGUST 1983.
0060 C
0061
0062       SUBROUTINE PYCMQR(NM,N,LOW,IGH,HR,HI,WR,WI,IERR)
0063
0064       INTEGER I,J,L,N,EN,LL,NM,IGH,ITN,ITS,LOW,LP1,ENM1,IERR
0065       DOUBLE PRECISION HR(4,4),HI(4,4),WR(4),WI(4)
0066       DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2,
0067      X       PYTHAG
0068
0069       IERR = 0
0070       IF (LOW .EQ. IGH) GOTO 130
0071 C     .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
0072       L = LOW + 1
0073 C
0074       DO 120 I = L, IGH
0075          LL = MIN0(I+1,IGH)
0076          IF (HI(I,I-1) .EQ. 0.0D0) GOTO 120
0077          NORM = PYTHAG(HR(I,I-1),HI(I,I-1))
0078          YR = HR(I,I-1) / NORM
0079          YI = HI(I,I-1) / NORM
0080          HR(I,I-1) = NORM
0081          HI(I,I-1) = 0.0D0
0082 C
0083          DO 100 J = I, IGH
0084             SI = YR * HI(I,J) - YI * HR(I,J)
0085             HR(I,J) = YR * HR(I,J) + YI * HI(I,J)
0086             HI(I,J) = SI
0087   100    CONTINUE
0088 C
0089          DO 110 J = LOW, LL
0090             SI = YR * HI(J,I) + YI * HR(J,I)
0091             HR(J,I) = YR * HR(J,I) - YI * HI(J,I)
0092             HI(J,I) = SI
0093   110    CONTINUE
0094 C
0095   120 CONTINUE
0096 C     .......... STORE ROOTS ISOLATED BY CBAL ..........
0097   130 DO 140 I = 1, N
0098          IF (I .GE. LOW .AND. I .LE. IGH) GOTO 140
0099          WR(I) = HR(I,I)
0100          WI(I) = HI(I,I)
0101   140 CONTINUE
0102 C
0103       EN = IGH
0104       TR = 0.0D0
0105       TI = 0.0D0
0106       ITN = 30*N
0107 C     .......... SEARCH FOR NEXT EIGENVALUE ..........
0108   150 IF (EN .LT. LOW) GOTO 320
0109       ITS = 0
0110       ENM1 = EN - 1
0111 C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
0112 C                FOR L=EN STEP -1 UNTIL LOW D0 -- ..........
0113   160 DO 170 LL = LOW, EN
0114          L = EN + LOW - LL
0115          IF (L .EQ. LOW) GOTO 180
0116          TST1 = DABS(HR(L-1,L-1)) + DABS(HI(L-1,L-1))
0117      X            + DABS(HR(L,L)) + DABS(HI(L,L))
0118          TST2 = TST1 + DABS(HR(L,L-1))
0119          IF (TST2 .EQ. TST1) GOTO 180
0120   170 CONTINUE
0121 C     .......... FORM SHIFT ..........
0122   180 IF (L .EQ. EN) GOTO 300
0123       IF (ITN .EQ. 0) GOTO 310
0124       IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GOTO 200
0125       SR = HR(EN,EN)
0126       SI = HI(EN,EN)
0127       XR = HR(ENM1,EN) * HR(EN,ENM1)
0128       XI = HI(ENM1,EN) * HR(EN,ENM1)
0129       IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GOTO 210
0130       YR = (HR(ENM1,ENM1) - SR) / 2.0D0
0131       YI = (HI(ENM1,ENM1) - SI) / 2.0D0
0132       CALL PYCSRT(YR**2-YI**2+XR,2.0D0*YR*YI+XI,ZZR,ZZI)
0133       IF (YR * ZZR + YI * ZZI .GE. 0.0D0) GOTO 190
0134       ZZR = -ZZR
0135       ZZI = -ZZI
0136   190 CALL PYCDIV(XR,XI,YR+ZZR,YI+ZZI,XR,XI)
0137       SR = SR - XR
0138       SI = SI - XI
0139       GOTO 210
0140 C     .......... FORM EXCEPTIONAL SHIFT ..........
0141   200 SR = DABS(HR(EN,ENM1)) + DABS(HR(ENM1,EN-2))
0142       SI = 0.0D0
0143 C
0144   210 DO 220 I = LOW, EN
0145          HR(I,I) = HR(I,I) - SR
0146          HI(I,I) = HI(I,I) - SI
0147   220 CONTINUE
0148 C
0149       TR = TR + SR
0150       TI = TI + SI
0151       ITS = ITS + 1
0152       ITN = ITN - 1
0153 C     .......... REDUCE TO TRIANGLE (ROWS) ..........
0154       LP1 = L + 1
0155 C
0156       DO 240 I = LP1, EN
0157          SR = HR(I,I-1)
0158          HR(I,I-1) = 0.0D0
0159          NORM = PYTHAG(PYTHAG(HR(I-1,I-1),HI(I-1,I-1)),SR)
0160          XR = HR(I-1,I-1) / NORM
0161          WR(I-1) = XR
0162          XI = HI(I-1,I-1) / NORM
0163          WI(I-1) = XI
0164          HR(I-1,I-1) = NORM
0165          HI(I-1,I-1) = 0.0D0
0166          HI(I,I-1) = SR / NORM
0167 C
0168          DO 230 J = I, EN
0169             YR = HR(I-1,J)
0170             YI = HI(I-1,J)
0171             ZZR = HR(I,J)
0172             ZZI = HI(I,J)
0173             HR(I-1,J) = XR * YR + XI * YI + HI(I,I-1) * ZZR
0174             HI(I-1,J) = XR * YI - XI * YR + HI(I,I-1) * ZZI
0175             HR(I,J) = XR * ZZR - XI * ZZI - HI(I,I-1) * YR
0176             HI(I,J) = XR * ZZI + XI * ZZR - HI(I,I-1) * YI
0177   230    CONTINUE
0178 C
0179   240 CONTINUE
0180 C
0181       SI = HI(EN,EN)
0182       IF (SI .EQ. 0.0D0) GOTO 250
0183       NORM = PYTHAG(HR(EN,EN),SI)
0184       SR = HR(EN,EN) / NORM
0185       SI = SI / NORM
0186       HR(EN,EN) = NORM
0187       HI(EN,EN) = 0.0D0
0188 C     .......... INVERSE OPERATION (COLUMNS) ..........
0189   250 DO 280 J = LP1, EN
0190          XR = WR(J-1)
0191          XI = WI(J-1)
0192 C
0193          DO 270 I = L, J
0194             YR = HR(I,J-1)
0195             YI = 0.0D0
0196             ZZR = HR(I,J)
0197             ZZI = HI(I,J)
0198             IF (I .EQ. J) GOTO 260
0199             YI = HI(I,J-1)
0200             HI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
0201   260       HR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
0202             HR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
0203             HI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
0204   270    CONTINUE
0205 C
0206   280 CONTINUE
0207 C
0208       IF (SI .EQ. 0.0D0) GOTO 160
0209 C
0210       DO 290 I = L, EN
0211          YR = HR(I,EN)
0212          YI = HI(I,EN)
0213          HR(I,EN) = SR * YR - SI * YI
0214          HI(I,EN) = SR * YI + SI * YR
0215   290 CONTINUE
0216 C
0217       GOTO 160
0218 C     .......... A ROOT FOUND ..........
0219   300 WR(EN) = HR(EN,EN) + TR
0220       WI(EN) = HI(EN,EN) + TI
0221       EN = ENM1
0222       GOTO 150
0223 C     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
0224 C                CONVERGED AFTER 30*N ITERATIONS ..........
0225   310 IERR = EN
0226   320 RETURN
0227       END