/*
  Copyright (C) 2011 Carlo de Falco, http://mox.polimi.it/~carlo

  Copyright (C) 2006-2009 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru/

  The program reads the matrix and RHS in the Matrix Market format then run MUMPS.

  This software is a copyrighted work licensed under the terms, described
  in the file "FREE_LICENSE". 
*/

#define F77_COMM_WORLD -987654
#define JOB_INIT -1 
#define JOB_SOLVE 6
#define JOB_END  -2


#include <fstream>
#include <iomanip>

#include "matrices.h"

extern "C"
{
#include "dmumps_c.h"
}

#include <mpi.h>

int main (int argc, char **argv)
{	

  // 0) Init MPI
  MPI::Init (argc, argv);
  int 
    rank = MPI::COMM_WORLD.Get_rank (),
    size = MPI::COMM_WORLD.Get_size ();

  
  // 1) Read matrix on the host
  SparseMatrix mat;                             
  int IsSymmetric = 0;
  if (rank == 0)
    {
      ifstream in (argv[1]);
      if (!in)
        {
          cout << "cannot open file " << argv[1] << endl;
          MPI::Finalize ();
          return 2;
        }

      cout << "reading matrix from file: " << argv[1] << " ..." << endl;
      mat.read (in);
      cout << "... done!" << endl;
      
      // if you know that you matrix is positive definite, IsSymmetric should be 1, then it is faster
      IsSymmetric = mat.IsSymmetric ? 2 : 0;
    }

  MPI::COMM_WORLD.Bcast (&IsSymmetric, 1, MPI::INT, 0);

  // 1a) setting up MUMPS parameters and initialize it
  DMUMPS_STRUC_C id;
  id.job          = JOB_INIT; 
  id.par          = 1;  
  id.sym          = IsSymmetric;
  id.comm_fortran = F77_COMM_WORLD;
  dmumps_c(&id);

  //streams    
  id.icntl[0] = 6; 
  id.icntl[1] = 6; 
  id.icntl[2] = 6; 
  id.icntl[3] = 1;

  id.icntl[21] = 1;

  //ordering AMD (0), or metis (5), or pord (4), or AMF (2), or QAMD (6)	
  id.icntl[6] = 0;

  // 2) converting matrix to the MUMPS representation 
  
  vector<int> irn;
  vector<int> jcn;
  vector<double> a;

  if (rank == 0)
    {      
      id.n = mat.NRows ();
      id.nz = mat.nnz;
      irn.resize (id.nz);
      jcn.resize (id.nz);
      a.resize (id.nz);
      int k = 0;
      for (int i = 0; i < id.n; ++i)
        for (map<int, double>::iterator j = mat[i].begin (); j != mat[i].end (); ++j)
          {
            irn[k] = i + 1;
            jcn[k] = (*j).first + 1;
            a[k] = (*j).second;
            ++k;
          }
      id.irn = &*irn.begin (); 
      id.jcn = &*jcn.begin ();
      id.a   = &*a.begin ();
    } 
  mat.clearMemory();
    
  
  // 3) Read RHS and set up MUMPS
  Matrix rhs;
  if (rank == 0)
    {
      ifstream in2(argv[2]);
      if (!in2)
        {
          cout << "cannot open file " << argv[2] << endl;
          return 4;
        }

      cout << "reading rhs from file: " << argv[2] << " ..." << endl;
      rhs.read (in2);
      cout << " done!" << endl;

      id.rhs = &*rhs.begin ();
      id.nrhs = rhs.NCols ();
      id.lrhs = id.n;
    }

  // 4) Solve system
  id.job = JOB_SOLVE;
  dmumps_c (&id);
  
  // 5) write the solution out
  //solution is in rhs (it is overwritten)

  if (rank == 0)
    {
      cout << "Saving solution in " << argv[2] << ".solution ..." << endl;
      ofstream outf ((string (argv[2]) + ".solution").c_str ());
      rhs.write(outf);
      cout << " done!" << endl;
    }

  // 6) clean up
  id.job = JOB_END; 
  dmumps_c (&id);

  MPI::Finalize ();
  return 0;	
}