/* F11BD_P0W_F C++ Header Example Program.
 *
 * Copyright 2020 Numerical Algorithms Group.
 * Mark 27.1, 2020.
 */

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

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

  cout << "F11BD_P0W_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  alpha;
  cin >> n;
  cin >> m;
  cin >> alpha;

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

  b    = new double [n];
  x    = new double [n];
  work = new double [lwork];

  double bb, b1, a, c;

  b1 = 12.0;
  a  = 1.0;
  c  = 1.0;
  bb = b1 - 2.0;

  ifail = 0;
  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;

  // Initialize rthe solver
  Integer iterm = 2, maxitn = 800, monit = 0, lwreq = lwork;
  double  sigmax = 0.0, anorm;
  double  tol = 1.0e-10;
  ifail = 0;
  f11bd_p0w_f_(ad_handle,"RGMRES","P","2","N",iterm,n,m,tol,maxitn,
               anorm,sigmax,monit,lwreq,work,lwork,ifail,6,1,1,1);

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

  while (irevcm != 4) {
    ifail = 0;
    f11be_p0w_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;
        }
      }
    }
  } 

  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 << x[i] << "     ";
    cout.width(13);cout << b[i] << endl;
  }

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

  return exit_status;
}