OGS
UnifiedMatrixSetters.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14
15#ifdef USE_PETSC
16
17// Global PETScMatrix/PETScVector //////////////////////////////////////////
18
19#include <numeric>
20
23
24namespace MathLib
25{
26void setVector(PETScVector& v, std::initializer_list<double> values)
27{
28 std::vector<double> const vals(values);
29 std::vector<PETScVector::IndexType> idcs(vals.size());
30 std::iota(idcs.begin(), idcs.end(), 0);
31
32 v.set(idcs, vals);
33}
34
37 double const value)
38{
39 v.set(index, value); // TODO handle negative indices
40}
41
42void setMatrix(PETScMatrix& m, std::initializer_list<double> values)
43{
44 m.setZero();
45 addToMatrix(m, values);
46}
47
48void setMatrix(PETScMatrix& m, Eigen::MatrixXd const& tmp)
49{
50 using IndexType = PETScMatrix::IndexType;
51
52 auto const rows = tmp.rows();
53 auto const cols = tmp.cols();
54
55 assert(rows == m.getNumberOfRows() && cols == m.getNumberOfColumns());
56
57 m.setZero();
58 std::vector<IndexType> row_idcs(rows);
59 std::vector<IndexType> col_idcs(cols);
60
61 std::iota(row_idcs.begin(), row_idcs.end(), 0);
62 std::iota(col_idcs.begin(), col_idcs.end(), 0);
63
64 // PETSc wants row-major
65 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
66 tmp_ = tmp;
67
68 m.add(row_idcs, col_idcs, tmp_);
69}
70
71void addToMatrix(PETScMatrix& m, std::initializer_list<double> values)
72{
73 using IndexType = PETScMatrix::IndexType;
74
75 auto const rows = m.getNumberOfRows();
76 auto const cols = m.getNumberOfColumns();
77
78 assert((IndexType)values.size() == rows * cols);
79
80 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp(
81 rows, cols);
82
83 auto it = values.begin();
84 for (IndexType r = 0; r < rows; ++r)
85 {
86 for (IndexType c = 0; c < cols; ++c)
87 {
88 tmp(r, c) = *(it++);
89 }
90 }
91
92 std::vector<IndexType> row_idcs(rows);
93 std::vector<IndexType> col_idcs(cols);
94
95 std::iota(row_idcs.begin(), row_idcs.end(), 0);
96 std::iota(col_idcs.begin(), col_idcs.end(), 0);
97
98 m.add(row_idcs, col_idcs, tmp);
99}
100
101} // namespace MathLib
102
103#else
104
105// Sparse global EigenMatrix/EigenVector
106// //////////////////////////////////////////
107
110
111namespace MathLib
112{
113void setVector(EigenVector& v_, std::initializer_list<double> values)
114{
115 auto& v(v_.getRawVector());
116 assert((std::size_t)v.size() == values.size());
117 auto it = values.begin();
118 for (std::size_t i = 0; i < values.size(); ++i)
119 {
120 v[i] = *(it++);
121 }
122}
123
124void setVector(EigenVector& v,
125 MatrixVectorTraits<EigenVector>::Index const index,
126 double const value)
127{
128 v.getRawVector()[index] = value;
129}
130
131void setMatrix(EigenMatrix& m, std::initializer_list<double> values)
132{
133 auto const rows = m.getNumberOfRows();
134 auto const cols = m.getNumberOfColumns();
135
136 assert(static_cast<EigenMatrix::IndexType>(values.size()) == rows * cols);
137 Eigen::MatrixXd tmp(rows, cols);
138
139 auto it = values.begin();
140 for (GlobalIndexType r = 0; r < rows; ++r)
141 {
142 for (GlobalIndexType c = 0; c < cols; ++c)
143 {
144 tmp(r, c) = *(it++);
145 }
146 }
147
148 m.getRawMatrix() = tmp.sparseView();
149}
150
151void setMatrix(EigenMatrix& m, Eigen::MatrixXd const& tmp)
152{
153 m.getRawMatrix() = tmp.sparseView();
154}
155
156void addToMatrix(EigenMatrix& m, std::initializer_list<double> values)
157{
158 auto const rows = m.getNumberOfRows();
159 auto const cols = m.getNumberOfColumns();
160
161 assert(static_cast<EigenMatrix::IndexType>(values.size()) == rows * cols);
162 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp(
163 rows, cols);
164
165 auto it = values.begin();
166 for (GlobalIndexType r = 0; r < rows; ++r)
167 {
168 for (GlobalIndexType c = 0; c < cols; ++c)
169 {
170 tmp(r, c) = *(it++);
171 }
172 }
173
174 m.getRawMatrix() += tmp.sparseView();
175}
176
177} // namespace MathLib
178
179#endif
GlobalMatrix::IndexType GlobalIndexType
Declaration of class PETScMatrix, which provides an interface to PETSc matrix routines.
Declaration of class PETScVector, which provides an interface to PETSc vector routines.
RawMatrixType::Index IndexType
Definition: EigenMatrix.h:32
Wrapper class for PETSc matrix routines for matrix.
Definition: PETScMatrix.h:32
PetscInt getNumberOfColumns() const
Get the number of columns.
Definition: PETScMatrix.h:76
void setZero()
Set all entries to zero.
Definition: PETScMatrix.h:95
void add(const PetscInt i, const PetscInt j, const PetscScalar value)
Add value to a single entry.
Definition: PETScMatrix.h:124
PetscInt getNumberOfRows() const
Get the number of rows.
Definition: PETScMatrix.h:74
Wrapper class for PETSc vector.
Definition: PETScVector.h:33
void addToMatrix(PETScMatrix &m, std::initializer_list< double > values)
void setVector(PETScVector &v, std::initializer_list< double > values)
static const double r
static const double v
void setMatrix(PETScMatrix &m, std::initializer_list< double > values)