//===========================================================================
//  CVS Information:                                                         
//                                                                           
//     $RCSfile: gmres.cpp,v $  $Revision: 1.2 $  $State: Exp $ 
//     $Author: llee $  $Date: 2001/10/18 16:57:11 $ 
//     $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
// Matrix 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.2  2001/10/18 16:57:11  llee
// *** empty log message ***
//
// Revision 1.1  2000/07/26 21:49:44  llee1
// change file extension from .cc to .cpp
//
// Revision 1.2  2000/07/19 19:58:00  llee1
// *** empty log message ***
//
// Revision 1.1  2000/07/17 15:44:06  llee1
// *** empty log message ***
//
// Revision 1.1  2000/03/19 21:17:37  llee1
// *** empty log message ***
//
// Revision 1.3  1999/11/30 19:50:00  llee1
// *** empty log message ***
//
// Revision 1.2  1999/11/29 01:41:02  llee1
// *** empty log message ***
//
// Revision 1.1  1999/11/29 00:44:26  llee1
// *** empty log message ***
//
//                                                                           
//===========================================================================

#if defined ITL_USING_MTL
#include "mtl/matrix.h"
#include "mtl/mtl.h"
#include "mtl/utils.h"

#include "itl/interface/mtl.h"
#else
#include "my_interface.h"
#endif

#include "itl/krylov/gmres.h"
#include <vector>
#include <iostream>
/*
  In thsi example, we show how to use GMRES algorithm, the output should be:

  teration 0: resid 1
  iteration 1: resid 0.480381
  iteration 2: resid 0.341156
  iteration 3: resid 0.234932
  iteration 4: resid 0.118403
  iteration 5: resid 8.04725e-16
  iteration 5: resid 6.18013e-16
  finished! error code = 0
  5 iterations
  6.18013e-16 final residual
    10      2      0      0      3   x           9.61541  =     100
     2     50      7      6      0   x           1.36768  =     100
     0      7     80      0      9   x           1.08868  =     100
     0      6      0    110     12   x          0.794106  =     100
     3      0      9     12    140   x          0.370188  =     100
*/
using namespace itl;
using namespace std;

typedef  double Type;

#if defined ITL_USING_MTL
typedef mtl::matrix< Type, mtl::rectangle<>,
                mtl::array< mtl::compressed<> >,
                mtl::row_major >::type Matrix;
typedef mtl::dense1D<Type> Vector;
#else
typedef std::vector<std::vector<Type> > Matrix;
typedef std::vector<Type> Vector;
#endif

int main()
{
  int max_iter = 50;
  //begin
#if defined ITL_USING_MTL
  Matrix A(5, 5);
  // Fill matrix...
  //end
  A(0, 0) = 10.0;
  A(0, 1) = 2.0;
  A(0, 4) = 3.0;
  A(1, 0) = 2.0;
  A(1, 1) = 50.0;
  A(1, 2) = 7.0;
  A(1, 3) = 6.0;
  A(2, 1) = 7.0;
  A(2, 2) = 80.0;
  A(2, 4) = 9.0;
  A(3, 1) = 6.;
  A(3, 3) = 110.;
  A(3, 4) = 12.;
  A(4, 0) = 3.0;
  A(4, 2) = 9.0;
  A(4, 3) = 12.; 
  A(4, 4) = 140.;
#else
  Matrix A(5, Vector(5, 0));
  A[0][0] = 10.;
  A[0][1] = 2.0;
  A[0][4] = 3.0;
  A[1][0] = 2.0;
  A[1][1] = 50.0;
  A[1][2] = 7.0;
  A[1][3] = 6.0;
  A[2][1] = 7.0;
  A[2][2] = 80.0;
  A[2][4] = 9.0;
  A[3][1] = 6.0;
  A[3][3] = 110.;
  A[3][4] = 12.;
  A[4][0] = 3; 
  A[4][2] = 9;
  A[4][3] = 12; 
  A[4][4] = 140.;
#endif
  //begin
  Vector x(nrows(A), Type(0));
  Vector b(nrows(A));
  for (Vector::iterator i=b.begin(); i!=b.end(); i++)
    *i = 100.;

  //SSOR preconditioner
  //SSOR<Matrix> precond(A);
  identity_preconditioner precond;

  Vector b2(b.size(), 0);
  solve(precond, b, b2);

  //iteration
  noisy_iteration<double> iter(b2, max_iter, 1e-12);
  int restart = 20;
  modified_gram_schmidt<Vector> orth(restart, size(x));

  gmres(A, x, b, precond(), restart,  iter, orth);
  //end

  //verify the result
  Vector b1(nrows(A));
  itl::mult(A, x, b1);
  itl::add(b1, scaled(b, -1.), b1);

  if ( itl::two_norm(b1) < 1.e-6) { //output
    for (int i=0; i<5; i++) {
      for (int j=0; j<5; j++) {
	cout.width(6);
#if defined ITL_USING_MTL
	cout << A(i,j) << " ";
#else
	cout << A[i][j] << " ";
#endif
      }
      cout << "  x  ";
      cout.width(16);
      cout << x[i] << "  =  ";
      cout.width(6);
      cout << b[i] << endl;
    }
  } else {
    cout << "Residual " << iter.resid() << endl;
  }

  return 0;
  }