/* E02BC_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>
using namespace std;
// typedef DCO_TYPE nagad_a1w_w_rtype
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 exit_status = 0;
cout << "E02BC_A1W_F C++ Header Example Program\n";
// Create AD handle
Integer ifail = 0;
void *ad_handle = 0;
x10aa_a1w_f_(ad_handle,ifail);
// Skip heading in data file
string mystr;
getline (cin, mystr);
Integer ncap, m;
cin >> ncap; cin >> m;
Integer ncap7 = ncap + 7;
if (m <= 0) {
printf("Invalid m.\n");
exit_status = 1;
return exit_status;
}
if (ncap <= 0) {
printf("Invalid ncap.\n");
exit_status = 1;
return exit_status;
}
// Allocate arrays
DCO_TYPE *c = 0, *lamda = 0;
double *dsdl = 0, *dsdc = 0;
c = new DCO_TYPE [ncap7];
lamda = new DCO_TYPE [ncap7];
dsdc = new double [2*m*ncap7];
dsdl = new double [2*m*ncap7];
// Read knots and spline coefficients from file
double tmp;
for (int j = 0; j < ncap7; j++) {
cin >> tmp;
lamda[j] = tmp;
}
for (int j = 0; j < ncap + 3; j++) {
cin >> tmp;
c[j] = tmp;
}
cout << "\n x Spline 1st deriv";
cout << " 2nd deriv 3rd deriv\n";
cout.setf(ios::scientific,ios::floatfield);
cout.setf(ios::right);
cout.precision(4);
DCO_TYPE s[4], x;
for (int i = 0; i < m; i++) {
cin >> tmp;
x = tmp;
for (Integer left = 1; left <= 2; left++) {
// Create AD tape
// nagad_a1w_ir_create();
DCO_MODE::global_tape=DCO_TAPE_TYPE::create();
for (int l = 0; l < ncap7; l++) {
// nagad_a1w_ir_register_variable(&lamda[l]);
DCO_MODE::global_tape->register_variable(lamda[l]);
}
for (int l = 0; l < ncap7; l++) {
// nagad_a1w_ir_register_variable(&c[l]);
DCO_MODE::global_tape->register_variable(c[l]);
}
e02bc_a1w_f_(ad_handle, ncap7, lamda, c, x, left, s, ifail);
if (ifail != 0) {
printf("\nError from e02bc_a1w_f_.\n%" NAG_IFMT "\n", ifail);
exit_status = 1;
return exit_status;
}
// double xv = x.value;
double xv = dco::value(x);
cout.width(11); cout << xv;
if (left == 1) {
cout << " Left";
} else {
cout << " Right";
}
for (int l = 0; l < 4; l++) {
// tmp = s[l].value;
tmp = dco::value(s[l]);
cout.width(12); cout << tmp;
}
// s[0].tangent = 1.0;
dco::derivative(s[0]) = 1.0;
// nagad_a1w_ir_interpret_adjoint(ifail);
DCO_MODE::global_tape->interpret_adjoint();
Integer indx = (left-1)*m*ncap7 + i*ncap7;
for (int l = 0; l < ncap7; l++) {
// tmp = lamda[l].tangent;
tmp = dco::derivative(lamda[l]);
dsdl[indx + l] = tmp;
// tmp = c[l].tangent;
tmp = dco::derivative(c[l]);
dsdc[indx + l] = tmp;
}
// nagad_a1w_ir_remove();
DCO_TAPE_TYPE::remove(DCO_MODE::global_tape);
cout << "\n";
}
cout << "\n";
}
cout << "\n Derivatives of spline values with respect to knots\n";
NagError fail;
INIT_FAIL(fail);
x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,ncap7,m,dsdl,ncap7,
"Left: ds/dlamda",0,&fail);
cout << "\n";
x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,ncap7,m,&dsdl[m*ncap7],
ncap7, "Right: ds/dlamda",0,&fail);
cout << "\n\n Derivatives with respect to spline coefficients\n";
x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,ncap7,m,dsdc,ncap7,
"Left: ds/dC",0,&fail);
cout << "\n";
x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,ncap7,m,&dsdc[m*ncap7],
ncap7, "Right: ds/dC",0,&fail);
x10ab_a1w_f_(ad_handle,ifail);
delete [] c;
delete [] lamda;
delete [] dsdc;
delete [] dsdl;
return exit_status;
}