Logo MTL4

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.
//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// 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>
#include <boost/numeric/itl/itl_fwd.hpp>
#include <boost/numeric/mtl/operation/resource.hpp>

namespace itl {

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

  Scalar rho(0), rho_1(0), alpha(0);
  Vector p(resource(x)), q(resource(x)), r(resource(x)), z(resource(x));
  
  r = b - A*x;
  while (! iter.finished(r)) {
      ++iter;
      z = solve(L, r);
      rho = dot(r, z);
    
      if (iter.first())
          p = z;
      else 
          p = z + (rho / rho_1) * p;
      
      q = A * p;
      alpha = rho / dot(p, q);
      
      x += alpha * p;
      r -= alpha * q;
      rho_1 = rho;
  }
  return iter;
}

template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB, 
           typename Preconditioner, typename RightPreconditioner, typename Iteration >
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b, 
       const Preconditioner& L, const RightPreconditioner&, Iteration& iter)
{
    return cg(A, x, b, L, iter);
}

} // 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 -- Gen. with rev. 7542 on 7 Apr 2011 by doxygen 1.5.9 -- © 2010 by SimuNova UG.