/* nag_real_general_eigensystem(f02bjc) Example Program
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nagf02.h>
#include <nagx02.h>

#define NMAX 8
#define TDA NMAX
#define TDB NMAX
#define TDZ NMAX

int main(void)
{
  Integer i, j, k, n, ip, iter[NMAX];
  Complex alfa[NMAX];
  double  beta[NMAX], eps1;
  double a[NMAX][TDA] , b[NMAX][TDB], z[NMAX][TDZ];
  Boolean matz;

  Vprintf("f02bjc Example Program Results\n");
  Vscanf("%*[^\n]"); /* Skip heading in data file */
  Vscanf("%ld", &n);
  if (n>0 && n<=NMAX)
    {
      for (i=0; i<n; ++i)
        for (j=0; j<n; ++j)
          Vscanf("%lf", &a[i][j]);
      for (i=0; i<n; ++i)
        for (j=0; j<n; ++j)
          Vscanf("%lf",&b[i][j]);
      matz = TRUE;
      eps1 = X02AJC;
      f02bjc(n, &a[0][0], (Integer)TDA, &b[0][0], (Integer)TDB, eps1,
             alfa, beta, matz, &z[0][0], (Integer)TDZ, iter,
             NAGERR_DEFAULT);
      ip = 0;
      for (i=0; i<n; ++i)
        {
          Vprintf("Eigensolution  %4ld\n",i+1);
          Vprintf("alfa[%ld].re %7.3f",i,alfa[i].re);
          Vprintf(" alfa[%ld].im %7.3f",i,alfa[i].im);
          Vprintf(" beta[%ld] %7.3f\n",i,beta[i]);
          if (beta[i] == 0.0)
            Vprintf("lambda is infinite");
          else
            if (alfa[i].im == 0.0)
              {
                Vprintf("lambda    %7.3f\n",alfa[i].re/beta[i]);
                Vprintf("Eigenvector\n");
                for (j=0; j<n; ++j)
                  Vprintf("%7.3f\n", z[j][i]);
              }
            else
              {
                Vprintf("lambda    %7.3f   %7.3f\n",
                        alfa[i].re/beta[i], alfa[i].im/beta[i] );
                Vprintf("Eigenvector\n");
                k = (Integer)pow((double)-1, (double)(ip+2));
                for (j=0; j<n; ++j)
                  {
                    Vprintf("%7.3f",z[j][i-ip]);
                    Vprintf("%7.3f\n",k*z[j][i-ip+1]);
                  }
                ip = 1-ip;
              }
        }
      Vprintf("Number of iterations (machine-dependent)\n");
      for (i=0; i<n; ++i)
        Vprintf("%2ld",iter[i]);
      Vprintf("\n");
      return EXIT_SUCCESS;
    }
  else
    {
      Vfprintf(stderr,"n is out of range: n = %4ld,\n",n);
      return EXIT_FAILURE;
    }
}