dc/d72/pycbal_8f_source.html

C*********************************************************************


C...PYCBAL

C...Auxiliary to PYEICG

C

C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE

C     CBALANCE, WHICH IS A COMPLEX VERSION OF BALANCE,

C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.

C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).

C

C     THIS SUBROUTINE BALANCES A COMPLEX MATRIX AND ISOLATES

C     EIGENVALUES WHENEVER POSSIBLE.

C

C     ON INPUT

C

C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL

C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM

C          DIMENSION STATEMENT.

C

C        N IS THE ORDER OF THE MATRIX.

C

C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,

C          RESPECTIVELY, OF THE COMPLEX MATRIX TO BE BALANCED.

C

C     ON OUTPUT

C

C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,

C          RESPECTIVELY, OF THE BALANCED MATRIX.

C

C        LOW AND IGH ARE TWO INTEGERS SUCH THAT AR(I,J) AND AI(I,J)

C          ARE EQUAL TO ZERO IF

C           (1) I IS GREATER THAN J AND

C           (2) J=1,...,LOW-1 OR I=IGH+1,...,N.

C

C        SCALE CONTAINS INFORMATION DETERMINING THE

C           PERMUTATIONS AND SCALING FACTORS USED.

C

C     SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH

C     HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED

C     WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS

C     OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J).  THEN

C        SCALE(J) = P(J),    FOR J = 1,...,LOW-1

C                 = D(J,J)       J = LOW,...,IGH

C                 = P(J)         J = IGH+1,...,N.

C     THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1,

C     THEN 1 TO LOW-1.

C

C     NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY.

C

C     THE ALGOL PROCEDURE EXC CONTAINED IN CBALANCE APPEARS IN

C     CBAL  IN LINE.  (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS

C     K,L HAVE BEEN REVERSED.)

C

C     ARITHMETIC IS REAL THROUGHOUT.

C

C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,

C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY

C

C     THIS VERSION DATED AUGUST 1983.

C


      SUBROUTINE pycbal(NM,N,AR,AI,LOW,IGH,SCALE)


      INTEGER i,j,k,l,m,n,jj,nm,igh,low,iexc

      DOUBLE PRECISION ar(4,4),ai(4,4),scale(4)

      DOUBLE PRECISION c,f,g,r,s,b2,radix

      LOGICAL noconv


      radix = 16.0d0

C

      b2 = radix * radix

      k = 1

      l = n

      goto 150

C     .......... IN-LINE PROCEDURE FOR ROW AND

C                COLUMN EXCHANGE ..........

  100 scale(m) = j

      IF (j .EQ. m) goto 130

C

      DO 110 i = 1, l

         f = ar(i,j)

         ar(i,j) = ar(i,m)

         ar(i,m) = f

         f = ai(i,j)

         ai(i,j) = ai(i,m)

         ai(i,m) = f

  110 CONTINUE

C

      DO 120 i = k, n

         f = ar(j,i)

         ar(j,i) = ar(m,i)

         ar(m,i) = f

         f = ai(j,i)

         ai(j,i) = ai(m,i)

         ai(m,i) = f

  120 CONTINUE

C

  130 IF(iexc.EQ.1) goto 140

      IF(iexc.EQ.2) goto 180

C     .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE

C                AND PUSH THEM DOWN ..........

  140 IF (l .EQ. 1) goto 320

      l = l - 1

C     .......... FOR J=L STEP -1 UNTIL 1 DO -- ..........

  150 DO 170 jj = 1, l

         j = l + 1 - jj

C

         DO 160 i = 1, l

            IF (i .EQ. j) goto 160

            IF (ar(j,i) .NE. 0.0d0 .OR. ai(j,i) .NE. 0.0d0) goto 170

  160    CONTINUE

C

         m = l

         iexc = 1

         goto 100

  170 CONTINUE

C

      goto 190

C     .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE

C                AND PUSH THEM LEFT ..........

  180 k = k + 1

C

  190 DO 210 j = k, l

C

         DO 200 i = k, l

            IF (i .EQ. j) goto 200

            IF (ar(i,j) .NE. 0.0d0 .OR. ai(i,j) .NE. 0.0d0) goto 210

  200    CONTINUE

C

         m = k

         iexc = 2

         goto 100

  210 CONTINUE

C     .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L ..........

      DO 220 i = k, l

  220 scale(i) = 1.0d0

C     .......... ITERATIVE LOOP FOR NORM REDUCTION ..........

  230 noconv = .false.

C

      DO 310 i = k, l

         c = 0.0d0

         r = 0.0d0

C

         DO 240 j = k, l

            IF (j .EQ. i) goto 240

            c = c + dabs(ar(j,i)) + dabs(ai(j,i))

            r = r + dabs(ar(i,j)) + dabs(ai(i,j))

  240    CONTINUE

C     .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW ..........

         IF (c .EQ. 0.0d0 .OR. r .EQ. 0.0d0) goto 310

         g = r / radix

         f = 1.0d0

         s = c + r

  250    IF (c .GE. g) goto 260

         f = f * radix

         c = c * b2

         goto 250

  260    g = r * radix

  270    IF (c .LT. g) goto 280

         f = f / radix

         c = c / b2

         goto 270

C     .......... NOW BALANCE ..........

  280    IF ((c + r) / f .GE. 0.95d0 * s) goto 310

         g = 1.0d0 / f

         scale(i) = scale(i) * f

         noconv = .true.

C

         DO 290 j = k, n

            ar(i,j) = ar(i,j) * g

            ai(i,j) = ai(i,j) * g

  290    CONTINUE

C

         DO 300 j = 1, l

            ar(j,i) = ar(j,i) * f

            ai(j,i) = ai(j,i) * f

  300    CONTINUE

C

  310 CONTINUE

C

      IF (noconv) goto 230

C

  320 low = k

      igh = l

      RETURN

      END