/* nag_hermitian_lin_eqn_mult_rhs (f04awc) Example Program.

 *

 * Copyright 1990 Numerical Algorithms Group.

 *

 * Mark 1, 1990.

 * Mark 8 revised, 2004.

 */



#include <nag.h>

#include <stdio.h>

#include <nag_stdlib.h>

#include <nagf01.h>

#include <nagf04.h>





#define COMPLEX(A) A.re, A.im

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

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

#define X(I,J) x[(I)*tdx + J]

int main(void)

{

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

  Integer exit_status=0, i, j, n, nrhs, tda, tdb, tdx;

  NagError fail;

  double *p=0;



  INIT_FAIL(fail);

  Vprintf("nag_hermitian_lin_eqn_mult_rhs (f04awc) Example Program Results\n");

  /* Skip heading in data file */

  Vscanf("%*[^\n]");

  Vscanf("%ld", &n);

  nrhs = 1;

  if (n>=1)

    {

      if ( !( p = NAG_ALLOC(n, double))  ||

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

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

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

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

      tda = n;

      tdb = nrhs;

      tdx = nrhs;

    }

  else

    {

      Vprintf("Invalid n.\n");

      exit_status = 1;

      return exit_status;

    }

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

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

      Vscanf(" ( %lf,  %lf ) ", COMPLEX(&A(i,j)));

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

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

      Vscanf(" ( %lf,  %lf ) ", COMPLEX(&B(i,j)));

  /* Factorize A */

  /* nag_complex_cholesky (f01bnc).

   * UU^H factorization of complex Hermitian positive-definite

   * matrix

   */

  nag_complex_cholesky(n, a, tda, p, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_complex_cholesky (f01bnc).\n%s\n", fail.message);

      exit_status = 1;

      goto END;

    }

  /* Solve A */

  /* nag_hermitian_lin_eqn_mult_rhs (f04awc).

   * Approximate solution of complex Hermitian

   * positive-definite simultaneous linear equations

   * (coefficient matrix already factorized by

   * nag_complex_cholesky (f01bnc))

   */

  nag_hermitian_lin_eqn_mult_rhs(n, nrhs, a, tda, p, b, tdb, x, tdx, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_hermitian_lin_eqn_mult_rhs (f04awc).\n%s\n",

              fail.message);

      exit_status = 1;

      goto END;

    }

  Vprintf("\nSolution\n\n");

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

    {

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

        Vprintf("  (%7.4f,%7.4f)", COMPLEX(X(i,j)));

      Vprintf("\n");

    }

 END:

  if (p) NAG_FREE(p);

  if (a) NAG_FREE(a);

  if (b) NAG_FREE(b);

  if (x) NAG_FREE(x);

  return exit_status;

}