cholesky_external.cc

一个通用的数学库

//===========================================================================
//  CVS Information:                                                         
//                                                                           
//     $RCSfile: cholesky_external.cc,v $  $Revision: 1.2 $  $State: Exp $ 
//     $Author: jsiek $  $Date: 2000/08/28 15:22:58 $ 
//     $Locker:  $ 
//---------------------------------------------------------------------------
//                                                                           
// DESCRIPTION                                                               
//                                                                           
//---------------------------------------------------------------------------
//                                                                           
// LICENSE AGREEMENT                                                         
// Copyright 1997, University of Notre Dame.
// Authors: Andrew Lumsdaine, Jeremy G. Siek
//
// This file is part of the Matrix 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: cholesky_external.cc,v $
// Revision 1.2  2000/08/28 15:22:58  jsiek
// got rid of iostream.h and replaced with iostream
//
// Revision 1.1.1.1  2000/07/12 21:51:43  jsiek
// stable mtl version
//
// Revision 1.2  1999/09/30 07:10:08  jsiek
// argc, argv unused
//
// Revision 1.1  1999/07/22 21:17:10  jsiek
// blah
//
// Revision 1.7  1999/04/19 17:45:02  jsiek
// blah
//
// Revision 1.6  1999/03/10 17:07:57  jsiek
// sparse -> compressed<> and compressed -> compressed<> name changes
//
// Revision 1.5  1999/03/04 17:36:41  llee1
// *** empty log message ***
//
// Revision 1.4  1999/02/02 22:07:20  arodrig6
// added and updated BLAS level I and II examples to reflect new style vectors and Matrix definitions
//
// Revision 1.1  1998/08/17 19:18:01  llee1
// new file
//
//                                                                           
//===========================================================================
#include "mtl/matrix.h"
#include "mtl/mtl.h"
#include "mtl/utils.h"


#include "itl/itl.h"
#include "itl/cholesky.h"
#include "itl/cg.h"



/*
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<upper>, 
               compressed<int, external>, 
               row_major>::type Matrix;
//end
int main (int , char* []) 
{
  int max_iter = 50;
  //begin
  const int nnz = 19;
  int indices[nnz];
  double values[nnz];
  int ptr[6];

  indices[0] = 0;
  indices[1] = 1;
  indices[2] = 4;
  indices[3] = 0;
  indices[4] = 1;
  indices[5] = 2;
  indices[6] = 3;
  indices[7] = 1;
  indices[8] = 3;
  indices[9] = 2;
  indices[10] = 4;
  indices[11] = 1;
  indices[12] = 2;
  indices[13] = 3;
  indices[14] = 4;
  indices[15] = 0;
  indices[16] = 2;
  indices[17] = 3;
  indices[18] = 4;

  values[0] = 1;
  values[1] = 2;
  values[2] = 3;
  values[3] = 2;
  values[4] = 5;
  values[5] = 7;
  values[6] = 6;
  values[7] = 7;
  values[8] = 10;
  values[9] = 8;
  values[10] =9;
  values[11] =6;
  values[12] =10;
  values[13] =11 ;
  values[14] = 12;
  values[15] = 3;
  values[16] = 9;
  values[17] = 12;
  values[18] = 14;

  ptr[0] = 0;
  ptr[1] = 3;
  ptr[2] = 7;
  ptr[3] = 11;
  ptr[4] = 15;
  ptr[5] = 19;

  Matrix A(5, 5, nnz, values, ptr, indices);

  //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.;
  //incomplete cholesky preconditioner
  cholesky<Matrix> precond(A);
  noisy_iteration<double> iter(b, max_iter, 1e-6);
  cg(A, x, b, precond(), iter);
  //end

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

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

  return 0;
  }