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