/* nag::ad::g02ab Tangent Example Program.
*/
#include <dco.hpp>
#include <iostream>
#include <nagad.h>
std::stringstream filecontent("4 \
2.0 -1.0 0.0 0.0 \
-1.0 2.0 -1.0 0.0 \
0.0 -1.0 2.0 -1.0 \
0.0 0.0 -1.0 2.0 \
100.0 20.0 20.0 20.0");
// Function which calls NAG AD routines.
template <typename T>
void func(std::vector<T> &g, std::vector<T> &w, std::vector<T> &x);
// Driver with tangent calls.
// Computes the Nearest Correlation Matrix X for an input matrix G.
// Also, computes the sum of all Hessian elements of d2X/dG2
void driver(const std::vector<double> &gv,
const std::vector<double> &wv,
std::vector<double> & xv,
double & d2x_dg2);
int main()
{
std::cout << " nag::ad::g02ab Tangent over Tangent Example Program Results\n";
Integer n;
filecontent >> n;
// Input Matrix G, whose NCM we want to compute
std::vector<double> gv(n * n);
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < n; ++j)
{
filecontent >> gv[i + j * n];
}
}
// The square roots of the diagonal elements of W
std::vector<double> wv(n);
for (int i = 0; i < n; i++)
{
filecontent >> wv[i];
}
// Output NCM X
std::vector<double> xv(n * n);
// Sum of all Hessian elements dX/dG
double d2x_dg2;
driver(gv, wv, xv, d2x_dg2);
std::cout.setf(std::ios::scientific, std::ios::floatfield);
std::cout.precision(5);
std::cout << "\n Nearest Correlation Matrix:\n\n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
std::cout.width(13);
std::cout << xv[i + j * n];
}
std::cout << std::endl;
}
std::cout << "\n Derivatives calculated: Second order tangents\n";
// Print derivatives of NCM X w.r.t. input matrix G
std::cout << "\n Sum of Hessian elements of NCM X w.r.t. input matrix G:\n";
std::cout << " sum_ij [d2X/dG2]_ij = " << d2x_dg2 << std::endl;
return 0;
}
// Driver with tangent calls.
// Computes the Nearest Correlation Matrix X for an input matrix G.
// Also, computes the sum of all Hessian elements of d2X/dG2
void driver(const std::vector<double> &gv,
const std::vector<double> &wv,
std::vector<double> & xv,
double & d2x_dg2)
{
using T = dco::gt1s<dco::gt1s<double>::type>::type;
// AD data types
std::vector<T> g(gv.size());
dco::passive_value(g) = gv;
dco::derivative(dco::value(g)) = std::vector<double>(gv.size(), 1.0);
dco::value(dco::derivative(g)) = std::vector<double>(gv.size(), 1.0);
std::vector<T> w(wv.size());
dco::passive_value(w) = wv;
// Variable to differentiate
std::vector<T> x(xv.size());
// Call the NAG AD Lib functions
func(g, w, x);
// Extract the computed NCM
xv = dco::passive_value(x);
// Sum of the tangents of x
d2x_dg2 = 0.0;
for (int i = 0; i < x.size(); i++)
{
d2x_dg2 += dco::derivative(dco::derivative(x[i]));
}
}
// Function which calls NAG AD routines.
template <typename T>
void func(std::vector<T> &g, std::vector<T> &w, std::vector<T> &x)
{
Integer n = sqrt(g.size());
Integer pdg = n;
Integer pdx = n;
// Set up method parameters
T errtol = 0.0, alpha = 0.02;
Integer maxits = 0;
Integer maxit = 0;
// Output variables
Integer iter, feval;
T nrmgrd;
// Set a copy of g, since it will be modified by the example
std::vector<T> g_tmp = g;
// Create AD configuration data object
Integer ifail = 0;
nag::ad::handle_t ad_handle;
// Routine for computing the NCM of G.
nag::ad::g02ab(ad_handle, g_tmp.data(), pdg, n, "B", alpha, w.data(), errtol,
maxits, maxit, x.data(), pdx, iter, feval, nrmgrd, ifail);
}