Introduction Krylov-Subspace Methods

The natural notation in MTL4 allows you to write Krylov-Subspace methods in the same way as in the mathematical literature. For instance, consider the conjugate gradient method as it is realized in the ITL version that is in the process of revision:

// Software License for MTL
// 
// Copyright (c) 2007 The Trustees of Indiana University.
//               2008 Dresden University of Technology and the Trustees of Indiana University.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
// 
// This file is part of the Matrix Template Library
// 
// See also license.mtl.txt in the distribution.

#ifndef ITL_CG_INCLUDE
#define ITL_CG_INCLUDE

#include <boost/numeric/mtl/concept/collection.hpp>


namespace itl {

template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB, 
           class Preconditioner, class Iteration >
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b, 
       const Preconditioner& M, Iteration& iter)
{
  typedef HilbertSpaceX TmpVec;
  typedef typename mtl::Collection<HilbertSpaceX>::value_type Scalar;

  Scalar rho, rho_1, alpha, beta;
  TmpVec p(size(x)), q(size(x)), r(size(x)), z(size(x));
  
  r = b - A*x;

  while (! iter.finished(r)) {
      z = solve(M, r);
      rho = dot(r, z);
    
      if (iter.first())
          p = z;
      else {
          beta = rho / rho_1;
          p = z + beta * p;
      }
      
      q = A * p;
      alpha = rho / dot(p, q);
      
      x += alpha * p;
      r -= alpha * q;
      rho_1 = rho;
      
      ++iter;
  }
  return iter.error_code();
}


} // namespace itl

#endif // ITL_CG_INCLUDE

If this iterative computation is performed with MTL4 operations on according objects the single statements are evaluated with expression templates providing equivalent performance as with algorithm-adapted loops. For a system with a million unknowns and a five-point-stencil as matrix (explicitly stored) about 10 iterations with a simple preconditioner can be performed in a second on a commodity PC.

Return to Triangular Solvers                                Table of Content                                Proceed to Using Predefined Linear Solvers






Introduction Krylov-Subspace Methods -- MTL 4 -- Peter Gottschling and Andrew Lumsdaine -- Generated on 19 May 2009 by Doxygen 1.5.5 -- Copyright 2007 by the Trustees of Indiana University.