feigen.cpp

Go to the documentation of this file.

//http://perso.orange.fr/jean-pierre.moreau/c_matrices.html
/* ----------------------- MODULE feigen.cpp ------------------------ *
 * Reference: "Numerical Algorithms with C By G. Engeln-Mueller and   *
 *             F. Uhlig, Springer-Verlag, 1996" [BIBLI 11].           *
 * ------------------------------------------------------------------ */
#include <basis.h>
#include <vmblock.h>


#define MAXIT 50                      /*  Maximal number of           */
                                      /*  iterations per eigenvalue   */

/*--------------------------------------------------------------------*
 * Auxiliary functions for  eigen                                     *
 *--------------------------------------------------------------------*/


static int balance       /* balance a matrix .........................*/
                   (int       n,      /* size of matrix ..............*/
                    REAL *    mat[],  /* matrix ......................*/
                    REAL      scal[], /* Scaling data ................*/
                    int *     low,    /* first relevant row index ....*/
                    int *     high,   /* last relevant row index .....*/
                    int       basis   /* base of computer numbers ....*/
                   )
/*====================================================================*
 *                                                                    *
 *  balance balances the matrix mat so that the rows with zero entries*
 *  off the diagonal are isolated and the remaining columns and rows  *
 *  are resized to have one norm close to 1.                          *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;  ( n > 0 )                                    *
 *               Dimension of mat                                     *
 *      mat      REAL   *mat[n];                                      *
 *               n x n input matrix                                   *
 *      basis    int basis;                                           *
 *               Base of number representaion in the given computer   *
 *               (see BASIS)                                          *
 *                                                                    *
 *   Output parameters:                                               *
 *   ==================                                               *
 *      mat      REAL   *mat[n];                                      *
 *               scaled matrix                                        *
 *      low      int *low;                                            *
 *      high     int *high;                                           *
 *               the rows 0 to low-1 and those from high to n-1       *
 *               contain isolated eigenvalues (only nonzero entry on  *
 *               the diagonal)                                        *
 *      scal     REAL   scal[];                                       *
 *               the vector scal contains the isolated eigenvalues in *
 *               the positions 0 to low-1 and high to n-1, its other  *
 *               components contain the scaling factors for           *
 *               transforming mat.                                    *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Macros:     SWAP, ABS                                            *
 *   =======                                                          *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Constants used:     TRUE, FALSE                                  *
 *   ==============                                                   *
 *                                                                    *
 *====================================================================*/
{
  register int i, j;
  int      iter, k, m;
  REAL     b2, r, c, f, g, s;

  b2 = (REAL) (basis * basis);
  m = 0;
  k = n - 1;

  do
  {
    iter = FALSE;
    for (j = k; j >= 0; j--)
    {
      for (r = ZERO, i = 0; i <= k; i++)
        if (i != j)  r += ABS (mat[j][i]);

      if (r == ZERO)
      {
        scal[k] = (REAL) j;
        if (j != k)
        {
          for (i = 0; i <= k; i++) SWAP (REAL, mat[i][j], mat[i][k])
          for (i = m; i < n; i++)  SWAP (REAL, mat[j][i], mat[k][i])
        }
        k--;
        iter = TRUE;
      }
    }   /* end of j */
  }   /* end of do  */
  while (iter);

  do
  {
    iter = FALSE;
    for (j = m; j <= k; j++)
    {
      for (c = ZERO, i = m; i <= k; i++)
        if (i != j) c += ABS (mat[i][j]);
      if (c == ZERO)
      {
        scal[m] = (REAL) j;
        if ( j != m )
        {
          for (i = 0; i <= k; i++) SWAP (REAL, mat[i][j], mat[i][m])
          for (i = m; i < n; i++)  SWAP (REAL, mat[j][i], mat[m][i])
        }
        m++;
        iter = TRUE;
      }
    }   /* end of j */
  }   /* end of do  */
  while (iter);

  *low = m;
  *high = k;
  for (i = m; i <= k; i++) scal[i] = ONE;

  do
  {
    iter = FALSE;
    for (i = m; i <= k; i++)
    {
      for (c = r = ZERO, j = m; j <= k; j++)
      if (j !=i)
      {
        c += ABS (mat[j][i]);
        r += ABS (mat[i][j]);
      }
      g = r / basis;
      f = ONE;
      s = c + r;

      while (c < g)
      {
        f *= basis;
        c *= b2;
      }

      g = r * basis;
      while (c >= g)
      {
        f /= basis;
        c /= b2;
      }

      if ((c + r) / f < (REAL)0.95 * s)
      {
        g = ONE / f;
        scal[i] *= f;
        iter = TRUE;
        for (j = m; j < n; j++ ) mat[i][j] *= g;
        for (j = 0; j <= k; j++ ) mat[j][i] *= f;
      }
    }
  }
  while (iter);

  return (0);
}


static int balback       /* reverse balancing ........................*/
                   (int     n,        /* Dimension of matrix .........*/
                    int     low,      /* first nonzero row ...........*/
                    int     high,     /* last nonzero row ............*/
                    REAL    scal[],   /* Scaling data ................*/
                    REAL *  eivec[]   /* Eigenvectors ................*/
                   )
/*====================================================================*
 *                                                                    *
 *  balback reverses the balancing of balance for the eigenvactors.   *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;  ( n > 0 )                                    *
 *               Dimension of mat                                     *
 *      low      int low;                                             *
 *      high     int high;   see balance                              *
 *      eivec    REAL   *eivec[n];                                    *
 *               Matrix of eigenvectors, as computed in  qr2          *
 *      scal     REAL   scal[];                                       *
 *               Scaling data from  balance                           *
 *                                                                    *
 *   Output parameter:                                                *
 *   ================                                                 *
 *      eivec    REAL   *eivec[n];                                    *
 *               Non-normalized eigenvectors of the original matrix   *
 *                                                                    *
 *   Macros :    SWAP()                                               *
 *   ========                                                         *
 *                                                                    *
 *====================================================================*/
{
  register int i, j, k;
  REAL s;

  for (i = low; i <= high; i++)
  {
    s = scal[i];
    for (j = 0; j < n; j++) eivec[i][j] *= s;
  }

  for (i = low - 1; i >= 0; i--)
  {
    k = (int) scal[i];
    if (k != i)
      for (j = 0; j < n; j++) SWAP (REAL, eivec[i][j], eivec[k][j])
  }

  for (i = high + 1; i < n; i++)
  {
    k = (int) scal[i];
    if (k != i)
      for (j = 0; j < n; j++) SWAP (REAL, eivec[i][j], eivec[k][j])
  }
  return (0);
}



static int elmhes       /* reduce matrix to upper Hessenberg form ....*/
                  (int       n,       /* Dimension of matrix .........*/
                   int       low,     /* first nonzero row ...........*/
                   int       high,    /* last nonzero row ............*/
                   REAL *    mat[],   /* input/output matrix .........*/
                   int       perm[]   /* Permutation vector ..........*/
                  )
/*====================================================================*
 *                                                                    *
 *  elmhes transforms the matrix mat to upper Hessenberg form.        *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;  ( n > 0 )                                    *
 *               Dimension of mat                                     *
 *      low      int low;                                             *
 *      high     int high; see  balance                               *
 *      mat      REAL   *mat[n];                                      *
 *               n x n matrix                                         *
 *                                                                    *
 *   Output parameter:                                                *
 *   =================                                                *
 *      mat      REAL   *mat[n];                                      *
 *               upper Hessenberg matrix; additional information on   *
 *               the transformation is stored in the lower triangle   *
 *      perm     int perm[];                                          *
 *               Permutation vector for elmtrans                      *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Macros:   SWAP, ABS                                              *
 *   =======                                                          *
 *                                                                    *
 *====================================================================*/
{
  register int i, j, m;
  REAL   x, y;

  for (m = low + 1; m < high; m++)
  {
    i = m;
    x = ZERO;
    for (j = m; j <= high; j++)
      if (ABS (mat[j][m-1]) > ABS (x))
      {
        x = mat[j][m-1];
        i = j;
      }

    perm[m] = i;
    if (i != m)
    {
      for (j = m - 1; j < n; j++) SWAP (REAL, mat[i][j], mat[m][j])
      for (j = 0; j <= high; j++) SWAP (REAL, mat[j][i], mat[j][m])
    }

    if (x != ZERO)
    {
      for (i = m + 1; i <= high; i++)
      {
        y = mat[i][m-1];
        if (y != ZERO)
        {
          y = mat[i][m-1] = y / x;
          for (j = m; j < n; j++) mat[i][j] -= y * mat[m][j];
          for (j = 0; j <= high; j++) mat[j][m] += y * mat[j][i];
        }
      } /* end i */
    }
  } /* end m */

  return (0);
}


static int elmtrans       /* copy to Hessenberg form .................*/
                    (int     n,       /* Dimension of matrix .........*/
                     int     low,     /* first nonzero row ...........*/
                     int     high,    /* last nonzero row ............*/
                     REAL *  mat[],   /* input matrix ................*/
                     int     perm[],  /* row permutations ............*/
                     REAL *  h[]      /* Hessenberg matrix ...........*/
                    )
/*====================================================================*
 *                                                                    *
 *  elmtrans copies the Hessenberg matrix stored in mat to h.         *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;  ( n > 0 )                                    *
 *               Dimension of  mat and eivec                          *
 *      low      int low;                                             *
 *      high     int high; see  balance                               *
 *      mat      REAL   *mat[n];                                      *
 *               n x n input matrix                                   *
 *      perm     int *perm;                                           *
 *               Permutation data from  elmhes                        *
 *                                                                    *
 *   Output parameter:                                                *
 *   ================                                                 *
 *      h        REAL   *h[n];                                        *
 *               Hessenberg matrix                                    *
 *                                                                    *
 *====================================================================*/
{
  register int k, i, j;

  for (i = 0; i < n; i++)
  {
    for (k = 0; k < n; k++) h[i][k] = ZERO;
    h[i][i] = ONE;
  }

  for (i = high - 1; i > low; i--)
  {
    j = perm[i];
    for (k = i + 1; k <= high; k++) h[k][i] = mat[k][i-1];
    if ( i != j )
    {
      for (k = i; k <= high; k++)
      {
        h[i][k] = h[j][k];
        h[j][k] = ZERO;
      }
      h[j][i] = ONE;
    }
  }

  return (0);
}


/* ------------------------------------------------------------------ */

static int orthes     /* reduce orthogonally to upper Hessenberg form */
                 (
                  int  n,                  /* Dimension of matrix     */
                  int  low,                /* [low,low]..[high,high]: */
                  int  high,               /* submatrix to be reduced */
                  REAL *mat[],             /* input/output matrix     */
                  REAL d[]                 /* reduction information   */
                 )                         /* error code              */

/***********************************************************************
* This function reduces matrix mat to upper Hessenberg form by         *
* Householder transformations. All details of the transformations are  *
* stored in the remaining triangle of the Hessenberg matrix and in     *
* vector d.                                                            *
*                                                                      *
* Input parameters:                                                    *
* =================                                                    *
* n        dimension of mat                                            *
* low  \   rows 0 to low-1 and high+1 to n-1 contain isolated          *
* high  >  eigenvalues, i. e. eigenvalues corresponding to             *
*      /   eigenvectors that are multiples of unit vectors             *
* mat      [0..n-1,0..n-1] matrix to be reduced                        *
*                                                                      *
* Output parameters:                                                   *
* ==================                                                   *
* mat      the desired Hessenberg matrix together with the first part  *
*          of the reduction information below the subdiagonal          *
* d        [0..n-1] vector with the remaining reduction information    *
*                                                                      *
* Return value:                                                        *
* =============                                                        *
* Error code. This can only be the value 0 here.                       *
*                                                                      *
* global names used:                                                   *
* ==================                                                   *
* REAL, MACH_EPS, ZERO, SQRT                                           *
*                                                                      *
* -------------------------------------------------------------------- *
* Literature: Numerical Mathematics 12 (1968), pages 359 and 360       *
***********************************************************************/

{
  int  i, j, m;    /* loop variables                                  */
  REAL s,          /* Euclidian norm sigma of the subdiagonal column  */
                   /* vector v of mat, that shall be reflected into a */
                   /* multiple of the unit vector e1 = (1,0,...,0)    */
                   /* (v = (v1,..,v(high-m+1))                        */
       x = ZERO,   /* first element of v in the beginning, then       */
                   /* summation variable in the actual Householder    */
                   /* transformation                                  */
       y,          /* sigma^2 in the beginning, then ||u||^2, with    */
                   /* u := v +- sigma * e1                            */
       eps;        /* tolerance for checking if the transformation is */
                   /* valid                                           */

  eps = (REAL)128.0 * MACH_EPS;

  for (m = low + 1; m < high; m++)
  {
    for (y = ZERO, i = high; i >= m; i--)
      x    = mat[i][m - 1],
      d[i] = x,
      y    = y + x * x;
    if (y <= eps)
      s = ZERO;
    else
    {
      s = (x >= ZERO) ? -SQRT(y) : SQRT(y);
      y    -= x * s;
      d[m] =  x - s;

      for (j = m; j < n; j++)               /* multiply mat from the  */
      {                                     /* left by  (E-(u*uT)/y)  */
        for (x = ZERO, i = high; i >= m; i--)
          x += d[i] * mat[i][j];
        for (x /= y, i = m; i <= high; i++)
          mat[i][j] -= x * d[i];
      }

      for (i = 0; i <= high; i++)           /* multiply mat from the  */
      {                                     /* right by  (E-(u*uT)/y) */
        for (x = ZERO, j = high; j >= m; j--)
          x += d[j] * mat[i][j];
        for (x /= y, j = m; j <= high; j++)
          mat[i][j] -= x * d[j];
      }
    }

    mat[m][m - 1] = s;
  }

  return 0;

}    /* --------------------------- orthes -------------------------- */


/* ------------------------------------------------------------------ */

static int orttrans       /* compute orthogonal transformation matrix */
                   (
                    int  n,      /* Dimension of matrix               */
                    int  low,    /* [low,low]..[high,high]: submatrix */
                    int  high,   /* affected by the reduction         */
                    REAL *mat[], /* Hessenberg matrix, reduction inf. */
                    REAL d[],    /* remaining reduction information   */
                    REAL *v[]    /* transformation matrix             */
                   )             /* error code                        */

/***********************************************************************
* compute the matrix v of accumulated transformations from the         *
* information left by the Householder reduction of matrix mat to upper *
* Hessenberg form below the Hessenberg matrix in mat and in the        *
* vector d. The contents of the latter are destroyed.                  *
*                                                                      *
* Input parameters:                                                    *
* =================                                                    *
* n        dimension of mat                                            *
* low  \   rows 0 to low-1 and high+1 to n-1 contain isolated          *
* high  >  eigenvalues, i. e. eigenvalues corresponding to             *
*      /   eigenvectors that are multiples of unit vectors             *
* mat      [0..n-1,0..n-1] matrix produced by `orthes' giving the      *
*          upper Hessenberg matrix and part of the information on the  *
*          orthogonal reduction                                        *
* d        [0..n-1] vector with the remaining information on the       *
*          orthogonal reduction to upper Hessenberg form               *
*                                                                      *
* Output parameters:                                                   *
* ==================                                                   *
* d        input vector destroyed by this function                     *
* v        [0..n-1,0..n-1] matrix defining the similarity reduction    *
*          to upper Hessenberg form                                    *
*                                                                      *
* Return value:                                                        *
* =============                                                        *
* Error code. This can only be the value 0 here.                       *
*                                                                      *
* global names used:                                                   *
* =================                                                    *
* REAL, ZERO, ONE                                                      *
*                                                                      *
* -------------------------------------------------------------------- *
* Literature: Numerical Mathematics 16 (1970), page 191                *
***********************************************************************/

{
  int  i, j, m;                        /* loop variables              */
  REAL x,                              /* summation variable in the   */
                                       /* Householder transformation  */
       y;                              /* sigma  respectively         */
                                       /* sigma * (v1 +- sigma)       */

  for (i = 0; i < n; i++)              /* form the unit matrix in v   */
  {
    for (j = 0; j < n; j++)
      v[i][j] = ZERO;
    v[i][i] = ONE;
  }

  for (m = high - 1; m > low; m--)     /* apply the transformations   */
  {                                    /* that reduced mat to upper   */
    y = mat[m][m - 1];                 /* Hessenberg form also to the */
                                       /* unit matrix in v. This      */
    if (y != ZERO)                     /* produces the desired        */
    {                                  /* transformation matrix in v. */
      y *= d[m];
      for (i = m + 1; i <= high; i++)
        d[i] = mat[i][m - 1];
      for (j = m; j <= high; j++)
      {
        for (x = ZERO, i = m; i <= high; i++)
          x += d[i] * v[i][j];
        for (x /= y, i = m; i <= high; i++)
          v[i][j] += x * d[i];
      }
    }
  }

  return 0;

}    /* -------------------------- orttrans ------------------------- */


static int hqrvec       /* compute eigenvectors ......................*/
                  (int     n,           /* Dimension of matrix .......*/
                   int     low,         /* first nonzero row .........*/
                   int     high,        /* last nonzero row ..........*/
                   REAL *  h[],         /* upper Hessenberg matrix ...*/
                   REAL    wr[],        /* Real parts of evalues .....*/
                   REAL    wi[],        /* Imaginary parts of evalues */
                   REAL *  eivec[]      /* Eigenvectors ..............*/
                  )
/*====================================================================*
 *                                                                    *
 *  hqrvec computes the eigenvectors for the eigenvalues found in hqr2*
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;  ( n > 0 )                                    *
 *               Dimension of  mat and eivec, number of eigenvalues.  *
 *      low      int low;                                             *
 *      high     int high; see  balance                               *
 *      h        REAL   *h[n];                                        *
 *               upper Hessenberg matrix                              *
 *      wr       REAL   wr[n];                                        *
 *               Real parts of the n eigenvalues.                     *
 *      wi       REAL   wi[n];                                        *
 *               Imaginary parts of the n eigenvalues.                *
 *                                                                    *
 *   Output parameter:                                                *
 *   ================                                                 *
 *      eivec    REAL   *eivec[n];                                    *
 *               Matrix, whose columns are the eigenvectors           *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      =  0     all ok                                               *
 *      =  1     h is the zero matrix.                                *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   function in use  :                                               *
 *   ==================                                               *
 *                                                                    *
 *      int   comdiv(): complex division                              *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Constants used  :    MACH_EPS                                    *
 *   =================                                                *
 *                                                                    *
 *   Macros :   SQR, ABS                                              *
 *   ========                                                         *
 *                                                                    *
 *====================================================================*/
{
  int   i, j, k;
  int   l, m, en, na;
  REAL  p, q, r = ZERO, s = ZERO, t, w, x, y, z = ZERO,
        ra, sa, vr, vi, norm;

  for (norm = ZERO, i = 0; i < n; i++)        /* find norm of h       */
  {
    for (j = i; j < n; j++) norm += ABS(h[i][j]);
  }

  if (norm == ZERO) return (1);               /* zero matrix          */

  for (en = n - 1; en >= 0; en--)             /* transform back       */
  {
    p = wr[en];
    q = wi[en];
    na = en - 1;
    if (q == ZERO)
    {
      m = en;
      h[en][en] = ONE;
      for (i = na; i >= 0; i--)
      {
        w = h[i][i] - p;
        r = h[i][en];
        for (j = m; j <= na; j++) r += h[i][j] * h[j][en];
        if (wi[i] < ZERO)
        {
          z = w;
          s = r;
        }
        else
        {
          m = i;
          if (wi[i] == ZERO)
            h[i][en] = -r / ((w != ZERO) ? (w) : (MACH_EPS * norm));
          else
          {  /* Solve the linear system:            */
             /* | w   x |  | h[i][en]   |   | -r |  */
             /* |       |  |            | = |    |  */
             /* | y   z |  | h[i+1][en] |   | -s |  */

             x = h[i][i+1];
             y = h[i+1][i];
             q = SQR (wr[i] - p) + SQR (wi[i]);
             h[i][en] = t = (x * s - z * r) / q;
             h[i+1][en] = ( (ABS(x) > ABS(z) ) ?
                                (-r -w * t) / x : (-s -y * t) / z);
          }
        }  /* wi[i] >= 0  */
      }  /*  end i     */
    }  /* end q = 0  */

    else if (q < ZERO)
    {
      m = na;
      if (ABS(h[en][na]) > ABS(h[na][en]))
      {
        h[na][na] = - (h[en][en] - p) / h[en][na];
        h[na][en] = - q / h[en][na];
      }
      else
       comdiv(-h[na][en], ZERO, h[na][na]-p, q, &h[na][na], &h[na][en]);

      h[en][na] = ONE;
      h[en][en] = ZERO;
      for (i = na - 1; i >= 0; i--)
      {
        w = h[i][i] - p;
        ra = h[i][en];
        sa = ZERO;
        for (j = m; j <= na; j++)
        {
          ra += h[i][j] * h[j][na];
          sa += h[i][j] * h[j][en];
        }

        if (wi[i] < ZERO)
        {
          z = w;
          r = ra;
          s = sa;
        }
        else
        {
          m = i;
          if (wi[i] == ZERO)
            comdiv (-ra, -sa, w, q, &h[i][na], &h[i][en]);
          else
          {

         /* solve complex linear system:                              */
         /* | w+i*q     x | | h[i][na] + i*h[i][en]  |   | -ra+i*sa | */
         /* |             | |                        | = |          | */
         /* |   y    z+i*q| | h[i+1][na]+i*h[i+1][en]|   | -r+i*s   | */

            x = h[i][i+1];
            y = h[i+1][i];
            vr = SQR (wr[i] - p) + SQR (wi[i]) - SQR (q);
            vi = TWO * q * (wr[i] - p);
            if (vr == ZERO && vi == ZERO)
              vr = MACH_EPS * norm *
                  (ABS (w) + ABS (q) + ABS (x) + ABS (y) + ABS (z));

            comdiv (x * r - z * ra + q * sa, x * s - z * sa -q * ra,
                    vr, vi, &h[i][na], &h[i][en]);
            if (ABS (x) > ABS (z) + ABS (q))
            {
              h[i+1][na] = (-ra - w * h[i][na] + q * h[i][en]) / x;
              h[i+1][en] = (-sa - w * h[i][en] - q * h[i][na]) / x;
            }
            else
              comdiv (-r - y * h[i][na], -s - y * h[i][en], z, q,
                                              &h[i+1][na], &h[i+1][en]);

          }   /* end wi[i] > 0  */
        }   /* end wi[i] >= 0  */
      }   /* end i            */
    }    /*  if q < 0        */
  }    /* end  en           */

  for (i = 0; i < n; i++)         /* Eigenvectors for the evalues for */
    if (i < low || i > high)      /* rows < low  and rows  > high     */
      for (k = i + 1; k < n; k++) eivec[i][k] = h[i][k];

  for (j = n - 1; j >= low; j--)
  {
    m = (j <= high) ? j : high;
    if (wi[j] < ZERO)
    {
      for (l = j - 1, i = low; i <= high; i++)
      {
        for (y = z = ZERO, k = low; k <= m; k++)
        {
          y += eivec[i][k] * h[k][l];
          z += eivec[i][k] * h[k][j];
        }

        eivec[i][l] = y;
        eivec[i][j] = z;
      }
    }
    else
      if (wi[j] == ZERO)
      {
        for (i = low; i <= high; i++)
        {
          for (z = ZERO, k = low; k <= m; k++)
            z += eivec[i][k] * h[k][j];
          eivec[i][j] = z;
        }
      }

  }  /*  end j  */

  return (0);
}


static int hqr2         /* compute eigenvalues .......................*/
                (int     vec,         /* switch for computing evectors*/
                 int     n,           /* Dimension of matrix .........*/
                 int     low,         /* first nonzero row ...........*/
                 int     high,        /* last nonzero row ............*/
                 REAL *  h[],         /* Hessenberg matrix ...........*/
                 REAL    wr[],        /* Real parts of eigenvalues ...*/
                 REAL    wi[],        /* Imaginary parts of evalues ..*/
                 REAL *  eivec[],     /* Matrix of eigenvectors ......*/
                 int     cnt[]        /* Iteration counter ...........*/
                )
/*====================================================================*
 *                                                                    *
 *  hqr2 computes the eigenvalues and (if vec != 0) the eigenvectors  *
 *  of an  n * n upper Hessenberg matrix.                             *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Control parameter:                                               *
 *   ==================                                               *
 *      vec      int vec;                                             *
 *       = 0     compute eigenvalues only                             *
 *       = 1     compute all eigenvalues and eigenvectors             *
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;  ( n > 0 )                                    *
 *               Dimension of  mat and eivec,                         *
 *               length of the real parts vector  wr and of the       *
 *               imaginary parts vector  wi of the eigenvalues.       *
 *      low      int low;                                             *
 *      high     int high; see  balance                               *
 *      h        REAL   *h[n];                                        *
 *               upper  Hessenberg matrix                             *
 *                                                                    *
 *   Output parameters:                                               *
 *   ==================                                               *
 *      eivec    REAL   *eivec[n];     ( bei vec = 1 )                *
 *               Matrix, which for vec = 1 contains the eigenvectors  *
 *               as follows  :                                        *
 *               For real eigebvalues the corresponding column        *
 *               contains the corresponding eigenvactor, while for    *
 *               complex eigenvalues the corresponding column contains*
 *               the real part of the eigenvactor with its imaginary  *
 *               part is stored in the subsequent column of eivec.    *
 *               The eigenvactor for the complex conjugate eigenvactor*
 *               is given by the complex conjugate eigenvactor.       *
 *      wr       REAL   wr[n];                                        *
 *               Real part of the n eigenvalues.                      *
 *      wi       REAL   wi[n];                                        *
 *               Imaginary parts of the eigenvalues                   *
 *      cnt      int cnt[n];                                          *
 *               vector of iterations used for each eigenvalue.       *
 *               For a complex conjugate eigenvalue pair the second   *
 *               entry is negative.                                   *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      =   0    all ok                                               *
 *      = 4xx    Iteration maximum exceeded when computing evalue xx  *
 *      =  99    zero  matrix                                         *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   functions in use  :                                              *
 *   ===================                                              *
 *                                                                    *
 *      int hqrvec(): reverse transform for eigenvectors              *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Constants used  :   MACH_EPS, MAXIT                              *
 *   =================                                                *
 *                                                                    *
 *   Macros  :   SWAP, ABS, SQRT                                      *
 *   =========                                                        *
 *                                                                    *
 *====================================================================*/
{
  int  i, j;
  int  na, en, iter, k, l, m;
  REAL p = ZERO, q = ZERO, r = ZERO, s, t, w, x, y, z;

  for (i = 0; i < n; i++)
    if (i < low || i > high)
    {
      wr[i] = h[i][i];
      wi[i] = ZERO;
      cnt[i] = 0;
    }

  en = high;
  t = ZERO;

  while (en >= low)
  {
    iter = 0;
    na = en - 1;

    for ( ; ; )
    {
      for (l = en; l > low; l--)             /* search for small      */
        if ( ABS(h[l][l-1]) <=               /* subdiagonal element   */
              MACH_EPS * (ABS(h[l-1][l-1]) + ABS(h[l][l])) )  break;

      x = h[en][en];
      if (l == en)                            /* found one evalue     */
      {
        wr[en] = h[en][en] = x + t;
        wi[en] = ZERO;
        cnt[en] = iter;
        en--;
        break;
      }

      y = h[na][na];
      w = h[en][na] * h[na][en];

      if (l == na)                            /* found two evalues    */
      {
        p = (y - x) * 0.5;
        q = p * p + w;
        z = SQRT (ABS (q));
        x = h[en][en] = x + t;
        h[na][na] = y + t;
        cnt[en] = -iter;
        cnt[na] = iter;
        if (q >= ZERO)
        {                                     /* real eigenvalues     */
          z = (p < ZERO) ? (p - z) : (p + z);
          wr[na] = x + z;
          wr[en] = s = x - w / z;
          wi[na] = wi[en] = ZERO;
          x = h[en][na];
          r = SQRT (x * x + z * z);

          if (vec)
          {
            p = x / r;
            q = z / r;
            for (j = na; j < n; j++)
            {
              z = h[na][j];
              h[na][j] = q * z + p * h[en][j];
              h[en][j] = q * h[en][j] - p * z;
            }

            for (i = 0; i <= en; i++)
            {
              z = h[i][na];
              h[i][na] = q * z + p * h[i][en];
              h[i][en] = q * h[i][en] - p * z;
            }

            for (i = low; i <= high; i++)
            {
              z = eivec[i][na];
              eivec[i][na] = q * z + p * eivec[i][en];
              eivec[i][en] = q * eivec[i][en] - p * z;
            }
          }  /* end if (vec) */
        }  /* end if (q >= ZERO) */
        else                                  /* pair of complex      */
        {                                     /* conjugate evalues    */
          wr[na] = wr[en] = x + p;
          wi[na] =   z;
          wi[en] = - z;
        }

        en -= 2;
        break;
      }  /* end if (l == na) */

      if (iter >= MAXIT)
      {
        cnt[en] = MAXIT + 1;
        return (en);                         /* MAXIT Iterations     */
      }

      if ( (iter != 0) && (iter % 10 == 0) )
      {
        t += x;
        for (i = low; i <= en; i++) h[i][i] -= x;
        s = ABS (h[en][na]) + ABS (h[na][en-2]);
        x = y = (REAL)0.75 * s;
        w = - (REAL)0.4375 * s * s;
      }

      iter ++;

      for (m = en - 2; m >= l; m--)
      {
        z = h[m][m];
        r = x - z;
        s = y - z;
        p = ( r * s - w ) / h[m+1][m] + h[m][m+1];
        q = h[m + 1][m + 1] - z - r - s;
        r = h[m + 2][m + 1];
        s = ABS (p) + ABS (q) + ABS (r);
        p /= s;
        q /= s;
        r /= s;
        if (m == l) break;
        if ( ABS (h[m][m-1]) * (ABS (q) + ABS (r)) <=
                 MACH_EPS * ABS (p)
                 * ( ABS (h[m-1][m-1]) + ABS (z) + ABS (h[m+1][m+1])) )
          break;
      }

      for (i = m + 2; i <= en; i++) h[i][i-2] = ZERO;
      for (i = m + 3; i <= en; i++) h[i][i-3] = ZERO;

      for (k = m; k <= na; k++)
      {
        if (k != m)             /* double  QR step, for rows l to en  */
        {                       /* and columns m to en                */
          p = h[k][k-1];
          q = h[k+1][k-1];
          r = (k != na) ? h[k+2][k-1] : ZERO;
          x = ABS (p) + ABS (q) + ABS (r);
          if (x == ZERO) continue;                  /*  next k        */
          p /= x;
          q /= x;
          r /= x;
        }
        s = SQRT (p * p + q * q + r * r);
        if (p < ZERO) s = -s;

        if (k != m) h[k][k-1] = -s * x;
          else if (l != m)
                 h[k][k-1] = -h[k][k-1];
        p += s;
        x = p / s;
        y = q / s;
        z = r / s;
        q /= p;
        r /= p;

        for (j = k; j < n; j++)               /* modify rows          */
        {
          p = h[k][j] + q * h[k+1][j];
          if (k != na)
          {
            p += r * h[k+2][j];
            h[k+2][j] -= p * z;
          }
          h[k+1][j] -= p * y;
          h[k][j]   -= p * x;
        }

        j = (k + 3 < en) ? (k + 3) : en;
        for (i = 0; i <= j; i++)              /* modify columns       */
        {
          p = x * h[i][k] + y * h[i][k+1];
          if (k != na)
          {
            p += z * h[i][k+2];
            h[i][k+2] -= p * r;
          }
          h[i][k+1] -= p * q;
          h[i][k]   -= p;
        }

        if (vec)      /* if eigenvectors are needed ..................*/
        {
          for (i = low; i <= high; i++)
          {
            p = x * eivec[i][k] + y * eivec[i][k+1];
            if (k != na)
            {
              p += z * eivec[i][k+2];
              eivec[i][k+2] -= p * r;
            }
            eivec[i][k+1] -= p * q;
            eivec[i][k]   -= p;
          }
        }
      }    /* end k          */

    }    /* end for ( ; ;) */

  }    /* while (en >= low)                      All evalues found    */

  if (vec)                                /* transform evectors back  */
    if (hqrvec (n, low, high, h, wr, wi, eivec)) return (99);
  return (0);
}


static int norm_1       /* normalize eigenvectors to have one norm 1 .*/
                  (int     n,       /* Dimension of matrix ...........*/
                   REAL *  v[],     /* Matrix with eigenvektors ......*/
                   REAL    wi[]     /* Imaginary parts of evalues ....*/
                  )
/*====================================================================*
 *                                                                    *
 *  norm_1 normalizes the one norm of the column vectors in v.        *
 *  (special attention to complex vectors in v  is given)             *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n; ( n > 0 )                                     *
 *               Dimension of matrix v                                *
 *      v        REAL   *v[];                                         *
 *               Matrix of eigenvectors                               *
 *      wi       REAL   wi[];                                         *
 *               Imaginary parts of the eigenvalues                   *
 *                                                                    *
 *   Output parameter:                                                *
 *   ================                                                 *
 *      v        REAL   *v[];                                         *
 *               Matrix with normalized eigenvectors                  *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      = 0      all ok                                               *
 *      = 1      n < 1                                                *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   functions used  :                                                *
 *   =================                                                *
 *      REAL   comabs():  complex absolute value                      *
 *      int    comdiv():  complex division                            *
 *                                                                    *
 *   Macros :   ABS                                                   *
 *   ========                                                         *
 *                                                                    *
 *====================================================================*/
{
  int  i, j;
  REAL maxi, tr, ti;

  if (n < 1) return (1);

  for (j = 0; j < n; j++)
  {
    if (wi[j] == ZERO)
    {
      maxi = v[0][j];
      for (i = 1; i < n; i++)
        if (ABS (v[i][j]) > ABS (maxi))  maxi = v[i][j];

      if (maxi != ZERO)
      {
        maxi = ONE / maxi;
        for (i = 0; i < n; i++) v[i][j] *= maxi;
      }
    }
    else
    {
      tr = v[0][j];
      ti = v[0][j+1];
      for (i = 1; i < n; i++)
        if ( comabs (v[i][j], v[i][j+1]) > comabs (tr, ti) )
        {
          tr = v[i][j];
          ti = v[i][j+1];
        }

      if (tr != ZERO || ti != ZERO)
        for (i = 0; i < n; i++)
          comdiv (v[i][j], v[i][j+1], tr, ti, &v[i][j], &v[i][j+1]);

      j++;                                          /* raise j by two */
    }
  }
  return (0);
}


int eigen               /* Compute all evalues/evectors of a matrix ..*/
          (
           int     vec,           /* switch for computing evectors ...*/
           int     ortho,         /* orthogonal Hessenberg reduction? */
           int     ev_norm,       /* normalize Eigenvectors? .........*/
           int     n,             /* size of matrix ..................*/
           REAL ** mat,           /* input matrix ....................*/
           REAL ** eivec,         /* Eigenvectors ....................*/
           REAL  * valre,         /* real parts of eigenvalues .......*/
           REAL  * valim,         /* imaginary parts of eigenvalues ..*/
           int   * cnt            /* Iteration counter ...............*/
          )
/*====================================================================*
 *                                                                    *
 *  The function  eigen  determines all eigenvalues and (if desired)  *
 *  all eigenvectors of a real square  n * n  matrix via the QR method*
 *  in the version of  Martin, Parlett, Peters, Reinsch and Wilkinson.*
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Literature:                                                      *
 *   ===========                                                      *
 *      1) Peters, Wilkinson: Eigenvectors of real and complex        *
 *         matrices by LR and QR triangularisations,                  *
 *         Num. Math. 16, p.184-204, (1970); [PETE70]; contribution   *
 *         II/15, p. 372 - 395 in [WILK71].                           *
 *      2) Martin, Wilkinson: Similarity reductions of a general      *
 *         matrix to Hessenberg form, Num. Math. 12, p. 349-368,(1968)*
 *         [MART 68]; contribution II,13, p. 339 - 358 in [WILK71].   *
 *      3) Parlett, Reinsch: Balancing a matrix for calculations of   *
 *         eigenvalues and eigenvectors, Num. Math. 13, p. 293-304,   *
 *         (1969); [PARL69]; contribution II/11, p.315 - 326 in       *
 *         [WILK71].                                                  *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Control parameters:                                              *
 *   ===================                                              *
 *      vec      int vec;                                             *
 *               call for eigen :                                     *
 *       = 0     compute eigenvalues only                             *
 *       = 1     compute all eigenvalues and eigenvectors             *
 *      ortho    flag that shows if transformation of mat to          *
 *               Hessenberg form shall be done orthogonally by        *
 *               `orthes' (flag set) or elementarily by `elmhes'      *
 *               (flag cleared). The Householder matrices used in     *
 *               orthogonal transformation have the advantage of      *
 *               preserving the symmetry of input matrices.           *
 *      ev_norm  flag that shows if Eigenvectors shall be             *
 *               normalized (flag set) or not (flag cleared)          *
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;  ( n > 0 )                                    *
 *               size of matrix, number of eigenvalues                *
 *      mat      REAL   *mat[n];                                      *
 *               matrix                                               *
 *                                                                    *
 *   Output parameters:                                               *
 *   ==================                                               *
 *      eivec    REAL   *eivec[n];     ( bei vec = 1 )                *
 *               matrix, if  vec = 1  this holds the eigenvectors     *
 *               thus :                                               *
 *               If the jth eigenvalue of the matrix is real then the *
 *               jth column is the corresponding real eigenvector;    *
 *               if the jth eigenvalue is complex then the jth column *
 *               of eivec contains the real part of the eigenvector   *
 *               while its imaginary part is in column j+1.           *
 *               (the j+1st eigenvector is the complex conjugate      *
 *               vector.)                                             *
 *      valre    REAL   valre[n];                                     *
 *               Real parts of the eigenvalues.                       *
 *      valim    REAL   valim[n];                                     *
 *               Imaginary parts of the eigenvalues                   *
 *      cnt      int cnt[n];                                          *
 *               vector containing the number of iterations for each  *
 *               eigenvalue. (for a complex conjugate pair the second *
 *               entry is negative.)                                  *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      =   0    all ok                                               *
 *      =   1    n < 1 or other invalid input parameter               *
 *      =   2    insufficient memory                                  *
 *      = 10x    error x from balance()                               *
 *      = 20x    error x from elmh()                                  *
 *      = 30x    error x from elmtrans()   (for vec = 1 only)         *
 *      = 4xx    error xx from hqr2()                                 *
 *      = 50x    error x from balback()    (for vec = 1 only)         *
 *      = 60x    error x from norm_1()     (for vec = 1 only)         *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Functions in use   :                                             *
 *   ===================                                              *
 *                                                                    *
 *   static int balance (): Balancing of an  n x n  matrix            *
 *   static int elmh ():    Transformation to upper Hessenberg form   *
 *   static int elmtrans(): intialize eigenvectors                    *
 *   static int hqr2 ():    compute eigenvalues/eigenvectors          *
 *   static int balback (): Reverse balancing to obtain eigenvectors  *
 *   static int norm_1 ():  Normalize eigenvectors                    *
 *                                                                    *
 *   void *vmalloc():       allocate vector or matrix                 *
 *   void vmfree():         free list of vectors and matrices         *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Constants used   :     NULL, BASIS                               *
 *   ===================                                              *
 *                                                                    *
 *====================================================================*/
{
  register int i;
  int      low, high, rc;
  REAL     *scale,
           *d = NULL;
  void     *vmblock;

  if (n < 1) return (1);                       /*  n >= 1 ............*/

  if (valre == NULL || valim == NULL || mat == NULL || cnt == NULL)
    return (1);

  for (i = 0; i < n; i++)
    if (mat[i] == NULL) return (1);

  for (i = 0; i < n; i++) cnt[i] = 0;

  if (n == 1)                                  /*  n = 1 .............*/
  {
    eivec[0][0] = ONE;
    valre[0]    = mat[0][0];
    valim[0]    = ZERO;
    return (0);
  }

  if (vec)
  {
    if (eivec == NULL) return (1);
    for (i = 0; i < n; i++)
      if (eivec[i] == NULL) return (1);
  }

  vmblock = vminit();
  scale = (REAL *)vmalloc(vmblock, VEKTOR, n, 0);
  if (! vmcomplete(vmblock))                 /* memory error         */
    return 2;

  if (vec && ortho)                          /* with Eigenvectors     */
  {                                          /* and orthogonal        */
                                             /* Hessenberg reduction? */
    d = (REAL *)vmalloc(vmblock, VEKTOR, n, 0);
    if (! vmcomplete(vmblock))
    {
      vmfree(vmblock);
      return 1;
    }
  }

                                            /* balance mat for nearly */
  rc = balance (n, mat, scale,              /* equal row and column   */
                  &low, &high, BASIS);      /* one norms              */
  if (rc)
  {
    vmfree(vmblock);
    return (100 + rc);
  }

  if (ortho)
    rc = orthes(n, low, high, mat, d);
  else
    rc = elmhes (n, low, high, mat, cnt);   /*  reduce mat to upper   */
  if (rc)                                   /*  Hessenberg form       */
  {
    vmfree(vmblock);
    return (200 + rc);
  }

  if (vec)                                  /*  initialize eivec      */
  {
    if (ortho)
      rc = orttrans(n, low, high, mat, d, eivec);
    else
      rc = elmtrans (n, low, high, mat, cnt, eivec);
    if (rc)
    {
      vmfree(vmblock);
      return (300 + rc);
    }
  }

  rc = hqr2 (vec, n, low, high, mat,        /*  execute Francis QR    */
             valre, valim, eivec, cnt);     /*  algorithm to obtain   */
  if (rc)                                   /*  eigenvalues           */
  {
    vmfree(vmblock);
    return (400 + rc);
  }

  if (vec)
  {
    rc = balback (n, low, high,             /*  reverse balancing if  */
                      scale, eivec);        /*  eigenvaectors are to  */
    if (rc)                                 /*  be determined         */
    {
      vmfree(vmblock);
      return (500 + rc);
    }
    if (ev_norm)
      rc = norm_1 (n, eivec, valim);        /* normalize eigenvectors */
    if (rc)
    {
      vmfree(vmblock);
      return (600 + rc);
    }
  }

  vmfree(vmblock);                          /* free buffers           */


  return (0);
}

/* ------------------------ END feigen.cpp -------------------------- */

---