/* nag_complex_lu_solve_mult_rhs (f04akc) Example Program.

 *

 * Copyright 1990 Numerical Algorithms Group.

 *

 * Mark 1A revised, (Oct 1990).

 * Mark 8 revised, 2004.

 */



#include <nag.h>

#include <stdio.h>

#include <nag_stdlib.h>

#include <nagf03.h>

#include <nagf04.h>



#define A(I,J) a[(I)*tda + J]

#define B(I,J) b[(I)*tdb + J]

int main(void)

{



  Complex *a=0, *b=0, det;

  Integer dete, exit_status=0, i, j, n, nrhs, *pivot=0, tda, tdb;

  NagError fail;



  INIT_FAIL(fail);



  Vprintf("nag_complex_lu_solve_mult_rhs (f04akc) Example Program Results\n");

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

  Vscanf("%ld",&n);

  nrhs = 1;

  if (n>=1)

    {

      if ( !( pivot = NAG_ALLOC(n, Integer))  ||

           !( a = NAG_ALLOC(n*n, Complex)) ||

           !( b = NAG_ALLOC(n*nrhs, Complex)) )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

      tda = n;

      tdb = nrhs;

    }

  else

    {

      Vprintf("Invalid n.\n");

      exit_status = 1;

      return exit_status;

    }

  for (i = 0; i < n; i++)

    for (j = 0; j < n; j++)

      Vscanf(" ( %lf,  %lf ) ", &A(i,j).re, &A(i,j).im);

  for (i = 0; i < n; i++)

    for (j = 0; j < nrhs; j++)

      Vscanf(" ( %lf,  %lf ) ", &B(i,j).re, &B(i,j).im);

  /* nag_complex_lu (f03ahc).

   * LU factorization and determinant of complex matrix

   */

  nag_complex_lu(n, a, tda, pivot, &det, &dete, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_complex_lu (f03ahc).\n%s\n", fail.message);

      exit_status = 1;

      goto END;

    }

  /* nag_complex_lu_solve_mult_rhs (f04akc).

   * Approximate solution of complex simultaneous linear

   * equations (coefficient matrix already factorized by

   * nag_complex_lu (f03ahc))

   */

  nag_complex_lu_solve_mult_rhs(n, nrhs, a, tda, pivot, b, tdb, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_complex_lu_solve_mult_rhs (f04akc).\n%s\n",

              fail.message);

      exit_status = 1;

      goto END;

    }

  Vprintf("Solution\n");

  for (i=0; i<n; i++)

    {

      for (j=0; j<nrhs; j++)

        Vprintf("(%7.3f, %7.3f) ", B(i,j).re, B(i,j).im);

      Vprintf("\n");

    }

 END:

  if (pivot) NAG_FREE(pivot);

  if (a) NAG_FREE(a);

  if (b) NAG_FREE(b);

  return exit_status;

}