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

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

typedef double DCO_BASE_TYPE;
typedef dco::ga1s<DCO_BASE_TYPE> DCO_MODE;
typedef DCO_MODE::type DCO_TYPE;
typedef DCO_MODE::tape_t DCO_TAPE_TYPE;

int main(void)
{
  /* Integer scalar and array declarations */
  Integer  exit_status = 0;
  Integer  feval, i, iter, j, maxit, maxits, n;
  Integer  pdg, pdx;

  /* Double scalar and array declarations */
  DCO_TYPE errtol, nrmgrd;
  DCO_TYPE *g = 0, *x = 0, *g_in = 0;
  double   *xr =0, *dxdg = 0;
  Integer  ifail;

  // Create AD tape
  DCO_MODE::global_tape=DCO_TAPE_TYPE::create();
  void     *ad_handle = 0;
  x10aa_a1w_f_(ad_handle,ifail);

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

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

  /* Read in the problem size */
#ifdef _WIN32
  scanf_s("%" NAG_IFMT "%*[^\n] ", &n);
#else
  scanf("%" NAG_IFMT "%*[^\n] ", &n);
#endif

  pdg = n;
  pdx = n;
  g = new DCO_TYPE [n * n];
  g_in = new DCO_TYPE [n * n];
  x = new DCO_TYPE [n * n];
  xr = new double [n * n];
  dxdg = new double [n * n];

  /* Read in the matrix g */
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
#ifdef _WIN32
      scanf_s("%lf", &dco::value(g_in[i+j*n]));
#else
      scanf("%lf", &dco::value(g_in[i+j*n]));
#endif
    }
  }
#ifdef _WIN32
  scanf_s("%*[^\n] ");
#else
  scanf("%*[^\n] ");
#endif

  /* Set up method parameters */
  errtol = 1.00e-7;
  maxits = 200;
  maxit = 10;


  for (i = 0; i < n * n; i ++) {
    DCO_MODE::global_tape->register_variable(g_in[i]);
    g[i] = g_in[i];
  }

  g02aa_a1w_f_(ad_handle, g, pdg, n, errtol, maxits, maxit, x, pdx,
                          iter, feval, nrmgrd, ifail);
  if (ifail != 0) {
    printf("Error from G02AA_A1W_F.%" NAG_IFMT "\n",
           ifail);
    exit_status = 1;
    goto END;
  }
  for (i=0; i < n*n; ++i) {
    xr[i] = dco::value(x[i]);
  }
  NagError  fail;
  INIT_FAIL(fail);
  x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,n,n,xr,n,
         "     Nearest Correlation Matrix",0,&fail);

  printf("\n\n Number of Newton steps taken:%11" NAG_IFMT "\n", iter);
  printf(" Number of function evaluations:%9" NAG_IFMT "\n", feval);
  if (nrmgrd > errtol)
    printf(" Norm of gradient of last Newton step: %12.3e\n",
           dco::value(nrmgrd));

  for (i = 0; i < n*n; i++) {
    dco::derivative(x[i]) += 1.0;
  }

  DCO_MODE::global_tape->interpret_adjoint();

  for (i = 0; i < n*n; ++i) {
    dxdg[i] = dco::derivative(g_in[i]);
  }
  printf("\n\n");
  x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,n,n,dxdg,n,
         "   (1,...,1)^T * dX/dG",0,&fail);
END:

  // Remove computational data object and tape
  x10ab_a1w_f_(ad_handle,ifail);
  DCO_TAPE_TYPE::remove(DCO_MODE::global_tape);

  delete [] g;
  delete [] g_in;
  delete [] x;
  delete [] xr;
  delete [] dxdg;

  return exit_status;
}