//===========================================================================
//  CVS Information:                                                         
//                                                                           
//     $RCSfile: iccg.cpp,v $  $Revision: 1.4 $  $State: Exp $ 
//     $Author: llee $  $Date: 2002/04/25 18:02:57 $ 
//     $Locker:  $ 
//---------------------------------------------------------------------------
//                                                                           
// DESCRIPTION                                                               
//                                                                           
//---------------------------------------------------------------------------
//                                                                           
// LICENSE AGREEMENT                                                         
// Copyright 1997, University of Notre Dame.
// Authors: Andrew Lumsdaine, Lie-Quan Lee, and 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: iccg.cpp,v $
// Revision 1.4  2002/04/25 18:02:57  llee
// *** empty log message ***
//
// Revision 1.3  2001/10/18 14:08:31  llee
// re-organize the directory structures
//
// Revision 1.2  2000/07/27 04:39:18  llee1
// *** empty log message ***
//
// Revision 1.1  2000/07/26 21:49:19  llee1
// change file extension from .cc to .cpp
//
// Revision 1.3  2000/07/26 21:34:12  llee1
// *** empty log message ***
//
// Revision 1.2  2000/07/17 15:44:04  llee1
// *** empty log message ***
//
//                                                                           
//===========================================================================
#include "mtl/matrix.h"
#include "mtl/mtl.h"
#include "mtl/utils.h"

#include <itl/interface/mtl.h>
#include "itl/preconditioner/modified_cholesky.h"
#include "itl/krylov/cg.h"

#include <complex>

/*
iteration 0: resid 1
iteration 1: resid 0.824265
iteration 2: resid 2.29125
iteration 3: resid 0.0475466
iteration 4: resid 0.469612
iteration 5: resid 9.22858e-14
finished! error code = 0
5 iterations
9.22858e-14 final residual
     1      2      0      0      3   x          0.355348  =       1
     2      5      7      6      0   x          0.184494  =       1
     0      7      8     10      9   x          0.017229  =       1
     0      6     10     11     12   x         -0.125628  =       1
     3      0      9     12     14   x          0.091888  =       1
*/

using namespace mtl;
using namespace itl;
typedef  double Type;
//begin
typedef matrix<Type, 
               symmetric<lower>, 
               array< compressed<> >, 
               column_major>::type Matrix;
//end
int main () 
{
  using std::cout;
  using std::endl;

  int max_iter = 50;
  //begin
  Matrix A(5);
  //end
  A(0, 0) = 10.0;
  A(1, 0) = 2.0;
  A(4, 0) = 3.0;
  A(1, 1) = 20.0;
  A(3, 1) = 6.0;
  A(2, 1) = 7.0;
  A(2, 2) = 30.0;
  A(4, 2) = 9.0;
  A(3, 2) = 10.;
  A(3, 3) = 40.;
  A(4, 3) = 12.;
  A(4, 4) = 50.;
  A(0, 1) = 2.0;
  A(0, 4) = 3.0;
  A(1, 3) = 6.0;
  A(1, 2) = 7.0;
  A(2, 4) = 9.0;
  A(2, 3) = 10.;
  A(3, 4) = 12.;
  //begin
  dense1D<Type> x(A.nrows(), Type(0));
  dense1D<Type> b(A.ncols());
  for (dense1D<Type>::iterator i=b.begin(); i!=b.end(); i++)
    *i = 1.;
  //inomplete cholesky preconditioner
  modified_cholesky<Matrix> precond(A);
  noisy_iteration<double> iter(b, max_iter, 1e-16);
  cg(A, x, b, precond(), iter);
  //end

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

  for (int ii=0; ii<5; ii++) {
    for (int j=0; j<5; j++) {
      cout.width(6);
      //!!! it is not recommended to use A(i,j) for large sparse matrices
      cout << A(ii, j) << " ";
    }
    cout << "  x  ";
    cout.width(16);
      cout << x[ii] << "  =  ";
      cout.width(6);
      cout << b[ii] << endl;
  }
  
  cout << "True Residual: " << itl::two_norm(b1) << endl;
  return 0;
  }