/* nag_dsyr2 (f16prc) 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, beta;

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



  /* Arrays */

  double *a=0, *x=0, *y=0;

  char nag_enum_arg[40];



  /* Nag Types */

  NagError fail;

  Nag_OrderType order;

  Nag_UploType uplo;

  Nag_MatrixType matrix;



#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_dsyr2 (f16prc) Example Program Results\n\n");



  /* Skip heading in data file */

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



  /* Read the problem dimension */

  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 scalar parameters */

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

  /* Read increment parameters */

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



  pda = n;



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

  ylen = MAX(1, 1 + (n - 1)*ABS(incy));



  if (n > 0)

    {



      /* Allocate memory */

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

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

           !(y = NAG_ALLOC(ylen, 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)

        {

          for (j = i; 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));

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

        }

    }

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

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

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

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



  /* nag_dsyr2(f16prc).

   * Rank two update of real symmetric matrix.

   *

   */

  nag_dsyr2(order, uplo, n, alpha, x, incx, y, incy, beta, a, pda,

            &fail);

  if (fail.code != NE_NOERROR)

    {

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

      exit_status = 1;

      goto END;

    }



  if (uplo == Nag_Upper)

    {

      matrix = Nag_UpperMatrix;

    }

  else

    {

       matrix = Nag_LowerMatrix;

    }

  /* Print updated matrix A */

  /* nag_gen_real_mat_print (x04cac).

   * Print real general matrix (easy-to-use)

   */

  nag_gen_real_mat_print(order, matrix, Nag_NonUnitDiag, n,

                         n, a, pda, "Updated Matrix A", 0, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_gen_real_mat_print (x04cac).\n%s\n",

              fail.message);

      exit_status = 1;

      goto END;

    }



 END:

  if (a) NAG_FREE(a);

  if (x) NAG_FREE(x);

  if (y) NAG_FREE(y);



  return exit_status;

}