/* F11BD_T1W_F C++ Header Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 * Mark 27, 2019.
 */

#include <dco.hpp>
#include <nag.h>
#include <nagx04.h>
#include <nagad.h>
#include <stdio.h>
#include <iostream>
using namespace std;

int main(void)
{
  int       exit_status = 0;
  void      *ad_handle = 0;
  Integer   ifail = 0;

  cout << "F11BD_T1W_F C++ Header Example Program Results\n\n";
  // Skip heading in data file
  string mystr;
  getline (cin, mystr);

  // Read problem size
  Integer n, m;
  double  alphar;
  cin >> n;
  cin >> m;
  cin >> alphar;

  // Allocate arrays containing A and its factorized form, B
  // and the solution X.
  nagad_t1w_w_rtype *b=0, *x=0, *work=0;
  double            *dx=0;
  Integer           lwork = 2*m*n + 1000;

  b    = new nagad_t1w_w_rtype [n];
  x    = new nagad_t1w_w_rtype [n];
  work = new nagad_t1w_w_rtype [lwork];
  dx   = new double [2*n];
  
  nagad_t1w_w_rtype alpha, b1, a, c;
  alpha = alphar;

  b1 = 12.0;
  a  = 1.0;
  c  = 1.0;

  // Create AD configuration data object
  ifail = 0;
  x10aa_t1w_f_(ad_handle,ifail);

  // Initialize rthe solver
  Integer           iterm = 2, maxitn = 800, monit = 0, lwreq = lwork;
  nagad_t1w_w_rtype sigmax = 0.0, anorm;
  nagad_t1w_w_rtype tol = 1.0e-10;
  double inc = 1.0, zero = 0.0;
  for (int k=0; k<2; ++k) {
    ifail = 0;
    f11bd_t1w_f_(ad_handle,"RGMRES","P","2","N",iterm,n,m,tol,maxitn,
                 anorm,sigmax,monit,lwreq,work,lwork,ifail,6,1,1,1);

    nagad_t1w_w_rtype bb = b1 - 2.0;
    if (k==0) {
      nagad_t1w_inc_derivative(&alpha,inc);
    } else {
      nagad_t1w_inc_derivative(&bb,inc);
    }

    for (int i=0; i<n; ++i) {
      b[i] = b1*(i+1);
      x[i] = 3.0;
    }
    b[n-1] = b[n-1] - (n+1);

    b[0] = b[0] + (b1-1.0)*alpha;
    for (int i=1; i<n-1; ++i) {
      b[i] = b[i] + b1*alpha;
    }
    b[n-1] = b[n-1] + (b1-1.0)*alpha;

    // Reverse communication call of solver
    Integer           irevcm = 0;
    nagad_t1w_w_rtype wgt[1];

    while (irevcm != 4) {
      ifail = 0;
      f11be_t1w_f_(ad_handle,irevcm,x,b,wgt,work,lwreq,ifail);
      if (irevcm != 4) {
        ifail = -1;
        if (irevcm == -1) {
          //  b = A^Tx
          b[0] = bb*x[0] + a*x[1];
          for (int i=1; i<n-1; ++i) {
            b[i] = c*x[i-1] + bb*x[i] + a*x[i+1];
          }
          b[n-1] = c*x[n-2] + bb*x[n-1];
        } 
        if (irevcm == 1) {
          // b = Ax
          b[0] = bb*x[0] + c*x[1];
          for (int i=1; i<n-1; ++i) {
            b[i] = a*x[i-1] + bb*x[i] + c*x[i+1];
          }
          b[n-1] = a*x[n-2] + bb*x[n-1];
        }
        if (irevcm == 2) {
          for (int i=0; i<n; ++i) {
            b[i] =  x[i]/bb;
          }
        }
      }
    } 
    if (k==0) {
      nagad_t1w_set_derivative(&alpha,zero);
      for (int j=0; j<n; j++) {
        dx[j] = nagad_t1w_get_derivative(x[j]);
      }
    } else {
      nagad_t1w_set_derivative(&bb,zero);
      for (int j=0; j<n; j++) {
        dx[n+j] = nagad_t1w_get_derivative(x[j]);
      }
    }
  }
  cout.setf(ios::scientific,ios::floatfield);
  cout.precision(2);

  cout << "  Solution vector   Residual vector\n";
  for (int i=0;i<n;++i) {
    cout.width(12);cout << nagad_t1w_get_value(x[i]) << "     ";
    cout.width(13);cout << nagad_t1w_get_value(b[i]) << endl;
  }

  cout << "\n\n Derivatives calculated: First order tangents\n";
  cout << " Computational mode    : algorithmic\n";
  cout << "\n Derivatives of X w.r.t alpha and bb:\n";
  
  // Print derivatives
  cout << endl;
  NagError fail;
  INIT_FAIL(fail);
  x04cac(Nag_ColMajor,Nag_GeneralMatrix,Nag_NonUnitDiag,n,2,dx,n,
         "      d/dalpha    d/dbb",0,&fail);

  // Remove computational data object
  ifail = 0;
  x10ab_t1w_f_(ad_handle,ifail);

  delete [] b;
  delete [] x;
  delete [] work;
  delete [] dx;

  return exit_status;
}