AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
Constraints.hpp
1#pragma once
2
3#include <algorithm>
4#include <vector>
5
6#include <boost/numeric/mtl/matrix/compressed2D.hpp>
7#include <boost/numeric/mtl/matrix/inserter.hpp>
8#include <boost/numeric/mtl/utility/property_map.hpp>
9#include <boost/numeric/mtl/utility/range_wrapper.hpp>
10
11#include <amdis/linearalgebra/Constraints.hpp>
12#include <amdis/linearalgebra/SymmetryStructure.hpp>
13#include <amdis/linearalgebra/mtl/MatrixBackend.hpp>
14#include <amdis/linearalgebra/mtl/VectorBackend.hpp>
15
16namespace AMDiS
17{
18 template <class T>
20 {
22 using Vector = MTLVector<T>;
23
24 template <class BitVector, class Associations>
25 static void periodicBC(Matrix& mat, Vector& sol, Vector& rhs, BitVector const& left, Associations const& left2right, bool setDiagonal = true)
26 {
27 SymmetryStructure const symmetry = mat.symmetry();
28 if (symmetry == SymmetryStructure::spd ||
29 symmetry == SymmetryStructure::symmetric ||
30 symmetry == SymmetryStructure::hermitian)
31 symmetricPeriodicBC(mat.matrix(), sol.vector(), rhs.vector(), left, left2right, setDiagonal);
32 else
33 unsymmetricPeriodicBC(mat.matrix(), sol.vector(), rhs.vector(), left, left2right, setDiagonal);
34 }
35
36
37 template <class Mat, class Vec, class BitVector, class Associations>
38 static void symmetricPeriodicBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& left, Associations const& left2right, bool setDiagonal = true)
39 {
40 error_exit("Not implemented");
41 }
42
43
44 template <class Value>
45 struct Triplet
46 {
47 std::size_t row, col;
48 Value value;
49 };
50
51 template <class Mat, class Vec, class BitVector, class Associations>
52 static void unsymmetricPeriodicBC(Mat& mat, Vec& sol, Vec& rhs, BitVector const& left, Associations const& left2right, bool setDiagonal = true)
53 {
54 std::vector<Triplet<typename Mat::value_type>> rowValues;
55 rowValues.reserve(left2right.size()*std::size_t(mat.nnz()/(0.9*num_rows(mat))));
56
57 // Define the property maps
58 // auto row = mtl::mat::row_map(mat);
59 auto col = mtl::mat::col_map(mat);
60 auto value = mtl::mat::value_map(mat);
61
62 // iterate over the matrix
63 std::size_t slotSize = 0;
64 for (auto r : mtl::rows_of(mat)) {
65 if (left[r.value()]) {
66 slotSize = std::max(slotSize, std::size_t(mat.nnz_local(r.value())));
67 std::size_t right = left2right.at(r.value());
68
69 for (auto i : mtl::nz_of(r)) {
70 rowValues.push_back({right,col(i),value(i)});
71 value(i, 0);
72 }
73 }
74 }
75
76 mtl::mat::inserter<Mat, mtl::update_plus<typename Mat::value_type> > ins(mat, 2*slotSize);
77 for (auto const& t : rowValues)
78 ins[t.row][t.col] += t.value;
79
80 for (std::size_t i = 0; i < mtl::size(left); ++i) {
81 if (left[i]) {
82 std::size_t j = left2right.at(i);
83 if (setDiagonal) {
84 ins[i][i] = 1;
85 ins[i][j] = -1;
86 }
87
88 rhs[j] += rhs[i];
89 rhs[i] = 0;
90
91 sol[j] = sol[i];
92 }
93 }
94 }
95
96 };
97
98} // end namespace AMDiS
The basic container that stores a base matrix.
Definition: MatrixBackend.hpp:25
SymmetryStructure symmetry() const
Symmetry of the matrix entries.
Definition: MatrixBackend.hpp:244
BaseMatrix & matrix()
Return a reference to the data-matrix matrix.
Definition: MatrixBackend.hpp:49
The basic container that stores a base vector data.
Definition: VectorBackend.hpp:20
BaseVector const & vector() const
Return the data-vector vector_.
Definition: VectorBackend.hpp:43
Definition: Constraints.hpp:14