pythiaeRHIC/pythia/pycrth.f

0001
0002 C*********************************************************************
0003
0004 C...PYCRTH
0005 C...Auxiliary to PYEICG.
0006 C
0007 C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
0008 C     THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968)
0009 C     BY MARTIN AND WILKINSON.
0010 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
0011 C
0012 C     GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE
0013 C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
0014 C     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
0015 C     UNITARY SIMILARITY TRANSFORMATIONS.
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        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
0030 C          RESPECTIVELY, OF THE COMPLEX INPUT MATRIX.
0031 C
0032 C     ON OUTPUT
0033 C
0034 C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
0035 C          RESPECTIVELY, OF THE HESSENBERG MATRIX.  INFORMATION
0036 C          ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION
0037 C          IS STORED IN THE REMAINING TRIANGLES UNDER THE
0038 C          HESSENBERG MATRIX.
0039 C
0040 C        ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE
0041 C          TRANSFORMATIONS.  ONLY ELEMENTS LOW THROUGH IGH ARE USED.
0042 C
0043 C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
0044 C
0045 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
0046 C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
0047 C
0048 C     THIS VERSION DATED AUGUST 1983.
0049 C
0050
0051       SUBROUTINE PYCRTH(NM,N,LOW,IGH,AR,AI,ORTR,ORTI)
0052
0053       INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW
0054       DOUBLE PRECISION AR(4,4),AI(4,4),ORTR(4),ORTI(4)
0055       DOUBLE PRECISION F,G,H,FI,FR,SCALE,PYTHAG
0056
0057       LA = IGH - 1
0058       KP1 = LOW + 1
0059       IF (LA .LT. KP1) GOTO 210
0060 C
0061       DO 200 M = KP1, LA
0062          H = 0.0D0
0063          ORTR(M) = 0.0D0
0064          ORTI(M) = 0.0D0
0065          SCALE = 0.0D0
0066 C     .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........
0067          DO 100 I = M, IGH
0068   100    SCALE = SCALE + DABS(AR(I,M-1)) + DABS(AI(I,M-1))
0069 C
0070          IF (SCALE .EQ. 0.0D0) GOTO 200
0071          MP = M + IGH
0072 C     .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
0073          DO 110 II = M, IGH
0074             I = MP - II
0075             ORTR(I) = AR(I,M-1) / SCALE
0076             ORTI(I) = AI(I,M-1) / SCALE
0077             H = H + ORTR(I) * ORTR(I) + ORTI(I) * ORTI(I)
0078   110    CONTINUE
0079 C
0080          G = DSQRT(H)
0081          F = PYTHAG(ORTR(M),ORTI(M))
0082          IF (F .EQ. 0.0D0) GOTO 120
0083          H = H + F * G
0084          G = G / F
0085          ORTR(M) = (1.0D0 + G) * ORTR(M)
0086          ORTI(M) = (1.0D0 + G) * ORTI(M)
0087          GOTO 130
0088 C
0089   120    ORTR(M) = G
0090          AR(M,M-1) = SCALE
0091 C     .......... FORM (I-(U*UT)/H) * A ..........
0092   130    DO 160 J = M, N
0093             FR = 0.0D0
0094             FI = 0.0D0
0095 C     .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
0096             DO 140 II = M, IGH
0097                I = MP - II
0098                FR = FR + ORTR(I) * AR(I,J) + ORTI(I) * AI(I,J)
0099                FI = FI + ORTR(I) * AI(I,J) - ORTI(I) * AR(I,J)
0100   140       CONTINUE
0101 C
0102             FR = FR / H
0103             FI = FI / H
0104 C
0105             DO 150 I = M, IGH
0106                AR(I,J) = AR(I,J) - FR * ORTR(I) + FI * ORTI(I)
0107                AI(I,J) = AI(I,J) - FR * ORTI(I) - FI * ORTR(I)
0108   150       CONTINUE
0109 C
0110   160    CONTINUE
0111 C     .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) ..........
0112          DO 190 I = 1, IGH
0113             FR = 0.0D0
0114             FI = 0.0D0
0115 C     .......... FOR J=IGH STEP -1 UNTIL M DO -- ..........
0116             DO 170 JJ = M, IGH
0117                J = MP - JJ
0118                FR = FR + ORTR(J) * AR(I,J) - ORTI(J) * AI(I,J)
0119                FI = FI + ORTR(J) * AI(I,J) + ORTI(J) * AR(I,J)
0120   170       CONTINUE
0121 C
0122             FR = FR / H
0123             FI = FI / H
0124 C
0125             DO 180 J = M, IGH
0126                AR(I,J) = AR(I,J) - FR * ORTR(J) - FI * ORTI(J)
0127                AI(I,J) = AI(I,J) + FR * ORTI(J) - FI * ORTR(J)
0128   180       CONTINUE
0129 C
0130   190    CONTINUE
0131 C
0132          ORTR(M) = SCALE * ORTR(M)
0133          ORTI(M) = SCALE * ORTI(M)
0134          AR(M,M-1) = -G * AR(M,M-1)
0135          AI(M,M-1) = -G * AI(M,M-1)
0136   200 CONTINUE
0137 C
0138   210 RETURN
0139       END