/* nag_dtrsv (f16pjc) Example Program.

 *

 * Copyright 2005 Numerical Algorithms Group.

 *

 * Mark 8, 2005.

 */



#include <stdio.h>

#include <nag.h>

#include <nag_stdlib.h>

#include <nagf16.h>

#include <nagx04.h>



int main(void)

{



  /* Scalars */

  double alpha;

  Integer exit_status, i, incx, j, n, pda, xlen;



  /* Arrays */

  double *a=0, *x=0;

  char nag_enum_arg[40];



  /* Nag Types */

  NagError fail;

  Nag_OrderType order;

  Nag_TransType trans;

  Nag_UploType uplo;

  Nag_DiagType diag;



#ifdef NAG_COLUMN_MAJOR

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

  order = Nag_ColMajor;

#else

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

  order = Nag_RowMajor;

#endif



  exit_status = 0;

  INIT_FAIL(fail);



  Vprintf( "nag_dtrsv (f16pjc) Example Program Results\n\n");



  /* Skip heading in data file */

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



  /* Read the problem dimensions */

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



  /* Read the uplo storage parameter */

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

  /* nag_enum_name_to_value(x04nac).

   * Converts NAG enum member name to value

   */

  uplo = nag_enum_name_to_value(nag_enum_arg);

  /* Read the transpose parameter */

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

  /* nag_enum_name_to_value(x04nac), see above. */

  trans = nag_enum_name_to_value(nag_enum_arg);

  /* Read the unit-diagonal parameter */

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

  /* nag_enum_name_to_value(x04nac), see above. */

  diag = nag_enum_name_to_value(nag_enum_arg);



  /* Read scalar parameters */

  Vscanf("%lf%*[^\n] ", &alpha);

  /* Read increment parameter */

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



  pda = n;

  xlen = MAX(1, 1 + (n - 1)*ABS(incx));



  if (n > 0)

    {

      /* Allocate memory */

      if ( !(a = NAG_ALLOC(pda*n, double)) ||

           !(x = NAG_ALLOC(xlen, double)))

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

    }

  else

    {

      Vprintf("Invalid n\n");

      exit_status = 1;

      return exit_status;

    }



  /* Input matrix A and vector x*/



  if (uplo == Nag_Upper)

    {

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

        {

          if (diag == Nag_NonUnitDiag)

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

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

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

        }

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

    }

  else

    {

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

        {

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

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

          if (diag == Nag_NonUnitDiag)

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

        }

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

    }

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

    Vscanf("%lf%*[^\n] ", &x[i]);



  /* nag_dtrsv(f16pjc).

   * Solution of real triangular system of linear equations.

   *

   */

  nag_dtrsv(order, uplo, trans, diag, n, alpha, a, pda, x, incx,

            &fail);

  if (fail.code != NE_NOERROR)

    {

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

      exit_status = 1;

      goto END;

    }



  /* Print output vector x */

  Vprintf("%s\n", " Solution x:");

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

    {

      Vprintf("%11f\n", x[i]);

    }



 END:

  if (a) NAG_FREE(a);

  if (x) NAG_FREE(x);



  return exit_status;

}