CoDiPack  2.2.0
A Code Differentiation Package
SciComp TU Kaiserslautern
Loading...
Searching...
No Matches
eigenLinearSystem.hpp
1/*
2 * CoDiPack, a Code Differentiation Package
3 *
4 * Copyright (C) 2015-2024 Chair for Scientific Computing (SciComp), University of Kaiserslautern-Landau
5 * Homepage: http://scicomp.rptu.de
6 * Contact: Prof. Nicolas R. Gauger (codi@scicomp.uni-kl.de)
7 *
8 * Lead developers: Max Sagebaum, Johannes Blühdorn (SciComp, University of Kaiserslautern-Landau)
9 *
10 * This file is part of CoDiPack (http://scicomp.rptu.de/software/codi).
11 *
12 * CoDiPack is free software: you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * CoDiPack is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty
19 * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
20 *
21 * See the GNU General Public License for more details.
22 * You should have received a copy of the GNU
23 * General Public License along with CoDiPack.
24 * If not, see <http://www.gnu.org/licenses/>.
25 *
26 * For other licensing options please contact us.
27 *
28 * Authors:
29 * - SciComp, University of Kaiserslautern-Landau:
30 * - Max Sagebaum
31 * - Johannes Blühdorn
32 * - Former members:
33 * - Tim Albring
34 */
35
36#pragma once
37
38#include "../../../config.h"
39
40#if CODI_EnableEigen
41
42 #include <Eigen/Eigen>
43 #include <vector>
44
46 #include "linearSystemHandler.hpp"
47
49namespace codi {
50
52 template<typename T_Type, template<typename> class T_Matrix, template<typename> class T_Vector>
54 public:
55
56 using Type = CODI_DD(T_Type, CODI_DEFAULT_LHS_EXPRESSION);
57
58 using Matrix = CODI_DD(CODI_T(T_Matrix<Type>),
59 CODI_T(Eigen::Matrix<Type, 2, 2>));
60 using Vector = CODI_DD(CODI_T(T_Vector<Type>),
61 CODI_T(Eigen::Matrix<Type, 2, 1>));
62
63 using Real = typename Type::Real;
64 using Identifier = typename Type::Identifier;
65
66 using MatrixReal = CODI_DD(CODI_T(T_Matrix<Real>),
67 CODI_T(Eigen::Matrix<Real, 2, 2>));
68 using VectorReal = CODI_DD(CODI_T(T_Vector<Real>),
69 CODI_T(Eigen::Matrix<Real, 2, 1>));
70
71 using MatrixIdentifier = CODI_DD(CODI_T(T_Matrix<Identifier>),
72 CODI_T(Eigen::Matrix<Identifier, 2, 2>));
73 using VectorIdentifier = CODI_DD(CODI_T(T_Vector<Identifier>),
74 CODI_T(Eigen::Matrix<Identifier, 2, 1>));
75 };
76
79 template<typename T_Type, template<typename> class T_Matrix, template<typename> class T_Vector>
80 struct EigenLinearSystem : public LinearSystemInterface<EigenLinearSystemTypes<T_Type, T_Matrix, T_Vector>> {
81 public:
83
84 using Type = typename InterfaceTypes::Type;
85
86 using Matrix = typename InterfaceTypes::Matrix;
89 using Vector = typename InterfaceTypes::Vector;
92
93 using Index = typename Matrix::Index;
94
95 /*******************************************************************************/
98
100 template<typename M>
102 return new MatrixReal(mat->rows(), mat->cols());
103 }
104
106 template<typename M>
108 return new MatrixIdentifier(mat->rows(), mat->cols());
109 }
110
112 template<typename V>
114 return new VectorReal(vec->size());
115 }
116
118 template<typename V>
120 return new VectorIdentifier(vec->size());
121 }
122
125 delete A_v;
126 }
127
130 delete A_id;
131 }
132
135 delete vec_v;
136 }
137
140 delete vec_id;
141 }
142
144 /*******************************************************************************/
147
149 template<typename Func, typename MatrixA, typename MatrixB>
150 void iterateMatrix(Func func, MatrixA* matA, MatrixB* matB) {
151 Index rows = matA->rows();
152 Index cols = matA->cols();
153
154 for (int i = 0; i < rows; i += 1) {
155 for (int j = 0; j < cols; j += 1) {
156 func(matA->coeffRef(i, j), matB->coeffRef(i, j));
157 }
158 }
159 }
160
162 template<typename Func, typename MatrixA, typename MatrixB, typename MatrixC>
163 void iterateMatrix(Func func, MatrixA* matA, MatrixB* matB, MatrixC* matC) {
164 Index rows = matA->rows();
165 Index cols = matA->cols();
166
167 for (int i = 0; i < rows; i += 1) {
168 for (int j = 0; j < cols; j += 1) {
169 func(matA->coeffRef(i, j), matB->coeffRef(i, j), matC->coeffRef(i, j));
170 }
171 }
172 }
173
175 template<typename Func, typename VectorA, typename VectorB>
176 void iterateVector(Func func, VectorA* vecA, VectorB* vecB) {
177 Index size = vecA->size();
178
179 for (int i = 0; i < size; i += 1) {
180 func(vecA->coeffRef(i), vecB->coeffRef(i));
181 }
182 }
183
185 template<typename Func, typename VectorA, typename VectorB, typename VectorC>
186 void iterateVector(Func func, VectorA* vecA, VectorB* vecB, VectorC* vecC) {
187 Index size = vecA->size();
188
189 for (int i = 0; i < size; i += 1) {
190 func(vecA->coeffRef(i), vecB->coeffRef(i), vecC->coeffRef(i));
191 }
192 }
193
195 template<typename Func, typename VectorA, typename VectorB, typename VectorC, typename VectorD>
196 void iterateVector(Func func, VectorA* vecA, VectorB* vecB, VectorC* vecC, VectorD* vecD) {
197 Index size = vecA->size();
198
199 for (int i = 0; i < size; i += 1) {
200 func(vecA->coeffRef(i), vecB->coeffRef(i), vecC->coeffRef(i), vecD->coeffRef(i));
201 }
202 }
203
205 /*******************************************************************************/
208
210 void solveSystem(MatrixReal const* A, VectorReal const* b, VectorReal* x);
211
213 /*******************************************************************************/
216
219 MatrixReal* A_v_trans = new MatrixReal(A_v->rows(), A_v->cols());
220 *A_v_trans = A_v->transpose();
221
222 return A_v_trans;
223 }
224
226 template<typename Func>
227 void iterateDyadic(Func func, MatrixIdentifier* mat_id, VectorReal* x_v, VectorReal* b_b) {
228 Index rows = mat_id->rows();
229 Index cols = mat_id->cols();
230
231 for (int i = 0; i < rows; i += 1) {
232 for (int j = 0; j < cols; j += 1) {
233 func(mat_id->coeffRef(i, j), x_v->coeffRef(j), b_b->coeffRef(i));
234 }
235 }
236 }
237
239 /*******************************************************************************/
242
244 void subtractMultiply(VectorReal* t, VectorReal const* b_d, MatrixReal const* A_d, VectorReal const* x) {
245 *t = *b_d - *A_d * *x;
246 }
247
249 };
250
253 template<typename T_Type, template<typename> class T_Matrix, template<typename> class T_Vector>
254 struct SparseEigenLinearSystem : public EigenLinearSystem<T_Type, T_Matrix, T_Vector> {
255 public:
257
258 using Type = typename InterfaceTypes::Type;
259
266
267 using Real = typename Type::Real;
268 using Identifier = typename Type::Identifier;
269
270 using Index = typename Matrix::Index;
271
272 private:
273
275 template<typename R, typename T, typename M>
276 R* cloneMatrix(M* mat) {
277 R* r = new R(mat->rows(), mat->cols());
278
279 std::vector<Eigen::Triplet<T>> entries(mat->nonZeros());
280
281 Index outerSize = mat->outerSize();
282
283 for (int k = 0; k < outerSize; ++k) {
284 typename M::InnerIterator iter(*mat, k);
285 while (iter) {
286 entries.push_back(Eigen::Triplet<T>(iter.row(), iter.col(), T()));
287 ++iter;
288 }
289 }
290 r->setFromTriplets(entries.begin(), entries.end());
291
292 return r;
293 }
294
295 public:
296
297 /*******************************************************************************/
300
302 template<typename M>
304 return cloneMatrix<MatrixReal, Real>(mat);
305 }
306
308 template<typename M>
310 return cloneMatrix<MatrixIdentifier, Identifier>(mat);
311 }
312
314 /*******************************************************************************/
317
319 template<typename Func, typename MatrixA, typename MatrixB>
320 void iterateMatrix(Func func, MatrixA* matA, MatrixB* matB) {
321 Index outerSize = matA->outerSize();
322
323 for (int k = 0; k < outerSize; ++k) {
324 typename MatrixA::InnerIterator iterA(*matA, k);
325 typename MatrixB::InnerIterator iterB(*matB, k);
326 while (iterA) {
327 func(iterA.valueRef(), iterB.valueRef());
328
329 ++iterA;
330 ++iterB;
331 }
332 }
333 }
334
336 template<typename Func, typename MatrixA, typename MatrixB, typename MatrixC>
337 void iterateMatrix(Func func, MatrixA* matA, MatrixB* matB, MatrixC* matC) {
338 Index outerSize = matA->outerSize();
339
340 for (int k = 0; k < outerSize; ++k) {
341 typename MatrixA::InnerIterator iterA(*matA, k);
342 typename MatrixB::InnerIterator iterB(*matB, k);
343 typename MatrixC::InnerIterator iterC(*matC, k);
344 while (iterA) {
345 func(iterA.valueRef(), iterB.valueRef(), iterC.valueRef());
346
347 ++iterA;
348 ++iterB;
349 ++iterC;
350 }
351 }
352 }
353
355 /*******************************************************************************/
358
360 template<typename Func>
361 void iterateDyadic(Func func, MatrixIdentifier* mat_id, VectorReal* x_v, VectorReal* b_b) {
362 Index outerSize = mat_id->outerSize();
363
364 for (int k = 0; k < outerSize; ++k) {
365 typename MatrixIdentifier::InnerIterator iterId(*mat_id, k);
366 while (iterId) {
367 func(iterId.valueRef(), x_v->coeffRef(iterId.col()), b_b->coeffRef(iterId.row()));
368
369 ++iterId;
370 }
371 }
372 }
373
375 };
376}
377
378#endif
#define CODI_DD(Type, Default)
Abbreviation for CODI_DECLARE_DEFAULT.
Definition macros.hpp:94
#define CODI_T(...)
Abbreviation for CODI_TEMPLATE.
Definition macros.hpp:111
CoDiPack - Code Differentiation Package.
Definition codi.hpp:90
Eigen definition for the LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:53
T_Type Type
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:56
T_Matrix< Real > MatrixReal
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:66
T_Vector< Type > Vector
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:60
typename Type::Real Real
See LhsExpressionInterface.
Definition eigenLinearSystem.hpp:63
T_Matrix< Identifier > MatrixIdentifier
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:71
T_Matrix< Type > Matrix
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:58
T_Vector< Real > VectorReal
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:68
typename Type::Identifier Identifier
See LhsExpressionInterface.
Definition eigenLinearSystem.hpp:64
T_Vector< Identifier > VectorIdentifier
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:73
Definition eigenLinearSystem.hpp:80
void solveSystem(MatrixReal const *A, VectorReal const *b, VectorReal *x)
Needs to be implemented by the user.
void iterateVector(Func func, VectorA *vecA, VectorB *vecB, VectorC *vecC)
Iterate over all elements in the vectors at the same time.
Definition eigenLinearSystem.hpp:186
typename InterfaceTypes::MatrixReal MatrixReal
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:87
void deleteMatrixReal(MatrixReal *A_v)
Delete a real matrix.
Definition eigenLinearSystem.hpp:124
typename Matrix::Index Index
Index of an Eigen matrix.
Definition eigenLinearSystem.hpp:93
void deleteMatrixIdentifier(MatrixIdentifier *A_id)
Delete an identifier matrix.
Definition eigenLinearSystem.hpp:129
typename InterfaceTypes::Matrix Matrix
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:86
typename InterfaceTypes::VectorIdentifier VectorIdentifier
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:91
void iterateVector(Func func, VectorA *vecA, VectorB *vecB, VectorC *vecC, VectorD *vecD)
Iterate over all elements in the vectors at the same time.
Definition eigenLinearSystem.hpp:196
void iterateVector(Func func, VectorA *vecA, VectorB *vecB)
Iterate over all elements in the vectors at the same time.
Definition eigenLinearSystem.hpp:176
MatrixReal * transposeMatrix(MatrixReal *A_v)
Create a transposed matrix.
Definition eigenLinearSystem.hpp:218
MatrixIdentifier * createMatrixIdentifier(M *mat)
Definition eigenLinearSystem.hpp:107
void deleteVectorReal(VectorReal *vec_v)
Delete a real Vector.
Definition eigenLinearSystem.hpp:134
typename InterfaceTypes::MatrixIdentifier MatrixIdentifier
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:88
typename InterfaceTypes::VectorReal VectorReal
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:90
typename InterfaceTypes::Type Type
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:84
void iterateDyadic(Func func, MatrixIdentifier *mat_id, VectorReal *x_v, VectorReal *b_b)
Definition eigenLinearSystem.hpp:227
void iterateMatrix(Func func, MatrixA *matA, MatrixB *matB)
Iterate over all elements in the matrices at the same time.
Definition eigenLinearSystem.hpp:150
void iterateMatrix(Func func, MatrixA *matA, MatrixB *matB, MatrixC *matC)
Iterate over all elements in the matrices at the same time.
Definition eigenLinearSystem.hpp:163
VectorIdentifier * createVectorIdentifier(V *vec)
Definition eigenLinearSystem.hpp:119
void deleteVectorIdentifier(VectorIdentifier *vec_id)
Delete an identifier vector.
Definition eigenLinearSystem.hpp:139
void subtractMultiply(VectorReal *t, VectorReal const *b_d, MatrixReal const *A_d, VectorReal const *x)
Computes t = b_d - A_d * x.
Definition eigenLinearSystem.hpp:244
MatrixReal * createMatrixReal(M *mat)
Definition eigenLinearSystem.hpp:101
typename InterfaceTypes::Vector Vector
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:89
VectorReal * createVectorReal(V *vec)
Definition eigenLinearSystem.hpp:113
Definition linearSystemInterface.hpp:109
Definition eigenLinearSystem.hpp:254
typename InterfaceTypes::VectorIdentifier VectorIdentifier
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:265
void iterateDyadic(Func func, MatrixIdentifier *mat_id, VectorReal *x_v, VectorReal *b_b)
Definition eigenLinearSystem.hpp:361
void iterateMatrix(Func func, MatrixA *matA, MatrixB *matB, MatrixC *matC)
Iterate over all elements in the matrices at the same time.
Definition eigenLinearSystem.hpp:337
MatrixReal * createMatrixReal(M *mat)
Definition eigenLinearSystem.hpp:303
typename Type::Real Real
See LhsExpressionInterface.
Definition eigenLinearSystem.hpp:267
typename Type::Identifier Identifier
See LhsExpressionInterface.
Definition eigenLinearSystem.hpp:268
typename InterfaceTypes::MatrixIdentifier MatrixIdentifier
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:262
typename InterfaceTypes::VectorReal VectorReal
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:264
typename InterfaceTypes::Type Type
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:258
typename InterfaceTypes::Vector Vector
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:263
typename InterfaceTypes::Matrix Matrix
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:260
MatrixIdentifier * createMatrixIdentifier(M *mat)
Definition eigenLinearSystem.hpp:309
typename InterfaceTypes::MatrixReal MatrixReal
See LinearSystemInterfaceTypes.
Definition eigenLinearSystem.hpp:261
void iterateMatrix(Func func, MatrixA *matA, MatrixB *matB)
Iterate over all elements in the matrices at the same time.
Definition eigenLinearSystem.hpp:320
typename Matrix::Index Index
Index of an Eigen matrix.
Definition eigenLinearSystem.hpp:270