/* nag_complex_gen_lin_solve (f04cac) Example Program.

 *

 * Copyright 2004 Numerical Algorithms Group.

 *

 * Mark 8, 2004.

 */



#include <stdio.h>

#include <nag.h>

#include <nag_stdlib.h>

#include <nagf04.h>

#include <nagx04.h>



int main(void)

{



  /* Scalars */

  double errbnd, rcond;

  Integer exit_status, i,  j, n, nrhs, pda, pdb;



  /* Arrays */

  char *clabs=0, *rlabs=0;

  Complex *a=0, *b=0;

  Integer *ipiv=0;



  /* Nag types */

  NagError fail;

  Nag_OrderType order;



#ifdef NAG_COLUMN_MAJOR

#define A(I,J) a[(J-1)*pda + I - 1]

#define B(I,J) b[(J-1)*pdb + I - 1]

  order = Nag_ColMajor;

#else

#define A(I,J) a[(I-1)*pda + J - 1]

#define B(I,J) b[(I-1)*pdb + J - 1]

  order = Nag_RowMajor;

#endif



  exit_status = 0;

  INIT_FAIL(fail);



  Vprintf("nag_complex_gen_lin_solve (f04cac) Example Program Results\n\n");



  /* Skip heading in data file */

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



  Vscanf("%ld%ld%*[^\n] ",   &n,  &nrhs);

  if (n>0 && nrhs>0)

    {

      /* Allocate memory */

      if ( !(clabs = NAG_ALLOC(2, char)) ||

           !(rlabs = NAG_ALLOC(2, char)) ||

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

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

           !(ipiv = NAG_ALLOC(8, Integer)) )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

#ifdef NAG_COLUMN_MAJOR

      pda = n;

      pdb = n;

#else

      pda = n;

      pdb = nrhs;

#endif



      /* Read A and B from data file */

    }

  else

    {

      Vprintf("%s\n",   "NMAX and/or NRHSMX too small");

      exit_status = 1;

      return exit_status;



    }



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

    {

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

        {

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

        }

    }

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



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

    {

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

        {

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

        }

    }

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





  /* Solve the equations AX = B for X */

  /* nag_complex_gen_lin_solve (f04cac).

   * Computes the solution and error-bound to a complex system

   * of linear equations

   */

  nag_complex_gen_lin_solve(order, n, nrhs, a, pda, ipiv, b, pdb, &rcond,

                            &errbnd, &fail);

  if (fail.code == NE_NOERROR)

    {



      /* Print solution, estimate of condition number and approximate */

      /* error bound */



      /* nag_gen_complx_mat_print_comp (x04dbc).

       * Print complex general matrix (comprehensive)

       */

      nag_gen_complx_mat_print_comp(order, Nag_GeneralMatrix, Nag_NonUnitDiag,

                                    n, nrhs, b, pdb, Nag_BracketForm, "%7.4f",

                                    "Solution", Nag_IntegerLabels, 0,

                                    Nag_IntegerLabels, 0, 80, 0, 0, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_gen_complx_mat_print_comp (x04dbc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      Vprintf("\n");

      Vprintf("%s\n%8s%9.1e\n", "Estimate of condition number", "", 1.0/rcond);

      Vprintf("\n\n");

      Vprintf("%s\n%8s%9.1e\n\n",

              "Estimate of error bound for computed solutions", "", errbnd);

    }

  if (fail.code == NE_RCOND)

    {



      /* Matrix is numerically singular.  Print estimate of */

      /* reciprocal of condition number and solution */



      Vprintf("\n");

      Vprintf("%s\n%8s%9.1e\n\n\n",

              "Estimate of reciprocal of condition number", "", rcond);

      /* nag_gen_complx_mat_print_comp (x04dbc), see above. */

      nag_gen_complx_mat_print_comp(order, Nag_GeneralMatrix, Nag_NonUnitDiag,

                                    n, nrhs, b, pdb, Nag_BracketForm, 0,

                                    "Solution", Nag_IntegerLabels, 0,

                                    Nag_IntegerLabels, 0, 80, 0, 0, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_gen_complx_mat_print_comp (x04dbc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }

    }

  if (fail.code == NE_SINGULAR)

    {



      /* The upper triangualr matrix U is exactly singular.  Print */

      /* details of factorization */

      Vprintf("\n");

      /* nag_gen_complx_mat_print_comp (x04dbc), see above. */

      nag_gen_complx_mat_print_comp(order, Nag_GeneralMatrix, Nag_NonUnitDiag,

                                    n, n, a, pda, Nag_BracketForm, "%7.4f",

                                    "Details of factorization",

                                    Nag_IntegerLabels, 0, Nag_IntegerLabels, 0,

                                    80, 0, 0, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_gen_complx_mat_print_comp (x04dbc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      /* Print pivot indices */

      Vprintf("\n");

      Vprintf("%s\n",   "Pivot indices");

      Vprintf("%1s", "");

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

        {

          Vprintf("%11ld%s", ipiv[i - 1], i%7 == 0 || i == n ?"\n":" ");

        }

      Vprintf("\n");

    }

 END:

  if (clabs) NAG_FREE(clabs);

  if (rlabs) NAG_FREE(rlabs);

  if (a) NAG_FREE(a);

  if (b) NAG_FREE(b);

  if (ipiv) NAG_FREE(ipiv);



  return exit_status;

}



#undef B

#undef A