项目作者: xiaohongchen1991

项目描述 :
A simple C++ library of Krylov subspace methods for solving linear systems
高级语言: C++
项目地址: git://github.com/xiaohongchen1991/krylov-solvers.git
创建时间: 2019-05-15T02:09:45Z
项目社区:https://github.com/xiaohongchen1991/krylov-solvers

开源协议:MIT License

下载


Release versions

version 1.0.0:

The library includes two solvers, gmres and conjudate gradient methods, and two preconditioners, identity and Jocobi preconditioners.

Quick start

This is a simple C++ library of some commonly used Krylov subspace methods for solving linear systems. The algorithms and APIs are designed to support boost.ublas containers.

Basics

  1. #include <gmres.hpp>
  2. namespace ublas=boost::numeric::ublas;

Create and initialize a ublas::matrix instance.

  1. using M = ublas::matrix<double>;
  2. M A(3, 3);
  3. A(0, 0) = 1.0;
  4. A(0, 1) = 2.0;
  5. A(0, 2) = 3.0;
  6. A(1, 0) = 4.0;
  7. A(1, 1) = 5.0;
  8. A(1, 2) = 6.0;
  9. A(2, 0) = 7.0;
  10. A(2, 1) = 8.0;
  11. A(2, 2) = 9.0;

Create and initialize a ublas::gmres instance.

  1. ublas::gmres<M> solver(A);

The function call operator is overloaded to solve the linear system Ax = b.

  1. using V = ublas::vector<double>;
  2. V b(3);
  3. b(0) = 14.0;
  4. b(1) = 32.0;
  5. b(2) = 50.0;
  6. V x(3, 0.0);
  7. auto result = solver(b, x);

To get the number of iterations and error, do

  1. int num_iter = std::get<0>(result);
  2. double error = std::get<1>(result);

Sparse matrix

The krylov-solvers library also supports boost.ublas sparse matrix containers.

  1. #include <conjugate_gradient.hpp>
  2. #include <boost/numeric/ublas/matrix_sparse.hpp>
  3. namespace ublas=boost::numeric::ublas;

Create a sparse matrix.

  1. using M = ublas::compressed_matrix<double>;
  2. M A(3, 3);
  3. A(0, 0) = -2.0;
  4. A(0, 1) = 1.0;
  5. A(1, 0) = 1.0;
  6. A(1, 1) = -2.0;
  7. A(1, 2) = 1.0;
  8. A(2, 1) = 1.0;
  9. A(2, 2) = -2.0;

Create a ublas::conjugate_gradient instance.

  1. ublas::conjugate_gradient<M> solver(A);

Solve linear system Ax = b.

  1. using V = ublas::vector<double>;
  2. V b(3);
  3. b(0) = 0.0;
  4. b(1) = 0.0;
  5. b(2) = -4.0;
  6. V x(3, 0.0);
  7. solver(b, x);

Configure solver parameters

Each solver class defines a nested struct “param”. One way to configure solver is to pass a param struct into its constructor. e.g.

  1. typename ublas::gmres<M>::param p;
  2. p.tol = 1e-8;
  3. ublas::gmres<M> solver(A, p);

The solver parameters can also be set by doing:

  1. solver.set_param(p);

Preconditioner

Two simple preconditioners are provided in “preconditioner.hpp” file, namely, identity preconditioner and Jacobi preconditioner. By default, identity preconditioner is used. To explicitly specify the preconditioner, do:

  1. using M = ublas::compressed_matrix<double>;
  2. using P = ublas::jacobi_precond<M>;
  3. ublas::conjugate_gradient<M, P> solver(A);

Customization

The Krylov subspace methods are iterative methods and do not require forming the matrix A explicitly. So the user can also initialize the solver by a functor which computes Ax. In a similar way, one can also specify a customized preconditioner. e.g.

  1. // functor to compute Ax
  2. auto Op = [&A](const auto& v){return ublas::prod(A, v);};
  3. // customized preconditioner
  4. auto PInv = [](const auto& v)->decltype(auto){return v;};
  5. ublas::conjugate_gradient<decltype(Op), decltype(PInv)> solver(Op, PInv);

Linking

This is a header only library.

Requirements

The library requires a C++ compiler that supports C++14 and boost.ublas. To build the test, add boost include directory into system path.

TODO list

  • Add bicgstab solver
  • Add ilu0 preconditioner