/* nag_complex_lu(f03ahc) Example Program
 *
 * Copyright 1990 Numerical Algorithms Group.
 *
 * Mark 1A revised, (Oct 1990).
 */

#include <nag.h>
#include <math.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nagf03.h>

int main(void)
{
#define NMAX 5
#define TDA NMAX

  Integer pivot[NMAX];
  Complex a[NMAX][TDA], det;
  Integer i, j, n, dete;
  static NagError fail;

  Vprintf("f03ahc Example Program Results\n");
  Vscanf("%*[^\n]"); /* Skip heading in data file */
  fail.print = TRUE;

  Vscanf("%ld",&n);
  if (n > 0 && n <= NMAX)
    {
      for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
          Vscanf(" ( %lf , %lf ) ", &a[i][j].re, &a[i][j].im);
      f03ahc(n, &a[0][0], (Integer)TDA, pivot, &det, &dete, &fail);
      if (fail.code!=NE_NOERROR)
        return EXIT_FAILURE;
      else
        {
          Vprintf("Array a after factorization\n");
          for (i=0; i<n; i++)
            {
              for (j=0; j<n; j++)
                Vprintf("(%7.3f, %7.3f) ", a[i][j].re, a[i][j].im);
              Vprintf("\n");
            }
          Vprintf("\nArray pivot\n");
          for (i=0; i<n; i++)
            Vprintf("%5ld",pivot[i]);
          Vprintf("\n\ndet.re = %7.4f, det.im = %7.4f, dete = %2ld.\n",
                  det.re, det.im, dete);
          det.re = ldexp(det.re, (int)dete);
          det.im = ldexp(det.im, (int)dete);
          Vprintf("\nValue of determinant = (%7.4f, %7.4f)\n", det.re, det.im);
        }
    }
  else
    {
      Vfprintf(stderr, "Error: n is out of range: n = %5ld\n", n);
      return EXIT_FAILURE;
    }
  return EXIT_SUCCESS;
}