/* G02AB_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 */
  Integer  exit_status = 0;
  DCO_TYPE alpha, errtol, nrmgrd;
  Integer  feval, i, iter, j, maxit, maxits, n, pdg, pdx;

  /* Arrays */
  char     nag_enum_arg[100];

  DCO_TYPE *g = 0, *w = 0, *x = 0, *g_in = 0;
  double   *xr = 0, *dxdg = 0;
  double   dx;

  Integer  ifail;

  /* Output preamble */
  printf("G02AB_T1W_F C++ Header Example Program Results\n\n");
  fflush(stdout);

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

  /* Skip heading in data file */
#ifdef _WIN32
  scanf_s("%*[^\n] ");
#else
  scanf("%*[^\n] ");
#endif
  /* Read in the problem size, opt and alpha */
#ifdef _WIN32
  scanf_s("%" NAG_IFMT "", &n);
#else
  scanf("%" NAG_IFMT "", &n);
#endif
#ifdef _WIN32
  scanf_s("%39s", nag_enum_arg, (unsigned)_countof(nag_enum_arg));
#else
  scanf("%39s", nag_enum_arg);
#endif

#ifdef _WIN32
  scanf_s("%lf%*[^\n]", &dco::value(alpha));
#else
  scanf("%lf%*[^\n]", &dco::value(alpha));
#endif

  pdg = n;
  pdx = n;

  g = new DCO_TYPE [n * n];
  g_in = new DCO_TYPE [n * n];
  w = new DCO_TYPE [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

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

  /* Use the defaults for errtol, maxits and maxit */
  errtol = 0.0;
  maxits = 0;
  maxit = 0;

  for (int i = 0; i < n * n; i ++) {
    dco::derivative(g_in[i]) = 0.0;
    g[i] = g_in[i];
  }

  NagError  fail;
  INIT_FAIL(fail);

  for (int k = 0; k < n; k++) {
    for (int l = 0; l < n; l++) {

      dco::derivative(g[k+l*n]) += 1.0;

      g02ab_t1w_f_(ad_handle, g, pdg, n, "B", alpha, w, errtol,
                   maxits, maxit, x, pdx, iter, feval,
                   nrmgrd, ifail, 100);
      if (ifail != 0) {
        printf("%" NAG_IFMT "\n", ifail);
        exit_status = 1;
        goto END;
      }

      if ( k == 0 && l == 0 ) {
        for (int i = 0; i < n * n; i ++)
          xr[i] = dco::value(x[i]);

        x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,n,n,xr,n,
               "     Nearest Correlation Matrix",0,&fail);

        printf("\nNumber of Newton steps taken: %11" NAG_IFMT "\n", iter);
        printf("Number of function evaluations: %9" NAG_IFMT "\n\n", feval);
        printf("alpha: %37.3f \n\n", dco::value(alpha));
        fflush(stdout);
      }

      dx = 0.0;
      for (i = 0; i < n*n; i++) {
          dx += dco::derivative(x[i]);
      }
      dxdg[k+l*n] = dx;

      for (int i = 0; i < n * n; i++)
        g[i] = 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
  x10ab_t1w_f_(ad_handle,ifail);

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

  return exit_status;
}