      subroutine balanc(nm,n,a,low,igh,scale)
c
      integer i,j,k,l,m,n,jj,nm,igh,low,iexc
      double precision a(nm,n),scale(n)
      double precision c,f,g,r,s,b2,radix
      logical noconv
c! purpose
c
c     this subroutine balances a real matrix and isolates
c     eigenvalues whenever possible.
c! calling sequence
c
c      subroutine balanc(nm,n,a,low,igh,scale)
c
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        a contains the input matrix to be balanced.
c
c     on output:
c
c        a contains the balanced matrix;
c
c        low and igh are two integers such that a(i,j)
c          is 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 balance appears in
c     balanc  in line.  (note that the algol roles of identifiers
c     k,l have been reversed.)
c
c! originator
c
c     this subroutine is a translation of the algol procedure 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     questions and comments should be directed to b. s. garbow,
c     applied mathematics division, argonne national laboratory
c
c     ------------------------------------------------------------------
c
c     :::::::::: radix is a machine dependent parameter specifying
c                the base of the machine floating point representation.
      data radix/2.0d+0/
c
      b2 = radix * radix
      k = 1
      l = n
      go to 100
c     :::::::::: in-line procedure for row and
c                column exchange ::::::::::
   20 scale(m) = j
      if (j .eq. m) go to 50
c
      do 30 i = 1, l
         f = a(i,j)
         a(i,j) = a(i,m)
         a(i,m) = f
   30 continue
c
      do 40 i = k, n
         f = a(j,i)
         a(j,i) = a(m,i)
         a(m,i) = f
   40 continue
c
   50 go to (80,130), iexc
c     :::::::::: search for rows isolating an eigenvalue
c                and push them down ::::::::::
   80 if (l .eq. 1) go to 280
      l = l - 1
c     :::::::::: for j=l step -1 until 1 do -- ::::::::::
  100 do 120 jj = 1, l
         j = l + 1 - jj
c
         do 110 i = 1, l
            if (i .eq. j) go to 110
            if (a(j,i) .ne. 0.0d+0) go to 120
  110    continue
c
         m = l
         iexc = 1
         go to 20
  120 continue
c
      go to 140
c     :::::::::: search for columns isolating an eigenvalue
c                and push them left ::::::::::
  130 k = k + 1
c
  140 do 170 j = k, l
c
         do 150 i = k, l
            if (i .eq. j) go to 150
            if (a(i,j) .ne. 0.0d+0) go to 170
  150    continue
c
         m = k
         iexc = 2
         go to 20
  170 continue
c     :::::::::: now balance the submatrix in rows k to l ::::::::::
      do 180 i = k, l
  180 scale(i) = 1.0d+0
c     :::::::::: iterative loop for norm reduction ::::::::::
  190 noconv = .false.
c
      do 270 i = k, l
         c = 0.0d+0
         r = 0.0d+0
c
         do 200 j = k, l
            if (j .eq. i) go to 200
            c = c + abs(a(j,i))
            r = r + abs(a(i,j))
  200    continue
c     :::::::::: guard against zero c or r due to underflow ::::::::::
         if (c .eq. 0.0d+0 .or. r .eq. 0.0d+0) go to 270
         g = r / radix
         f = 1.0d+0
         s = c + r
  210    if (c .ge. g) go to 220
         f = f * radix
         c = c * b2
         go to 210
  220    g = r * radix
  230    if (c .lt. g) go to 240
         f = f / radix
         c = c / b2
         go to 230
c     :::::::::: now balance ::::::::::
  240    if ((c + r) / f .ge. 0.950d+0 * s) go to 270
         g = 1.0d+0 / f
         scale(i) = scale(i) * f
         noconv = .true.
c
         do 250 j = k, n
  250    a(i,j) = a(i,j) * g
c
         do 260 j = 1, l
  260    a(j,i) = a(j,i) * f
c
  270 continue
c
      if (noconv) go to 190
c
  280 low = k
      igh = l
      return
c     :::::::::: last card of balanc ::::::::::
      end
