pythiaeRHIC/pythia/pycbal.f

0001
0002 C*********************************************************************
0003
0004 C...PYCBAL
0005 C...Auxiliary to PYEICG
0006 C
0007 C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE
0008 C     CBALANCE, WHICH IS A COMPLEX VERSION OF BALANCE,
0009 C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
0010 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
0011 C
0012 C     THIS SUBROUTINE BALANCES A COMPLEX MATRIX AND ISOLATES
0013 C     EIGENVALUES WHENEVER POSSIBLE.
0014 C
0015 C     ON INPUT
0016 C
0017 C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
0018 C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
0019 C          DIMENSION STATEMENT.
0020 C
0021 C        N IS THE ORDER OF THE MATRIX.
0022 C
0023 C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
0024 C          RESPECTIVELY, OF THE COMPLEX MATRIX TO BE BALANCED.
0025 C
0026 C     ON OUTPUT
0027 C
0028 C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
0029 C          RESPECTIVELY, OF THE BALANCED MATRIX.
0030 C
0031 C        LOW AND IGH ARE TWO INTEGERS SUCH THAT AR(I,J) AND AI(I,J)
0032 C          ARE EQUAL TO ZERO IF
0033 C           (1) I IS GREATER THAN J AND
0034 C           (2) J=1,...,LOW-1 OR I=IGH+1,...,N.
0035 C
0036 C        SCALE CONTAINS INFORMATION DETERMINING THE
0037 C           PERMUTATIONS AND SCALING FACTORS USED.
0038 C
0039 C     SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH
0040 C     HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED
0041 C     WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS
0042 C     OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J).  THEN
0043 C        SCALE(J) = P(J),    FOR J = 1,...,LOW-1
0044 C                 = D(J,J)       J = LOW,...,IGH
0045 C                 = P(J)         J = IGH+1,...,N.
0046 C     THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1,
0047 C     THEN 1 TO LOW-1.
0048 C
0049 C     NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY.
0050 C
0051 C     THE ALGOL PROCEDURE EXC CONTAINED IN CBALANCE APPEARS IN
0052 C     CBAL  IN LINE.  (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS
0053 C     K,L HAVE BEEN REVERSED.)
0054 C
0055 C     ARITHMETIC IS REAL THROUGHOUT.
0056 C
0057 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
0058 C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
0059 C
0060 C     THIS VERSION DATED AUGUST 1983.
0061 C
0062
0063       SUBROUTINE PYCBAL(NM,N,AR,AI,LOW,IGH,SCALE)
0064
0065       INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC
0066       DOUBLE PRECISION AR(4,4),AI(4,4),SCALE(4)
0067       DOUBLE PRECISION C,F,G,R,S,B2,RADIX
0068       LOGICAL NOCONV
0069
0070       RADIX = 16.0D0
0071 C
0072       B2 = RADIX * RADIX
0073       K = 1
0074       L = N
0075       GOTO 150
0076 C     .......... IN-LINE PROCEDURE FOR ROW AND
0077 C                COLUMN EXCHANGE ..........
0078   100 SCALE(M) = J
0079       IF (J .EQ. M) GOTO 130
0080 C
0081       DO 110 I = 1, L
0082          F = AR(I,J)
0083          AR(I,J) = AR(I,M)
0084          AR(I,M) = F
0085          F = AI(I,J)
0086          AI(I,J) = AI(I,M)
0087          AI(I,M) = F
0088   110 CONTINUE
0089 C
0090       DO 120 I = K, N
0091          F = AR(J,I)
0092          AR(J,I) = AR(M,I)
0093          AR(M,I) = F
0094          F = AI(J,I)
0095          AI(J,I) = AI(M,I)
0096          AI(M,I) = F
0097   120 CONTINUE
0098 C
0099   130 IF(IEXC.EQ.1) GOTO 140
0100       IF(IEXC.EQ.2) GOTO 180
0101 C     .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE
0102 C                AND PUSH THEM DOWN ..........
0103   140 IF (L .EQ. 1) GOTO 320
0104       L = L - 1
0105 C     .......... FOR J=L STEP -1 UNTIL 1 DO -- ..........
0106   150 DO 170 JJ = 1, L
0107          J = L + 1 - JJ
0108 C
0109          DO 160 I = 1, L
0110             IF (I .EQ. J) GOTO 160
0111             IF (AR(J,I) .NE. 0.0D0 .OR. AI(J,I) .NE. 0.0D0) GOTO 170
0112   160    CONTINUE
0113 C
0114          M = L
0115          IEXC = 1
0116          GOTO 100
0117   170 CONTINUE
0118 C
0119       GOTO 190
0120 C     .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE
0121 C                AND PUSH THEM LEFT ..........
0122   180 K = K + 1
0123 C
0124   190 DO 210 J = K, L
0125 C
0126          DO 200 I = K, L
0127             IF (I .EQ. J) GOTO 200
0128             IF (AR(I,J) .NE. 0.0D0 .OR. AI(I,J) .NE. 0.0D0) GOTO 210
0129   200    CONTINUE
0130 C
0131          M = K
0132          IEXC = 2
0133          GOTO 100
0134   210 CONTINUE
0135 C     .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L ..........
0136       DO 220 I = K, L
0137   220 SCALE(I) = 1.0D0
0138 C     .......... ITERATIVE LOOP FOR NORM REDUCTION ..........
0139   230 NOCONV = .FALSE.
0140 C
0141       DO 310 I = K, L
0142          C = 0.0D0
0143          R = 0.0D0
0144 C
0145          DO 240 J = K, L
0146             IF (J .EQ. I) GOTO 240
0147             C = C + DABS(AR(J,I)) + DABS(AI(J,I))
0148             R = R + DABS(AR(I,J)) + DABS(AI(I,J))
0149   240    CONTINUE
0150 C     .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW ..........
0151          IF (C .EQ. 0.0D0 .OR. R .EQ. 0.0D0) GOTO 310
0152          G = R / RADIX
0153          F = 1.0D0
0154          S = C + R
0155   250    IF (C .GE. G) GOTO 260
0156          F = F * RADIX
0157          C = C * B2
0158          GOTO 250
0159   260    G = R * RADIX
0160   270    IF (C .LT. G) GOTO 280
0161          F = F / RADIX
0162          C = C / B2
0163          GOTO 270
0164 C     .......... NOW BALANCE ..........
0165   280    IF ((C + R) / F .GE. 0.95D0 * S) GOTO 310
0166          G = 1.0D0 / F
0167          SCALE(I) = SCALE(I) * F
0168          NOCONV = .TRUE.
0169 C
0170          DO 290 J = K, N
0171             AR(I,J) = AR(I,J) * G
0172             AI(I,J) = AI(I,J) * G
0173   290    CONTINUE
0174 C
0175          DO 300 J = 1, L
0176             AR(J,I) = AR(J,I) * F
0177             AI(J,I) = AI(J,I) * F
0178   300    CONTINUE
0179 C
0180   310 CONTINUE
0181 C
0182       IF (NOCONV) GOTO 230
0183 C
0184   320 LOW = K
0185       IGH = L
0186       RETURN
0187       END