/* F08GA_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 */
double eerrbd, eps;
Integer exit_status = 0, i, j, n;
/* Arrays */
DCO_TYPE *ap = 0, *dummy = 0, *w = 0, *work = 0, *ap_in = 0;
double *wr = 0, *dwda = 0;
Integer ifail;
std::string uplo = "U";
#define AP(I, J) ap_in[J * (J - 1) / 2 + I - 1]
// Create AD tape
DCO_MODE::global_tape=DCO_TAPE_TYPE::create();
void *ad_handle = 0;
x10aa_a1w_f_(ad_handle,ifail);
printf("F08GA_A1W_F C++ Example Program Results\n\n");
/* Skip heading in data file */
#ifdef _WIN32
scanf_s("%*[^\n]");
#else
scanf("%*[^\n]");
#endif
#ifdef _WIN32
scanf_s("%" NAG_IFMT "%*[^\n]", &n);
#else
scanf("%" NAG_IFMT "%*[^\n]", &n);
#endif
/* Allocate memory */
ap = new DCO_TYPE [n * (n + 1) / 2];
ap_in = new DCO_TYPE [n * (n + 1) / 2];
dummy = new DCO_TYPE [1];
w = new DCO_TYPE [n];
work = new DCO_TYPE [3*n];
wr = new double [n];
dwda = new double [n*n];
/* Read the upper triangular part of the matrix A from data file */
for (i = 1; i <= n; ++i) {
for (j = i; j <= n; ++j) {
AP(i,j) = 0.0;
#ifdef _WIN32
scanf_s("%lf", &dco::value(AP(i, j)));
#else
scanf("%lf", &dco::value(AP(i, j)));
#endif
}
}
#ifdef _WIN32
scanf_s("%*[^\n]");
#else
scanf("%*[^\n]");
#endif
for (int i = 0; i < n; i ++) {
DCO_MODE::global_tape->register_variable(ap_in[(i*(i+3))/2]);
}
for (int i = 0; i < (n*(n+1))/2; i ++) {
ap[i] = ap_in[i];
}
f08ga_a1w_f_(ad_handle, "N", "U", n, ap, w, dummy, 1, work, ifail, 1, 1);
if (ifail != 0) {
printf("Error from F08GA_A1W_F.\n%" NAG_IFMT " ", ifail);
exit_status = 1;
goto END;
}
/* Print solution */
for (i = 0; i < n; ++i) {
wr[i] = dco::value(w[i]);
}
printf("Eigenvalues\n");
for (j = 0; j < n; ++j)
printf("%8.4f%s", wr[j], (j + 1) % 8 == 0 ? "\n" : " ");
printf("\n");
/* Get the machine precision, eps, using nag_machine_precision (X02AJC)
* and compute the approximate error bound for the computed eigenvalues.
* Note that for the 2-norm, ||A|| = max {|w[i]|, i=0..n-1}, and since
* the eigenvalues are in ascending order ||A|| = max( |w[0]|, |w[n-1]|).
*/
eps = x02ajf_();
eerrbd = eps * MAX(fabs(wr[0]), fabs(wr[n - 1]));
/* Print the approximate error bound for the eigenvalues */
printf("\nError estimate for the eigenvalues\n");
printf("%11.1e\n", eerrbd);
for (i = 0; i < n; i ++) {
dco::a1w::global_ir->zero_adjoints();
dco::derivative(w[i]) += 1.0;
DCO_MODE::global_tape->interpret_adjoint();
for (int j = 0; j < n; j++)
dwda[i+j*n] = dco::derivative(ap_in[j*(j+3)/2]);
}
printf("\nDerivatives of eigenvalues w.r.t. diagonals of A\n");
NagError fail;
INIT_FAIL(fail);
x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,n,n,dwda,n,
" dW_i/dA_jj",0,&fail);
END:
delete [] ap;
delete [] ap_in;
delete [] dummy;
delete [] w;
delete [] work;
delete [] wr;
delete [] dwda;
// Remove computational data object and tape
x10ab_a1w_f_(ad_handle,ifail);
DCO_TAPE_TYPE::remove(DCO_MODE::global_tape);
return exit_status;
}
#undef AP