pythiaeRHIC/pythia/pycmq2.f

0001
0002 C*********************************************************************
0003
0004 C...PYCMQ2
0005 C...Auxiliary to PYEICG.
0006 C
0007 C     THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
0008 C     ALGOL PROCEDURE  COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS
0009 C     AND WILKINSON.
0010 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(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 AND EIGENVECTORS
0015 C     OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR
0016 C     METHOD.  THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX
0017 C     CAN ALSO BE FOUND IF  CORTH  HAS BEEN USED TO REDUCE
0018 C     THIS GENERAL MATRIX TO HESSENBERG FORM.
0019 C
0020 C     ON INPUT
0021 C
0022 C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
0023 C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
0024 C          DIMENSION STATEMENT.
0025 C
0026 C        N IS THE ORDER OF THE MATRIX.
0027 C
0028 C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
0029 C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
0030 C          SET LOW=1, IGH=N.
0031 C
0032 C        ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
0033 C          FORMATIONS USED IN THE REDUCTION BY  CORTH, IF PERFORMED.
0034 C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.  IF THE EIGENVECTORS
0035 C          OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND
0036 C          ORTI(J) TO 0.0D0 FOR THESE ELEMENTS.
0037 C
0038 C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
0039 C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
0040 C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER
0041 C          INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE
0042 C          REDUCTION BY  CORTH, IF PERFORMED.  IF THE EIGENVECTORS OF
0043 C          THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE
0044 C          ARBITRARY.
0045 C
0046 C     ON OUTPUT
0047 C
0048 C        ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI
0049 C          HAVE BEEN DESTROYED.
0050 C
0051 C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
0052 C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR
0053 C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
0054 C          FOR INDICES IERR+1,...,N.
0055 C
0056 C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
0057 C          RESPECTIVELY, OF THE EIGENVECTORS.  THE EIGENVECTORS
0058 C          ARE UNNORMALIZED.  IF AN ERROR EXIT IS MADE, NONE OF
0059 C          THE EIGENVECTORS HAS BEEN FOUND.
0060 C
0061 C        IERR IS SET TO
0062 C          ZERO       FOR NORMAL RETURN,
0063 C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
0064 C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
0065 C
0066 C     CALLS PYCDIV FOR COMPLEX DIVISION.
0067 C     CALLS PYCSRT FOR COMPLEX SQUARE ROOT.
0068 C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
0069 C
0070 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
0071 C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
0072 C
0073 C     THIS VERSION DATED OCTOBER 1989.
0074 C
0075 C  MESHED OVERFLOW CONTROL WITH VECTORS OF ISOLATED ROOTS (10/19/89 BSG)
0076 C  MESHED OVERFLOW CONTROL WITH TRIANGULAR MULTIPLY (10/30/89 BSG)
0077 C
0078
0079       SUBROUTINE PYCMQ2(NM,N,LOW,IGH,ORTR,ORTI,HR,HI,WR,WI,ZR,ZI,IERR)
0080
0081       INTEGER I,J,K,L,M,N,EN,II,JJ,LL,NM,NN,IGH,IP1,
0082      X        ITN,ITS,LOW,LP1,ENM1,IEND,IERR
0083       DOUBLE PRECISION HR(4,4),HI(4,4),WR(4),WI(4),ZR(4,4),ZI(4,4),
0084      X       ORTR(4),ORTI(4)
0085       DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2,
0086      X       PYTHAG
0087
0088       IERR = 0
0089 C     .......... INITIALIZE EIGENVECTOR MATRIX ..........
0090       DO 110 J = 1, N
0091 C
0092          DO 100 I = 1, N
0093             ZR(I,J) = 0.0D0
0094             ZI(I,J) = 0.0D0
0095   100    CONTINUE
0096          ZR(J,J) = 1.0D0
0097   110 CONTINUE
0098 C     .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS
0099 C                FROM THE INFORMATION LEFT BY CORTH ..........
0100       IEND = IGH - LOW - 1
0101       IF (IEND.LT.0) GOTO 220
0102       IF (IEND.EQ.0) GOTO 170
0103 C     .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
0104       DO 160 II = 1, IEND
0105          I = IGH - II
0106          IF (ORTR(I) .EQ. 0.0D0 .AND. ORTI(I) .EQ. 0.0D0) GOTO 160
0107          IF (HR(I,I-1) .EQ. 0.0D0 .AND. HI(I,I-1) .EQ. 0.0D0) GOTO 160
0108 C     .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ..........
0109          NORM = HR(I,I-1) * ORTR(I) + HI(I,I-1) * ORTI(I)
0110          IP1 = I + 1
0111 C
0112          DO 120 K = IP1, IGH
0113             ORTR(K) = HR(K,I-1)
0114             ORTI(K) = HI(K,I-1)
0115   120    CONTINUE
0116 C
0117          DO 150 J = I, IGH
0118             SR = 0.0D0
0119             SI = 0.0D0
0120 C
0121             DO 130 K = I, IGH
0122                SR = SR + ORTR(K) * ZR(K,J) + ORTI(K) * ZI(K,J)
0123                SI = SI + ORTR(K) * ZI(K,J) - ORTI(K) * ZR(K,J)
0124   130       CONTINUE
0125 C
0126             SR = SR / NORM
0127             SI = SI / NORM
0128 C
0129             DO 140 K = I, IGH
0130                ZR(K,J) = ZR(K,J) + SR * ORTR(K) - SI * ORTI(K)
0131                ZI(K,J) = ZI(K,J) + SR * ORTI(K) + SI * ORTR(K)
0132   140       CONTINUE
0133 C
0134   150    CONTINUE
0135 C
0136   160 CONTINUE
0137 C     .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
0138   170 L = LOW + 1
0139 C
0140       DO 210 I = L, IGH
0141          LL = MIN0(I+1,IGH)
0142          IF (HI(I,I-1) .EQ. 0.0D0) GOTO 210
0143          NORM = PYTHAG(HR(I,I-1),HI(I,I-1))
0144          YR = HR(I,I-1) / NORM
0145          YI = HI(I,I-1) / NORM
0146          HR(I,I-1) = NORM
0147          HI(I,I-1) = 0.0D0
0148 C
0149          DO 180 J = I, N
0150             SI = YR * HI(I,J) - YI * HR(I,J)
0151             HR(I,J) = YR * HR(I,J) + YI * HI(I,J)
0152             HI(I,J) = SI
0153   180    CONTINUE
0154 C
0155          DO 190 J = 1, LL
0156             SI = YR * HI(J,I) + YI * HR(J,I)
0157             HR(J,I) = YR * HR(J,I) - YI * HI(J,I)
0158             HI(J,I) = SI
0159   190    CONTINUE
0160 C
0161          DO 200 J = LOW, IGH
0162             SI = YR * ZI(J,I) + YI * ZR(J,I)
0163             ZR(J,I) = YR * ZR(J,I) - YI * ZI(J,I)
0164             ZI(J,I) = SI
0165   200    CONTINUE
0166 C
0167   210 CONTINUE
0168 C     .......... STORE ROOTS ISOLATED BY CBAL ..........
0169   220 DO 230 I = 1, N
0170          IF (I .GE. LOW .AND. I .LE. IGH) GOTO 230
0171          WR(I) = HR(I,I)
0172          WI(I) = HI(I,I)
0173   230 CONTINUE
0174 C
0175       EN = IGH
0176       TR = 0.0D0
0177       TI = 0.0D0
0178       ITN = 30*N
0179 C     .......... SEARCH FOR NEXT EIGENVALUE ..........
0180   240 IF (EN .LT. LOW) GOTO 430
0181       ITS = 0
0182       ENM1 = EN - 1
0183 C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
0184 C                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
0185   250 DO 260 LL = LOW, EN
0186          L = EN + LOW - LL
0187          IF (L .EQ. LOW) GOTO 270
0188          TST1 = DABS(HR(L-1,L-1)) + DABS(HI(L-1,L-1))
0189      X            + DABS(HR(L,L)) + DABS(HI(L,L))
0190          TST2 = TST1 + DABS(HR(L,L-1))
0191          IF (TST2 .EQ. TST1) GOTO 270
0192   260 CONTINUE
0193 C     .......... FORM SHIFT ..........
0194   270 IF (L .EQ. EN) GOTO 420
0195       IF (ITN .EQ. 0) GOTO 550
0196       IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GOTO 290
0197       SR = HR(EN,EN)
0198       SI = HI(EN,EN)
0199       XR = HR(ENM1,EN) * HR(EN,ENM1)
0200       XI = HI(ENM1,EN) * HR(EN,ENM1)
0201       IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GOTO 300
0202       YR = (HR(ENM1,ENM1) - SR) / 2.0D0
0203       YI = (HI(ENM1,ENM1) - SI) / 2.0D0
0204       CALL PYCSRT(YR**2-YI**2+XR,2.0D0*YR*YI+XI,ZZR,ZZI)
0205       IF (YR * ZZR + YI * ZZI .GE. 0.0D0) GOTO 280
0206       ZZR = -ZZR
0207       ZZI = -ZZI
0208   280 CALL PYCDIV(XR,XI,YR+ZZR,YI+ZZI,XR,XI)
0209       SR = SR - XR
0210       SI = SI - XI
0211       GOTO 300
0212 C     .......... FORM EXCEPTIONAL SHIFT ..........
0213   290 SR = DABS(HR(EN,ENM1)) + DABS(HR(ENM1,EN-2))
0214       SI = 0.0D0
0215 C
0216   300 DO 310 I = LOW, EN
0217          HR(I,I) = HR(I,I) - SR
0218          HI(I,I) = HI(I,I) - SI
0219   310 CONTINUE
0220 C
0221       TR = TR + SR
0222       TI = TI + SI
0223       ITS = ITS + 1
0224       ITN = ITN - 1
0225 C     .......... REDUCE TO TRIANGLE (ROWS) ..........
0226       LP1 = L + 1
0227 C
0228       DO 330 I = LP1, EN
0229          SR = HR(I,I-1)
0230          HR(I,I-1) = 0.0D0
0231          NORM = PYTHAG(PYTHAG(HR(I-1,I-1),HI(I-1,I-1)),SR)
0232          XR = HR(I-1,I-1) / NORM
0233          WR(I-1) = XR
0234          XI = HI(I-1,I-1) / NORM
0235          WI(I-1) = XI
0236          HR(I-1,I-1) = NORM
0237          HI(I-1,I-1) = 0.0D0
0238          HI(I,I-1) = SR / NORM
0239 C
0240          DO 320 J = I, N
0241             YR = HR(I-1,J)
0242             YI = HI(I-1,J)
0243             ZZR = HR(I,J)
0244             ZZI = HI(I,J)
0245             HR(I-1,J) = XR * YR + XI * YI + HI(I,I-1) * ZZR
0246             HI(I-1,J) = XR * YI - XI * YR + HI(I,I-1) * ZZI
0247             HR(I,J) = XR * ZZR - XI * ZZI - HI(I,I-1) * YR
0248             HI(I,J) = XR * ZZI + XI * ZZR - HI(I,I-1) * YI
0249   320    CONTINUE
0250 C
0251   330 CONTINUE
0252 C
0253       SI = HI(EN,EN)
0254       IF (SI .EQ. 0.0D0) GOTO 350
0255       NORM = PYTHAG(HR(EN,EN),SI)
0256       SR = HR(EN,EN) / NORM
0257       SI = SI / NORM
0258       HR(EN,EN) = NORM
0259       HI(EN,EN) = 0.0D0
0260       IF (EN .EQ. N) GOTO 350
0261       IP1 = EN + 1
0262 C
0263       DO 340 J = IP1, N
0264          YR = HR(EN,J)
0265          YI = HI(EN,J)
0266          HR(EN,J) = SR * YR + SI * YI
0267          HI(EN,J) = SR * YI - SI * YR
0268   340 CONTINUE
0269 C     .......... INVERSE OPERATION (COLUMNS) ..........
0270   350 DO 390 J = LP1, EN
0271          XR = WR(J-1)
0272          XI = WI(J-1)
0273 C
0274          DO 370 I = 1, J
0275             YR = HR(I,J-1)
0276             YI = 0.0D0
0277             ZZR = HR(I,J)
0278             ZZI = HI(I,J)
0279             IF (I .EQ. J) GOTO 360
0280             YI = HI(I,J-1)
0281             HI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
0282   360       HR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
0283             HR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
0284             HI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
0285   370    CONTINUE
0286 C
0287          DO 380 I = LOW, IGH
0288             YR = ZR(I,J-1)
0289             YI = ZI(I,J-1)
0290             ZZR = ZR(I,J)
0291             ZZI = ZI(I,J)
0292             ZR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
0293             ZI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
0294             ZR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
0295             ZI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
0296   380    CONTINUE
0297 C
0298   390 CONTINUE
0299 C
0300       IF (SI .EQ. 0.0D0) GOTO 250
0301 C
0302       DO 400 I = 1, EN
0303          YR = HR(I,EN)
0304          YI = HI(I,EN)
0305          HR(I,EN) = SR * YR - SI * YI
0306          HI(I,EN) = SR * YI + SI * YR
0307   400 CONTINUE
0308 C
0309       DO 410 I = LOW, IGH
0310          YR = ZR(I,EN)
0311          YI = ZI(I,EN)
0312          ZR(I,EN) = SR * YR - SI * YI
0313          ZI(I,EN) = SR * YI + SI * YR
0314   410 CONTINUE
0315 C
0316       GOTO 250
0317 C     .......... A ROOT FOUND ..........
0318   420 HR(EN,EN) = HR(EN,EN) + TR
0319       WR(EN) = HR(EN,EN)
0320       HI(EN,EN) = HI(EN,EN) + TI
0321       WI(EN) = HI(EN,EN)
0322       EN = ENM1
0323       GOTO 240
0324 C     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND
0325 C                VECTORS OF UPPER TRIANGULAR FORM ..........
0326   430 NORM = 0.0D0
0327 C
0328       DO 440 I = 1, N
0329 C
0330          DO 440 J = I, N
0331             TR = DABS(HR(I,J)) + DABS(HI(I,J))
0332             IF (TR .GT. NORM) NORM = TR
0333   440 CONTINUE
0334 C
0335       IF (N .EQ. 1 .OR. NORM .EQ. 0.0D0) GOTO 560
0336 C     .......... FOR EN=N STEP -1 UNTIL 2 DO -- ..........
0337       DO 500 NN = 2, N
0338          EN = N + 2 - NN
0339          XR = WR(EN)
0340          XI = WI(EN)
0341          HR(EN,EN) = 1.0D0
0342          HI(EN,EN) = 0.0D0
0343          ENM1 = EN - 1
0344 C     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
0345          DO 490 II = 1, ENM1
0346             I = EN - II
0347             ZZR = 0.0D0
0348             ZZI = 0.0D0
0349             IP1 = I + 1
0350 C
0351             DO 450 J = IP1, EN
0352                ZZR = ZZR + HR(I,J) * HR(J,EN) - HI(I,J) * HI(J,EN)
0353                ZZI = ZZI + HR(I,J) * HI(J,EN) + HI(I,J) * HR(J,EN)
0354   450       CONTINUE
0355 C
0356             YR = XR - WR(I)
0357             YI = XI - WI(I)
0358             IF (YR .NE. 0.0D0 .OR. YI .NE. 0.0D0) GOTO 470
0359                TST1 = NORM
0360                YR = TST1
0361   460          YR = 0.01D0 * YR
0362                TST2 = NORM + YR
0363                IF (TST2 .GT. TST1) GOTO 460
0364   470       CONTINUE
0365             CALL PYCDIV(ZZR,ZZI,YR,YI,HR(I,EN),HI(I,EN))
0366 C     .......... OVERFLOW CONTROL ..........
0367             TR = DABS(HR(I,EN)) + DABS(HI(I,EN))
0368             IF (TR .EQ. 0.0D0) GOTO 490
0369             TST1 = TR
0370             TST2 = TST1 + 1.0D0/TST1
0371             IF (TST2 .GT. TST1) GOTO 490
0372             DO 480 J = I, EN
0373                HR(J,EN) = HR(J,EN)/TR
0374                HI(J,EN) = HI(J,EN)/TR
0375   480       CONTINUE
0376 C
0377   490    CONTINUE
0378 C
0379   500 CONTINUE
0380 C     .......... END BACKSUBSTITUTION ..........
0381 C     .......... VECTORS OF ISOLATED ROOTS ..........
0382       DO 520 I = 1, N
0383          IF (I .GE. LOW .AND. I .LE. IGH) GOTO 520
0384 C
0385          DO 510 J = I, N
0386             ZR(I,J) = HR(I,J)
0387             ZI(I,J) = HI(I,J)
0388   510    CONTINUE
0389 C
0390   520 CONTINUE
0391 C     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
0392 C                VECTORS OF ORIGINAL FULL MATRIX.
0393 C                FOR J=N STEP -1 UNTIL LOW DO -- ..........
0394       DO 540 JJ = LOW, N
0395          J = N + LOW - JJ
0396          M = MIN0(J,IGH)
0397 C
0398          DO 540 I = LOW, IGH
0399             ZZR = 0.0D0
0400             ZZI = 0.0D0
0401 C
0402             DO 530 K = LOW, M
0403                ZZR = ZZR + ZR(I,K) * HR(K,J) - ZI(I,K) * HI(K,J)
0404                ZZI = ZZI + ZR(I,K) * HI(K,J) + ZI(I,K) * HR(K,J)
0405   530       CONTINUE
0406 C
0407             ZR(I,J) = ZZR
0408             ZI(I,J) = ZZI
0409   540 CONTINUE
0410 C
0411       GOTO 560
0412 C     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
0413 C                CONVERGED AFTER 30*N ITERATIONS ..........
0414   550 IERR = EN
0415   560 RETURN
0416       END