/* nag_eigen_real_gen_quad (f02jcc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 24, 2013.
 */

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

int main(void)
{
  /* Integer scalar and array declarations */
  Integer           i, iwarn, j, pda, pdb, pdc, pdvl, pdvr, n;
  Integer           exit_status = 0;

  /* Nag Types */
  NagError          fail;
  Nag_ScaleType     scal;
  Nag_LeftVecsType  jobvl;
  Nag_RightVecsType jobvr;
  Nag_CondErrType   sense;

  /* Double scalar and array declarations */
  double            bmax, inf, tmp;
  double            tol = 0.0;
  double            *a = 0, *alphai= 0, *alphar= 0, *b = 0, *beta= 0, *bevl= 0;
  double            *bevr= 0, *c = 0, *ei= 0, *er= 0, *s= 0, *vl = 0, * vr = 0;

  /* Character scalar declarations */
  char              cjobvl[40], cjobvr[40], cscal[40], csense[40];

  /* Initialise the error structure */
  INIT_FAIL(fail);

  printf("nag_eigen_real_gen_quad (f02jcc) Example Program Results\n\n");
  fflush(stdout);

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Read in the problem size, scaling and output required */
  scanf("%ld%39s%39s%*[^\n] ", &n, cscal, csense);
  scal = (Nag_ScaleType) nag_enum_name_to_value(cscal);
  sense = (Nag_CondErrType) nag_enum_name_to_value(csense);
  scanf("%39s%39s%*[^\n] ",cjobvl,cjobvr);
  jobvl = (Nag_LeftVecsType) nag_enum_name_to_value(cjobvl);
  jobvr = (Nag_RightVecsType) nag_enum_name_to_value(cjobvr);

  pda = n;
  pdb = n;
  pdc = n;
  pdvl = n;
  pdvr = n;

  if (!(a = NAG_ALLOC(n*pda, double)) ||
      !(b = NAG_ALLOC(n*pdb, double)) ||
      !(c = NAG_ALLOC(n*pdc, double)) ||
      !(alphai = NAG_ALLOC(2*n, double)) ||
      !(alphar = NAG_ALLOC(2*n, double)) ||
      !(beta = NAG_ALLOC(2*n, double)) ||
      !(ei = NAG_ALLOC(2*n, double)) ||
      !(er = NAG_ALLOC(2*n, double)) ||
      !(vl = NAG_ALLOC(2*n*pdvl, double)) ||
      !(vr = NAG_ALLOC(2*n*pdvr, double)) ||
      !(s = NAG_ALLOC(2*n, double)) ||
      !(bevr = NAG_ALLOC(2*n, double)) ||
      !(bevl = NAG_ALLOC(2*n, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Read in the matrix A */
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      scanf("%lf",&a[j*pda+i]);
  scanf("%*[^\n] ");
      
  /* Read in the matrix B */
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      scanf("%lf",&b[j*pdb+i]);
  scanf("%*[^\n] ");

  /* Read in the matrix C */
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      scanf("%lf",&c[j*pdc+i]);
  scanf("%*[^\n] ");


  /* nag_eigen_real_gen_quad (f02jcc): Solve the quadratic eigenvalue problem */
  nag_eigen_real_gen_quad(scal, jobvl, jobvr, sense, tol, n, a, pda, b, pdb, c,
                          pdc, alphar, alphai, beta, vl, pdvl, vr, pdvr, s,
                          bevl, bevr, &iwarn, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_eigen_real_gen_quad (f02jcc).\n%s\n", fail.message);
    exit_status = -1;
    goto END;
  } else if (iwarn!=0) {
    printf("Warning from nag_eigen_real_gen_quad (f02jcc).");
    printf("  iwarn = %ld\n", iwarn);
  }

  /* Infinity */
  inf = X02ALC;

  /* Display eigenvalues */
  for (j = 0; j < 2*n; j++) {
    if (beta[j] >= 1.0) {
      er[j] = alphar[j]/beta[j];
      ei[j] = alphai[j]/beta[j];
    } else {
      tmp = inf*beta[j];
      if ((fabs(alphar[j])<tmp) && (fabs(alphai[j])<tmp)) {
        er[j] = alphar[j]/beta[j];
        ei[j] = alphai[j]/beta[j];
      } else {
        er[j] = inf;
        ei[j] = 0.0;
      }
    }
    if (er[j]<inf) {
      printf("Eigenvalue(%3ld) = (%11.4e, %11.4e)\n",j+1,er[j],ei[j]);
    } else {
      printf("Eigenvalue(%3ld) is infinte\n",j+1);
    }
  }

  if (jobvr == Nag_RightVecs) {
    printf("\n");
    fflush(stdout);
    /* x04cac: Print out the right eigenvectors */
    nag_gen_real_mat_print (Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag,
                            n, 2*n, vr, pdvr, "Right eigenvectors (matrix VR)",
                            NULL, &fail);
  }

  if (jobvl == Nag_LeftVecs && fail.code == NE_NOERROR) {
    printf("\n");
    fflush(stdout);
    /* x04cac: Print out the left eigenvectors */
    nag_gen_real_mat_print (Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, 
                            n, 2*n, vl, pdvl, "Left eigenvectors (matrix VL)",
                            NULL, &fail);
  }
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  /* Display condition numbers and errors, as required */
  if (sense==Nag_CondOnly || sense==Nag_CondBackErrLeft || 
      sense==Nag_CondBackErrRight || sense==Nag_CondBackErrBoth) {
    printf("\n");
    printf("Eigenvalue Condition numbers\n");
    for (j = 0 ; j < 2*n; j++)
      printf("%2ld  %11.4e\n", j+1, s[j]);
  }

  if (sense==Nag_BackErrRight || sense==Nag_BackErrBoth || 
      sense==Nag_CondBackErrRight || sense==Nag_CondBackErrBoth) {
    /* nag_dmax_val (f16jnc).
     * Get maximum value (bmax) and location of that value (j) of bevr.
     */
    nag_dmax_val(2*n, bevr, 1, &j, &bmax, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_dmax_val (f16jnc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    printf("\n");
    printf("Backward errors for eigenvalues and right eigenvectors\n");
    if (bmax<10.0*X02AJC) {
      printf("  All errors are less than 10 times machine precision\n");
    } else {
      for (j = 0; j < 2*n; j++) 
        printf("%11.4e\n", bevr[j]);
    }
  }

  if (sense==Nag_BackErrLeft || sense==Nag_BackErrBoth || 
      sense==Nag_CondBackErrLeft || sense==Nag_CondBackErrBoth) {
    /* nag_dmax_val (f16jnc).
     * Get maximum value (bmax) and location of that value (j) of bevl.
     */
    nag_dmax_val(2*n, bevl, 1, &j, &bmax, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_dmax_val (f16jnc).\n%s\n", fail.message);
      exit_status = 2;
      goto END;
    }
    printf("\n");
    printf("Backward errors for eigenvalues and left eigenvectors\n");
    if (bmax<10.0*X02AJC) {
      printf("  All errors are less than 10 times machine precision\n");
    } else {
      for (j = 0; j < 2*n; j++) 
        printf("%11.4e\n", bevl[j]);
    }
  }

 END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(alphai);
  NAG_FREE(alphar);
  NAG_FREE(beta);
  NAG_FREE(ei);
  NAG_FREE(er);
  NAG_FREE(vl);
  NAG_FREE(vr);
  NAG_FREE(s);
  NAG_FREE(bevr);
  NAG_FREE(bevl);

  return(exit_status);
}