Balance(M) = Return a balanced matrix with same eigenvalues

Balance(M) = Return a balanced matrix with same eigenvalues:

#include 
#define RADIX 2.0
void balance(float **a, int n)
// Given a matrix a[1..n][1..n], this routine replaces it by a balanced matrix
// with identical eigenvalues. A symmetric matrix is already balanced and is
// unaffected by this procedure. The parameter RADIX should be the machine's
// floating-point radix.
{
   int last, j, i;
   float s, r, g, f, c, sqrdx;
   sqrdx = RADIX * RADIX;
   last = 0;
   while (last == 0)
   {
      last = 1;
      for (i = 1; i <= n; i++) // Calculate row and column norms.
      {
         r = c= 0.0;
         for (j = 1; j <= n; j++)
            if (j != i)
            {
               c += fabs(a[j][i]);
               r += fabs(a[i][j]);
            }
         if (c && r) // If both are nonzero,
         {
            g = r / RADIX;
            f = 1.0;
            s = c + r;
            while (c < g) // find the integer power of the machine radix that
            {             // comes closest to balancing the matrix.
               f *= RADIX;
               c *= sqrdx;
            }
            g = r * RADIX;
            while (c > g)
            {
               f /= RADIX;
               c /= sqrdx;
            }
            if ((c + r) / f < 0.95*s)
            {
               last = 0;
               g = 1.0 / f;
               for (j = 1; j <= n; j++) // Apply similarity transformation.
                  a[i][j] *= g;
               for (j = 1; j <= n; j++)
                  a[j][i] *= f;
            }
         }
      }
   }
} // balance

See: Reduction of a General Matrix to Hessenberg Form -- From Numerical Recipes in C

Return to Matrix and Polynomial Computations
Return to Harry's Home Page

---

This page accessed times since October 8, 2006.
Page created by: hjsmithh@sbcglobal.net
Changes last made on Monday, 06-Aug-07 20:22:39 PDT

---