/* G02AB_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)
{
/* 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;
Integer ifail;
// Create AD tape
DCO_MODE::global_tape=DCO_TAPE_TYPE::create();
void *ad_handle = 0;
x10aa_a1w_f_(ad_handle,ifail);
/* Output preamble */
printf("G02AB_A1W_F ");
printf("C++ Header Example Program Results\n\n");
fflush(stdout);
/* 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_MODE::global_tape->register_variable(g_in[i]);
g[i] = g_in[i];
}
g02ab_a1w_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;
}
for (int 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("\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);
for (i = 0; i < n*n; i++) {
dco::derivative(x[i]) += 1.0;
}
DCO_MODE::global_tape->interpret_adjoint();
for (int 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 [] w;
delete [] x;
delete [] g_in;
delete [] xr;
delete [] dxdg;
return exit_status;
}