//===========================================================================
//  CVS Information:                                                         
//                                                                           
//     $RCSfile: qmr.cpp,v $  $Revision: 1.3 $  $State: Exp $ 
//     $Author: llee $  $Date: 2001/10/18 15:44:32 $ 
//     $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: qmr.cpp,v $
// Revision 1.3  2001/10/18 15:44:32  llee
// compile and test for examples
//
// Revision 1.2  2001/07/05 22:28:57  llee1
// gcc 3.0 fix
//
// Revision 1.1  2000/07/26 21:49:26  llee1
// change file extension from .cc to .cpp
//
// Revision 1.1  2000/07/25 22:57:53  llee1
// *** empty log message ***
//
// Revision 1.2  2000/07/17 15:44:04  llee1
// *** empty log message ***
//
// Revision 1.1  1999/11/21 20:47:07  lums
// *** empty log message ***
//
//===========================================================================

#include <itl/interface/blas.h>
#include <itl/krylov/qmr.h>

/*
  In thsi example, we show how to use QMR algorithm, the output should be:

iteration 0: resid 2.23607
iteration 1: resid 0.920575
iteration 2: resid 0.798734
iteration 3: resid 0.635402
iteration 4: resid 0.723989
iteration 5: resid 7.84055e-14
finished! error code = 0
5 iterations
7.84055e-14 is actual final residual.
3.5064e-14 is actual relative tolerance achieved.
Relative tol: 1e-06  Absolute tol: 0
Residual 7.7393e-14
*/

using namespace itl;

typedef  double Type;

typedef sparse_matrix<Type> Matrix;

typedef std::vector<Type> Vector;

int main()
{
  int max_iter = 50;

  std::vector<double> val(14);
  std::vector<int> ind(14);
  std::vector<int> ptr(6);
  for (int i = 0; i<14;i++) val[i] = i+1;
//   Matrix A(5, 5);
//   A(0, 0) = 1.0;
  ptr[0] = 0;
  ind[0] = 0;
//   A(0, 1) = 2.0;
  ind[1] = 1;
//   A(0, 4) = 3.0;
  ind[2] = 4;
  ptr[1] = 3;
//   A(1, 0) = 4.0;
  ind[3] = 0;
//   A(1, 1) = 5.0;
  ind[4] = 1;
//   A(1, 3) = 6.0;
  ind[5] = 3;
  ptr[2] = 6;
//   A(2, 1) = 7.0;
  ind[6] = 1;
//   A(2, 2) = 8.0;
  ind[7] = 2;
//   A(2, 4) = 9.0;
  ind[8] = 4;
  ptr[3] = 9;
//   A(3, 2) = 10.;
  ind[9] = 2;
//   A(3, 3) = 11.;
  ind[10] = 3;
//   A(3, 4) = 12.;
  ind[11] = 4;
  ptr[4] = 12;
//   A(4, 2) = 13.; 
  ind[12] = 2;
//   A(4, 4) = 14.;
  ind[13] = 4;
  ptr[5] = 14;

  //begin
  Matrix A(5, 5, 14, &val[0], &ind[0], &ptr[0]);
  //end

  //begin
  Vector x(A.nrows(), Type(0));
  Vector b(A.ncols());
  for (Vector::iterator i=b.begin(); i!=b.end(); i++)
    *i = 1.;

  identity_preconditioner precond;

  //iteration
  noisy_iteration<double> iter(b, max_iter, 1e-6);
  //qmr algorithm
  qmr(A, x, b, precond, precond, iter);
  //bicgstab(A, x, b, precond, iter);
  //bicg(A, x, b, precond, iter);
  //cgs(A, x, b, precond, iter);
  //tfqmr(A, x, b, precond, precond, iter);
  //end

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

  std::cout << "Residual " << itl::two_norm(b1) << std::endl;
    
  return 0;
}