/* F08GA_T1W_F C++ Header Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 * Mark 27, 2019.
 */

#include <nag.h>
#include <dco.hpp>
#include <nagad.h>
#include <stdio.h>
#include <math.h>
#include <iostream>
#include <string>

typedef double DCO_BASE_TYPE;
typedef dco::gt1s<DCO_BASE_TYPE> DCO_MODE;
typedef DCO_MODE::type DCO_TYPE;


int main(void)
{
  /* Scalars */
  double      eerrbd, eps;
  Integer     exit_status = 0, i, j, n;
  /* Arrays */
  DCO_TYPE    *ap = 0, *dummy = 0, *w = 0, *work = 0, *ap_in = 0;
  Integer     ifail;
  std::string uplo = "U";
  double      *wr = 0, *dwda = 0;

#define AP(I, J) ap_in[J * (J - 1) / 2 + I - 1]

  printf("F08GA_T1W_F C++ Example Program Results\n\n");

  void    *ad_handle = 0;
  ifail = 0;
  x10aa_t1w_f_(ad_handle,ifail);

  /* Skip heading in data file */
#ifdef _WIN32
  scanf_s("%*[^\n]");
#else
  scanf("%*[^\n]");
#endif
#ifdef _WIN32
  scanf_s("%" NAG_IFMT "%*[^\n]", &n);
#else
  scanf("%" NAG_IFMT "%*[^\n]", &n);
#endif



  /* Allocate memory */
  ap = new DCO_TYPE [n * (n + 1) / 2];
  ap_in = new DCO_TYPE [n * (n + 1) / 2];
  dummy = new DCO_TYPE [1];
  w = new DCO_TYPE [n];
  work = new DCO_TYPE [3*n];
  wr = new double [n];
  dwda = new double [n*n];

  /* Read the upper or lower triangular part of the matrix A from data file */
  for (i = 1; i <= n; ++i) {
    for (j = i; j <= n; ++j) {
#ifdef _WIN32
      scanf_s("%lf", &dco::value(AP(i, j)));
#else
      scanf("%lf", &dco::value(AP(i, j)));
#endif
    }
  }
#ifdef _WIN32
  scanf_s("%*[^\n]");
#else
  scanf("%*[^\n]");
#endif

  for (int i = 0; i < (n*(n+1))/2; i ++) {
    dco::derivative(ap_in[i]) = 0.0;
    ap[i] = ap_in[i];
  }

  for (int j = 0; j < n; j ++) {
    dco::derivative(ap[j*(j+3)/2]) = 1.0;

    f08ga_t1w_f_(ad_handle, "N", "U", n, ap, w, dummy, 1, work, ifail, 1, 1);
    if (ifail != 0) {
      printf("Error from f08ga_t1w_f_.\n%" NAG_IFMT " ", ifail);
      exit_status = 1;
      goto END;
    }

    if (j == 0) {
      for (i = 0; i < n; ++i) {
        wr[i] = dco::value(w[i]);
      }
    }

    for (int i = 0; i < n; i++)
      dwda[i+j*n] = dco::derivative(w[i]);

    for (int i = 0;i < (n*(n+1))/2; i++)
      ap[i] = ap_in[i];
  }

  /* Print solution */
  printf("Eigenvalues\n");
  for (j = 0; j < n; ++j)
    printf("%8.4f%s", wr[j], (j + 1) % 8 == 0 ? "\n" : " ");
  printf("\n");

  /* Get the machine precision, eps, using nag_machine_precision (X02AJC)
   * and compute the approximate error bound for the computed eigenvalues.
   * Note that for the 2-norm, ||A|| = max {|w[i]|, i=0..n-1}, and since
   * the eigenvalues are in ascending order ||A|| = max( |w[0]|, |w[n-1]|).
   */
  eps = X02AJC;
  eerrbd = eps * MAX(fabs(wr[0]), fabs(wr[n - 1]));

  /* Print the approximate error bound for the eigenvalues */
  printf("\nError estimate for the eigenvalues\n");
  printf("%11.1e\n", dco::value(eerrbd));
  
  printf("\nDerivatives of eigenvalues w.r.t. diagonals of A\n");
  NagError  fail;
  INIT_FAIL(fail);
  x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,n,n,dwda,n,
         "  dW_i/dA_jj",0,&fail);

END:
  delete [] ap;
  delete [] ap_in;
  delete [] dummy;
  delete [] w;
  delete [] work;
  delete [] wr;
  delete [] dwda;

  // Remove computational data object
  x10ab_t1w_f_(ad_handle,ifail);


  return exit_status;
}

#undef AP