Cadabra Computer algebra system for field theory problems
Linear.hh
Go to the documentation of this file.
1
2 #pragma once
3
4 #include <numeric>
5 #include <vector>
6 #include <boost/numeric/ublas/matrix.hpp>
7 #include <boost/numeric/ublas/vector.hpp>
8 #include <boost/numeric/ublas/matrix_proxy.hpp>
9 #include <boost/numeric/ublas/lu.hpp>
10 #include "Storage.hh"
11
12 namespace linear {
13  // bool gaussian_elimination(const std::vector<std::vector<multiplier_t> >&, const std::vector<multiplier_t>& );
15
16  // Class for solving equations of the form Ax = y for x. First call
17  // factorize with the matrix A and then call solve as many times as needed
18  template <typename T>
19  struct Solver
20  {
21  public:
22  using matrix_type = boost::numeric::ublas::matrix<T>;
23  using vector_type = boost::numeric::ublas::vector<T>;
24
25  // Initialise the solver with the matrix A
26  bool factorize(const matrix_type& A_);
27
28  // Solve for Ax = y
29  vector_type solve(const vector_type& y);
30
31  private:
33  boost::numeric::ublas::vector<size_t> P;
35  size_t N;
36  };
37
38  template <typename T>
40  {
41  assert(A_.size1() == A_.size2());
42  N = A_.size1();
43  A = A_;
44  P.resize(N);
45
46  // Bring swap and abs into namespace
47  using namespace std;
48  using namespace boost::numeric::ublas;
49
50  std::iota(P.begin(), P.end(), 0);
51
52  for (size_t i = 0; i < N; ++i) {
53  T maxA = 0;
54  size_t imax = i;
55  for (size_t k = i; k < N; ++k) {
56  T absA = abs(A(k, i));
57  if (absA > maxA) {
58  maxA = absA;
59  imax = k;
60  }
61  }
62
63  if (imax != i) {
64  swap(P(i), P(imax)); //pivoting P
65  swap(row(A, i), row(A, imax)); //pivoting rows of A
66  }
67
68  if (A(i, i) == 0)
69  return false;
70
71  for (size_t j = i + 1; j < N; ++j) {
72  A(j, i) /= A(i, i);
73  for (size_t k = i + 1; k < N; ++k)
74  A(j, k) -= A(j, i) * A(i, k);
75  }
76  }
77
78  return true;
79  }
80
81
82  template <typename T>
84  {
85  x.resize(y.size());
86  for (size_t i = 0; i < N; ++i) {
87  x(i) = y(P(i));
88  for (size_t k = 0; k < i; ++k)
89  x(i) -= A(i, k) * x(k);
90  }
91
92  for (size_t i = N; i > 0; --i) {
93  for (size_t k = i; k < N; ++k)
94  x(i - 1) -= A(i - 1, k) * x(k);
95  x(i - 1) = x(i - 1) / A(i - 1, i - 1);
96  }
97
98  return x;
99  }
100
101  }
102
Definition: Linear.hh:12
int k
Definition: passing.cc:4
Definition: Linear.hh:20
boost::numeric::ublas::vector< size_t > P
Definition: Linear.hh:33
boost::numeric::ublas::matrix< T > matrix_type
Definition: Linear.hh:22
boost::numeric::ublas::vector< T > vector_type
Definition: Linear.hh:23
matrix_type A
Definition: Linear.hh:32
size_t N
Definition: Linear.hh:35
bool factorize(const matrix_type &A_)
Definition: Linear.hh:39
vector_type solve(const vector_type &y)
Definition: Linear.hh:83
vector_type x
Definition: Linear.hh:34