/*
  This source file is adapted from feigen.c that comes with the book
  Numeric Algorithm with C by Frank Uhlig et al. I cannot find the
  license of the original source codes. I release my modifications under
  the MIT license. The modified version requires C99 as it uses complex
  numbers. I may modify the code if this is a concern.
 */

/* The MIT License

   Copyright (c) 1996 Frank Uhlig et al.
                 2009 Genome Research Ltd (GRL).

   Permission is hereby granted, free of charge, to any person obtaining
   a copy of this software and associated documentation files (the
   "Software"), to deal in the Software without restriction, including
   without limitation the rights to use, copy, modify, merge, publish,
   distribute, sublicense, and/or sell copies of the Software, and to
   permit persons to whom the Software is furnished to do so, subject to
   the following conditions:

   The above copyright notice and this permission notice shall be
   included in all copies or substantial portions of the Software.

   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   SOFTWARE.
*/

/* Contact: Heng Li <lh3lh3@gmail.com> */


/*.FE{C 7.8}
     {QR Algorithm}
     {Eigenvalues and Eigenvectors of a Matrix via the
      QR Algorithm}*/

/*.BE*/
/* ------------------------- MODULE feigen.c ------------------------ */

#include <math.h>
#include <float.h> /* for DBL_EPSILON */
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>

#define REAL   double
#define TRUE   1
#define FALSE  0
#define ZERO   0.
#define ONE    1.
#define TWO    2.
#define VEKTOR 0
#define MACH_EPS DBL_EPSILON

#define ABS(x) (((x) >= 0.)? (x) : -(x))
#define SQRT(x) sqrt(x)
#define SQR(x) ((x) * (x))
#define SWAP(typ, a, b) { typ t; t = (a); (a) = (b); (b) = t; }
#define BASIS basis()

typedef struct {
	int n, max;
	REAL *mem;
} vmblock_t;

static void *vminit(void)
{
	return (vmblock_t*)calloc(1, sizeof(vmblock_t));
}
static int vmcomplete(void *vmblock)
{
	return 1;
}
static void vmfree(void *vmblock)
{
	vmblock_t *vmb = (vmblock_t*)vmblock;
	free(vmb->mem); free(vmblock);
}
static void *vmalloc(void *vmblock, int typ, size_t zeilen, size_t spalten)
{
	vmblock_t *vmb = (vmblock_t*)vmblock;
	double *ret = 0;
	if (typ == 0) {
		if (vmb->n + zeilen > vmb->max) {
			vmb->max = vmb->n + zeilen;
			vmb->mem = (REAL*)realloc(vmb->mem, vmb->max * sizeof(REAL));
		}
		ret = vmb->mem + vmb->n;
		vmb->n += zeilen;
	}
	return ret;
}

static int basis(void) /* find basis used for computer number representation */
/*.IX{basis}*/

/***********************************************************************
* Find the basis for representing numbers on the computer, if not      *
* already done.                                                        *
*                                                                      *
* global names used:                                                   *
* ==================                                                   *
* REAL, ONE, TWO                                                       *
***********************************************************************/

{
  REAL x,
       eins,
       b;

  x = eins = b = ONE;

  while ((x + eins) - x == eins)
    x *= TWO;
  while ((x + b) == x)
    b *= TWO;


  return (int)((x + b) - x);
}


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


/*--------------------------------------------------------------------*
 * Aux functions for  eigen                                           *
 *--------------------------------------------------------------------*/


static int balance       /* balance a matrix .........................*/
/*.IX{balance}*/
                   (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 ........................*/
/*.IX{balback}*/
                   (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 ....*/
/*.IX{elmhes}*/
                  (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 .................*/
/*.IX{elmtrans}*/
                    (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 */
/*.IX{orthes}*/
                 (
                  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 */
/*.IX{orttrans}*/
                   (
                    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.                       *
*                                                                      *
* benutzte globale Namen:                                              *
* =======================                                              *
* 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 ......................*/
/*.IX{hqrvec}*/
                  (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]); */
		  REAL complex c;
		  c = -h[na][en] / (h[na][na]-p + q * I);
		  h[na][na] = creal(c); h[na][en] = cimag(c);
	  }

      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]); */
			  REAL complex c;
			  c = (-ra - sa * I) / (w + q * I);
			  h[i][na] = creal(c); h[i][en] = cimag(c);
          } 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]); */
				REAL complex c;
				c = (x * r - z * ra + q * sa + I * (x * s - z * sa -q * ra)) / (vr + I * vi);
				h[i][na] = creal(c); h[i][en] = cimag(c);
			}
            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]); */
				REAL complex c;
				c = (-r - y * h[i][na] + I * (-s - y * h[i][en])) / (z + I * q);
				h[i+1][na] = creal(c); h[i+1][en] = cimag(c);
			}

          }   /* 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 .......................*/
/*.IX{hqr2}*/
                (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 .*/
/*.IX{norm\unt 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) ) */
		  if (cabs(v[i][j] + I * v[i][j+1]) > cabs(tr + I * 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]); */
			  REAL complex c;
			  c = (v[i][j] + I * v[i][j+1]) / (tr + I * ti);
			  v[i][j] = creal(c); v[i][j+1] = cimag(c);
		  }

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

/*.BA*/

static int eigen        /* Compute all evalues/evectors of a matrix ..*/
/*.IX{eigen}*/
          (
           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.*
 *                                                                    *
 *====================================================================*
.BE*)
 *                                                                    *
 *   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                               *
 *   ===================                                              *
 *                                                                    *
 *====================================================================*/
{
  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.c -------------------------- */

/* _a[0..n^2-1] is a real general matrix. On return, evalr store the
   real part of eigenvalues and evali the imgainary part. If _evec is
   not a NULL pointer, eigenvectors will be stored there. */
int n_eigeng(double *_a, int n, double *evalr, double *evali, double *_evec)
{
	double **a, **evec = 0;
	int i, j, *cnt;
	a = (double**)calloc(n, sizeof(void*));
	if (_evec) evec = (double**)calloc(n, sizeof(void*));
	cnt = (int*)calloc(n, sizeof(int));
	for (i = 0; i < n; ++i) {
		a[i] = _a + i * n;
		if (_evec) evec[i] = _evec + i * n;
	}
	eigen(_evec? 1 : 0, 0, 1, n, a, evec, evalr, evali, cnt);
	if (_evec) {
		double tmp;
		for (j = 0; j < n; ++j) {
			tmp = 0.0;
			for (i = 0; i < n; ++i) tmp += SQR(evec[i][j]);
			tmp = SQRT(tmp);
			for (i = 0; i < n; ++i) evec[i][j] /= tmp;
		}
	}
	free(a); free(evec); free(cnt);
	return 0;
}

#ifdef LH3_EIGENG_MAIN
int main(void)
{
	static double   a[5][5] = {{1.0, 6.0, -3.0, -1.0, 7.0},
							   {8.0, -15.0, 18.0, 5.0, 4.0},
							   {-2.0, 11.0, 9.0, 15.0, 20.0},
							   {-13.0, 2.0, 21.0, 30.0, -6.0},
							   {17.0, 22.0, -5.0, 3.0, 6.0}};

	static double b[5][5]={ {10.0,1.0,2.0,3.0,4.0},
							{1.0,9.0,-1.0,2.0,-3.0},
							{2.0,-1.0,7.0,3.0,-5.0},
							{3.0,2.0,3.0,12.0,-1.0},
							{4.0,-3.0,-5.0,-1.0,15.0}};

	double *mem, *evalr, *evali, *evec;
	int i, j;
	mem = (double*)calloc(5 * 5 + 10, sizeof(double));
	evec = mem;
	evalr = evec + 25;
	evali = evalr + 5;
//	n_eigeng(a[0], 5, evalr, evali, evec);
	n_eigeng(a[0], 5, evalr, evali, 0);
	for (i = 0; i < 5; ++i) {
		printf("%le + %le J\n", evalr[i], evali[i]);
	}

	printf("\n\n");
	n_eigeng(b[0], 5, evalr, evali, evec);
	for (i = 0; i <= 4; i++)
		printf("%13.7e + %13.7e J\n", evalr[i], evali[i]);
	printf("\n");
	for (i = 0; i <= 4; i++) {
		for (j = 0; j <= 4; j++)
			printf("%12.6e ", evec[i*5 + j]);
		printf("\n");
	}
	printf("\n");

	free(mem);
	return 0;
}
#endif