Hess(M) = Return upper Hessenberg matrix, same eigenvalues

Hess(M) = Return upper Hessenberg matrix, same eigenvalues:

#include 
#define SWAP(g, h) {y = (g); (g) = (h); (h) = y;}
void hess(float **a, int n)
// Reduction to Hessenberg form by the elimination method. The real,
// nonsymmetric matrix a[1..n][1..n] is replaced by an upper Hessenberg matrix
// with identical eigenvalues. Recommended, but not required, is that this
// routine be preceded by balance. On output, the Hessenberg matrix is in
// elements a[i][j] with i <= j + 1. Elements with i > j+1 are zero.
{
   int m, j, i;
   float y, x;
   for (m = 2; m < n; m++) // m is called r + 1 in the text.
   {
      X = 0.0;
      i = m;
      for (j = m; j <= n; j++) // Find the pivot.
      {
         if (fabs(a[j][m - 1]) > fabs(x))
         {
            x = a[j][m - 1];
            i = j;
         }
      }
      if (i != m) // Interchange rows and columns.
      {
         for (j = m - 1; j <= n; j++)
            SWAP(a[i][j], a[m][j])
         for (j = 1; j <= n; ++)
            SWAP(a[j][i], a[j][m])
      }
      if (x != 0.0) // Carry out the elimination.
      {
         for (I = m + 1; I <= n; i++)
         {
            if ((y = a[i][m - 1]) != 0.0)
            {
               y /= x;
               a[i][m-1] = y;
               for (j = m; j <= n; j++)
                  a[i][j] -= y * a[m][j];
               for (j = 1; j <= n; j++)
                  a[j][m] += y * a[j][i];
            }
         }
      }
   }
   for (j = 1; j <= n - 2; j++)
      for (i = j + 2; i <= n; i++)
         a[i, j] = 0.0;
} // hess

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:45 PDT
1