//===========================================================================
//  CVS Information:                                                         
//                                                                           
//     $RCSfile: gmres.cpp,v $  $Revision: 1.7 $  $State: Exp $ 
//     $Author: llee $  $Date: 2002/06/15 15:49:47 $ 
//     $Locker:  $ 
//---------------------------------------------------------------------------
//                                                                           
// DESCRIPTION                                                               
//                                                                           
//---------------------------------------------------------------------------
//                                                                           
// LICENSE AGREEMENT                                                         
// Copyright 1997, University of Notre Dame.
// Authors: Andrew Lumsdaine, Lie-Quan Lee
//
// This file is part of the Iterative Template Library
//
// You should have received a copy of the License Agreement for the
// Iterative Template Library along with the software;  see the
// file LICENSE.  If not, contact Office of Research, University of Notre
// Dame, Notre Dame, IN  46556.
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//---------------------------------------------------------------------------
//                                                                           
// REVISION HISTORY:                                                         
//                                                                           
// $Log: gmres.cpp,v $
// Revision 1.7  2002/06/15 15:49:47  llee
// *** empty log message ***
//
// Revision 1.6  2002/06/04 03:02:34  llee
// *** empty log message ***
//


#include "mtl/matrix.h"
#include "mtl/mtl.h"
#include "mtl/utils.h"

#include "itl/interface/mtl.h"
//#include "itl/interface/blas.h"
#include "itl/preconditioner/ilu.h"
#include "itl/krylov/gmres.h"

#include "laplacian.h"
#include "timer.hpp"
#include "parser.h"

using namespace mtl;
using namespace itl;

typedef  double Type;

//begin
typedef matrix< Type, 
                rectangle<>, 
  compressed<>,
  //	        array< compressed<> >, 
                row_major >::type Matrix;
typedef matrix< Type, 
                rectangle<>, 
  	        array< compressed<> >, 
                row_major >::type MMatrix;


int main (int argc, char* argv[]) 
{
  using std::cout;
  using std::endl;
  using std::parser;

  if ( argc == 1 ) {
    cout << argv[0] << " --help will get your the usage of flags." << endl;
  }

  parser myparser(argc, argv);

  myparser.register_flag("--mx", 1, "number of points in x axis in the mesh");
  myparser.register_flag("--my", 1, "number of points in y axis in the mesh");
  myparser.register_flag("--max_it", 1, "maxmal iteration allowed");  
  myparser.register_flag("--ksp_atol", 1, "Absolute tolerance in KSP iteration");  
  myparser.register_flag("--ksp_rtol", 1, "Relative tolerance in KSP iteration");  

  myparser.help();
  int mx, my, max_it;
  double ksp_atol, ksp_rtol;
 
  myparser.get("--mx", mx, 16);
  myparser.get("--my", my, mx);
  myparser.get("--max_it", max_it, 4000);
  myparser.get("--ksp_atol", ksp_atol, 1.0e-12);
  myparser.get("--ksp_rtol", ksp_rtol, 1.0e-6);
 
#ifdef TEST_3D
  const int N = mx * my * mx;

  Matrix A(N, N);
  {
    MMatrix AA(N, N);
    generate_laplacian_3D(AA, mx);
    mtl::copy(AA, A);
  }
#else
  const int N = mx * my;

  Matrix A(N, N);
  {
    MMatrix AA(N, N);
    generate_laplacian_2D(AA, mx, my);
    mtl::copy(AA, A);
  }
#endif

  dense1D<Type> x(A.nrows(), Type(1));
  dense1D<Type> b(A.ncols(), 0);

  itl::mult(A, x, b);

  for (dense1D<Type>::iterator i=x.begin(); i!=x.end(); i++)
    *i = 0;

  //SSOR<Matrix> precond(A);
  //SSOR<Matrix>::Precond p = precond();
  ILU<Matrix> precond(A);
  ILU<Matrix>::Precond p = precond();

  //gmres needs the preconditioned b to pass into iter object.
  dense1D<Type> b2(A.ncols());
  solve(p, b, b2); 

  basic_iteration<double> iter(b2, max_it, ksp_rtol, ksp_atol);

  int restart = 20;
  modified_gram_schmidt<dense1D<Type> > orth(restart, size(x));
  //classical_gram_schmidt< dense1D<Type> > orth(restart, size(x));
  //gmres algorithm
  Timer timer;
  timer.start();
  gmres(A, x, b, p, restart, iter, orth); 
  timer.stop();
  //end

  //verify the result
  dense1D<Type> b1(A.ncols());
  itl::mult(A, x, b1);
  itl::add(b1, itl::scaled(b, -1.), b1);
  cout << "number of iterations: " << iter.iterations() 
       << " Residual: " << iter.resid() << endl;
  cout << " True Residual: " << itl::two_norm(b1) << endl;
  cout << "Time used: " << timer.clock_time() <<endl;
  return 0;
}