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 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 this subroutine balances a real 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 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 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 c ------------------------------------------------------------------ c radix = 16.0d0 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.0d0) 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.0d0) 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.0d0 c .......... iterative loop for norm reduction .......... 190 noconv = .false. c do 270 i = k, l c = 0.0d0 r = 0.0d0 c do 200 j = k, l if (j .eq. i) go to 200 c = c + dabs(a(j,i)) r = r + dabs(a(i,j)) 200 continue c .......... guard against zero c or r due to underflow .......... if (c .eq. 0.0d0 .or. r .eq. 0.0d0) go to 270 g = r / radix f = 1.0d0 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.95d0 * s) go to 270 g = 1.0d0 / 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 end