/* F11DB_P0W_F C++ Header Example Program.
 *
 * Copyright 2020 Numerical Algorithms Group.
 * Mark 27.1, 2020.
 */
#include <dco_light.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 << "F11DB_P0W_F C++ Header Example Program Results\n\n";
  // Skip heading in data file
  string mystr;
  getline (cin, mystr);

  // Read order of matrix and number of nonzero entries
  Integer           n, nnz;
  cin >> n;
  cin >> nnz;

  Integer la = 2*nnz;
  Integer liwork = 7*n + 2;
  double  *a=0, *x=0, *y=0;
  Integer *icol=0, *idiag=0, *ipivp=0, *ipivq=0, *irow=0;
  Integer *istr=0, *iwork=0;

  a     = new double [la];
  x     = new double [n];
  y     = new double [n];
  icol  = new Integer [la];
  idiag = new Integer [n];
  ipivp = new Integer [n];
  ipivq = new Integer [n];
  irow  = new Integer [la];
  istr  = new Integer [n+1];
  iwork = new Integer [liwork];

  // Read the matrix A
  for (int i=0; i<nnz; i++) {
    cin >> a[i] >> irow[i] >> icol[i];
  }
      
  // Read the vector y
  for (int i=0; i<n; i++) {
    cin >> y[i];
  }


  ifail = 0;
  // Calculate LU factorization
  Integer lfill = -1;
  double  dtol;
  Integer nnzc, npivm;
  dtol = 0.0;
  ifail  = 0;
  f11da_p0w_f_(ad_handle,n,nnz,a,la,irow,icol,lfill,dtol,"C","N",
               ipivp,ipivq,istr,idiag,nnzc,npivm,iwork,liwork,ifail,1,1);

  // Check value of npivm
  if (npivm>0) {
    cout << " Factorization is not complete" << endl;
  } else {

    //       Solve P L D U x = y
    ifail = 0;
    f11db_p0w_f_(ad_handle,"N",n,a,la,irow,icol,ipivp,ipivq,istr,
                 idiag,"C",y,x,ifail,1,1);

    // Output results
    cout.setf(ios::scientific,ios::floatfield);
    cout.precision(4);
    cout << "  Solution vector" << endl;
    for (int i=0; i<n; ++i) {
      cout.width(12);cout << x[i] << "     ";
    }
  }
  
  delete [] a;
  delete [] x;
  delete [] y;
  delete [] icol;
  delete [] idiag;
  delete [] ipivp;
  delete [] ipivq;
  delete [] irow;
  delete [] istr;
  delete [] iwork;

  return exit_status;
}