Skip to content

Commit 96115a6

Browse files
AlePaluAlePalu
committed
[unstable]: reimplementing testing suite, linear_algebra Matrix and Vector API tested
1 parent e85edd4 commit 96115a6

File tree

105 files changed

+584
-277376
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

105 files changed

+584
-277376
lines changed

‎fdaPDE/linear_algebra.h‎

Lines changed: 63 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -20,40 +20,40 @@
2020
// clang-format off
2121

2222
// include Eigen linear algebra library
23-
#include <Eigen/Eigen>
24-
#define __FDAPDE_HAS_EIGEN__
25-
26-
namespace fdapde {
27-
namespace internals {
28-
29-
// define basic eigen traits
30-
template <typename XprType> struct is_eigen_dense_xpr {
31-
static constexpr bool value =
32-
std::is_base_of<Eigen::MatrixBase<std::decay_t<XprType>>, std::decay_t<XprType>>::value;
33-
};
34-
template <typename XprType> constexpr bool is_eigen_dense_xpr_v = is_eigen_dense_xpr<XprType>::value;
35-
template <typename XprType> class is_eigen_dense_vec {
36-
private:
37-
using XprType_ = std::decay_t<XprType>;
38-
static constexpr bool check_() {
39-
if constexpr (is_eigen_dense_xpr_v<XprType_>) {
40-
return XprType_::IsVectorAtCompileTime;
41-
}
42-
return false;
43-
}
44-
public:
45-
static constexpr bool value = check_();
46-
};
47-
template <typename XprType> constexpr bool is_eigen_dense_vec_v = is_eigen_dense_vec<XprType>::value;
48-
49-
template <typename XprType> struct is_eigen_sparse_xpr {
50-
static constexpr bool value =
51-
std::is_base_of_v<Eigen::SparseMatrixBase<std::decay_t<XprType>>, std::decay_t<XprType>>;
52-
};
53-
template <typename XprType> constexpr bool is_eigen_sparse_xpr_v = is_eigen_sparse_xpr<XprType>::value;
54-
55-
} // namespace internals
56-
} // namespace fdapde
23+
// #include <Eigen/Eigen>
24+
// #define __FDAPDE_HAS_EIGEN__
25+
26+
// namespace fdapde {
27+
// namespace internals {
28+
29+
// // define basic eigen traits
30+
// template <typename XprType> struct is_eigen_dense_xpr {
31+
// static constexpr bool value =
32+
// std::is_base_of<Eigen::MatrixBase<std::decay_t<XprType>>, std::decay_t<XprType>>::value;
33+
// };
34+
// template <typename XprType> constexpr bool is_eigen_dense_xpr_v = is_eigen_dense_xpr<XprType>::value;
35+
// template <typename XprType> class is_eigen_dense_vec {
36+
// private:
37+
// using XprType_ = std::decay_t<XprType>;
38+
// static constexpr bool check_() {
39+
// if constexpr (is_eigen_dense_xpr_v<XprType_>) {
40+
// return XprType_::IsVectorAtCompileTime;
41+
// }
42+
// return false;
43+
// }
44+
// public:
45+
// static constexpr bool value = check_();
46+
// };
47+
// template <typename XprType> constexpr bool is_eigen_dense_vec_v = is_eigen_dense_vec<XprType>::value;
48+
49+
// template <typename XprType> struct is_eigen_sparse_xpr {
50+
// static constexpr bool value =
51+
// std::is_base_of_v<Eigen::SparseMatrixBase<std::decay_t<XprType>>, std::decay_t<XprType>>;
52+
// };
53+
// template <typename XprType> constexpr bool is_eigen_sparse_xpr_v = is_eigen_sparse_xpr<XprType>::value;
54+
55+
// } // namespace internals
56+
// } // namespace fdapde
5757

5858
// include required modules
5959
#include "utility.h"
@@ -80,12 +80,10 @@ namespace internals {
8080

8181
// detects whether XprType represents a static sized or dynamic sized expression
8282
template <typename XprType> struct is_dynamic_sized {
83-
private:
84-
using XprTypeClean = std::decay_t<XprType>;
85-
public:
86-
static constexpr bool value = XprTypeClean::Rows == Dynamic || XprTypeClean::Cols == Dynamic;
83+
static constexpr bool value = std::decay_t<XprType>::Rows == Dynamic || std::decay_t<XprType>::Cols == Dynamic;
8784
};
8885
template <typename XprType> static constexpr bool is_dynamic_sized_v = is_dynamic_sized<XprType>::value;
86+
8987

9088
// if XprType has its NestAsRef bit set, sets type member type to XprType&, otherwise just repeats XprType
9189
template <typename XprType, bool has_ref_bit> struct ref_select_impl;
@@ -105,37 +103,44 @@ template <typename XprType> using ref_select_t = typename ref_select<XprType>::t
105103
} // namespace internals
106104
}
107105

106+
// matrix types
107+
108+
#include "src/linear_algebra/unary_op.h"
109+
#include "src/linear_algebra/binary_op.h"
110+
#include "src/linear_algebra/vectorwise_op.h"
111+
#include "src/linear_algebra/block.h"
112+
113+
108114
#include "src/linear_algebra/matrix.h"
109115
#include "src/linear_algebra/diagonal.h"
110116
#include "src/linear_algebra/triangular.h"
111-
112-
// special matrices
113117
#include "src/linear_algebra/orthogonal.h"
114118
#include "src/linear_algebra/symmetric.h"
115-
#include "src/linear_algebra/skew.h"
116-
#include "src/linear_algebra/permutation.h"
117-
#include "src/linear_algebra/spd.h"
119+
// #include "src/linear_algebra/skew.h"
120+
// #include "src/linear_algebra/permutation.h"
121+
// #include "src/linear_algebra/spd.h"
118122

119-
#include "src/linear_algebra/matrix_expr.h"
123+
// expression template system
124+
#include "src/linear_algebra/expr.h"
120125

121126
// algorithms
122127
#include "src/linear_algebra/evd.h"
123-
#include "src/linear_algebra/partial_piv_lu.h"
124-
#include "src/linear_algebra/qr.h"
128+
// #include "src/linear_algebra/partial_piv_lu.h"
129+
// #include "src/linear_algebra/qr.h"
125130

126131
// eigen support
127-
#include "src/linear_algebra/eigen/utility.h"
128-
129-
#include "src/linear_algebra/eigen/eigen_helper.h"
130-
#include "src/linear_algebra/eigen/fspai.h"
131-
#include "src/linear_algebra/eigen/kronecker.h"
132-
#include "src/linear_algebra/eigen/lumping.h"
133-
#include "src/linear_algebra/eigen/sparse_block_matrix.h"
134-
#include "src/linear_algebra/eigen/woodbury.h"
135-
// randomized linear algebra
136-
#include "src/linear_algebra/eigen/rsi.h"
137-
#include "src/linear_algebra/eigen/rbki.h"
138-
#include "src/linear_algebra/eigen/rp_chol.h"
132+
// #include "src/linear_algebra/eigen/utility.h"
133+
134+
// #include "src/linear_algebra/eigen/eigen_helper.h"
135+
// #include "src/linear_algebra/eigen/fspai.h"
136+
// #include "src/linear_algebra/eigen/kronecker.h"
137+
// #include "src/linear_algebra/eigen/lumping.h"
138+
// #include "src/linear_algebra/eigen/sparse_block_matrix.h"
139+
// #include "src/linear_algebra/eigen/woodbury.h"
140+
// // randomized linear algebra
141+
// #include "src/linear_algebra/eigen/rsi.h"
142+
// #include "src/linear_algebra/eigen/rbki.h"
143+
// #include "src/linear_algebra/eigen/rp_chol.h"
139144

140145
// clang-format on
141146

‎fdaPDE/src/linear_algebra/binary_op.h‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -14,10 +14,10 @@
1414
// You should have received a copy of the GNU General Public License
1515
// along with this program. If not, see <http://www.gnu.org/licenses/>.
1616

17-
#ifndef __FDAPDE_MATRIX_BINARY_OP_H__
18-
#define __FDAPDE_MATRIX_BINARY_OP_H__
17+
#ifndef __FDAPDE_LINALG_BINARY_OP_H__
18+
#define __FDAPDE_LINALG_BINARY_OP_H__
1919

20-
#include "../header_check.h"
20+
#include "header_check.h"
2121

2222
namespace fdapde {
2323

@@ -102,7 +102,7 @@ namespace internals {
102102

103103
template <typename Scalar> struct matrix_coeff_mult_t {
104104
constexpr explicit matrix_coeff_mult_t(Scalar x) noexcept : x_(x) { }
105-
template <typename Scalar_> constexpr auto operator()(const Scalar_& y) { return x_ * y; }
105+
template <typename Scalar_> constexpr auto operator()(const Scalar_& y) const { return x_ * y; }
106106
private:
107107
Scalar x_;
108108
};
@@ -278,4 +278,4 @@ constexpr MatrixKroneckerProductOp<LhsXprType, RhsXprType> kron(
278278

279279
} // namespace fdapde
280280

281-
#endif // __FDAPDE_MATRIX_BINARY_OP_H__
281+
#endif // __FDAPDE_LINALG_BINARY_OP_H__

‎fdaPDE/src/linear_algebra/block.h‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -14,10 +14,10 @@
1414
// You should have received a copy of the GNU General Public License
1515
// along with this program. If not, see <http://www.gnu.org/licenses/>.
1616

17-
#ifndef __FDAPDE_MATRIX_BLOCK_H__
18-
#define __FDAPDE_MATRIX_BLOCK_H__
17+
#ifndef __FDAPDE_LINALG_BLOCK_H__
18+
#define __FDAPDE_LINALG_BLOCK_H__
1919

20-
#include "../header_check.h"
20+
#include "header_check.h"
2121

2222
namespace fdapde {
2323

@@ -127,10 +127,10 @@ class MatrixBlock : public MatrixExpr<BlockRows_, BlockCols_, MatrixBlock<BlockR
127127
}
128128
private:
129129
int start_row_ = 0, start_col_ = 0;
130-
std::conditional_t<Rows == Dynamic, int, empty_t> block_rows_, block_cols_;
130+
int block_rows_ = 0, block_cols_ = 0;
131131
XprTypeNested xpr_;
132132
};
133133

134134
} // namespace fdapde
135135

136-
#endif // __FDAPDE_MATRIX_BLOCK_H__
136+
#endif // __FDAPDE_LINALG_BLOCK_H__

‎fdaPDE/src/linear_algebra/diagonal.h‎

Lines changed: 41 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -14,10 +14,10 @@
1414
// You should have received a copy of the GNU General Public License
1515
// along with this program. If not, see <http://www.gnu.org/licenses/>.
1616

17-
#ifndef __FDAPDE_DIAGONAL_MATRIX_H__
18-
#define __FDAPDE_DIAGONAL_MATRIX_H__
17+
#ifndef __FDAPDE_LINALG_DIAGONAL_H__
18+
#define __FDAPDE_LINALG_DIAGONAL_H__
1919

20-
#include "../header_check.h"
20+
#include "header_check.h"
2121

2222
namespace fdapde {
2323

@@ -86,8 +86,43 @@ struct Diagonal : public MatrixExpr<Rows_, Cols_, Diagonal<Rows_, Cols_, XprType
8686
XprTypeNested xpr_;
8787
};
8888

89-
template <typename Scalar, int DiagonalSize> class DiagonalMatrix;
89+
template <int Rows_, int Cols_, typename XprType> struct DiagonalMatrixExpr;
9090

91+
namespace internals {
92+
93+
// class wrapping an expression of the diagonal coefficients to a square matrix. internal usage only
94+
template <int Rows_, int Cols_, typename DiagonalXprType>
95+
struct diagonal_wrapper : DiagonalMatrixExpr<Rows_, Cols_, diagonal_wrapper<Rows_, Cols_, DiagonalXprType>> {
96+
using Base = DiagonalMatrixExpr<Rows_, Cols_, diagonal_wrapper<Rows_, Cols_, DiagonalXprType>>;
97+
using DiagonalXprTypeNested = internals::ref_select_t<const DiagonalXprType>;
98+
using Scalar = typename DiagonalXprType::Scalar;
99+
static constexpr int Rows = Rows_;
100+
static constexpr int Cols = Cols_;
101+
static constexpr int NestAsRef = 0;
102+
static constexpr int ReadOnly = 1;
103+
104+
template <typename XprType>
105+
requires(std::is_constructible_v<DiagonalXprTypeNested, XprType>)
106+
constexpr diagonal_wrapper(XprType&& xpr) : Base(), xpr_(std::forward<XprType>(xpr)) { }
107+
constexpr Scalar operator()(int i, int j) const {
108+
fdapde_constexpr_assert(i >= 0 && i < Base::size_ && j >= 0 && j < Base::size_);
109+
return i == j ? xpr_[i] : Scalar(0);
110+
}
111+
private:
112+
DiagonalXprTypeNested xpr_;
113+
};
114+
115+
// helper cast function
116+
template <typename XprType> auto diagonal_cast(XprType&& xpr) {
117+
using XprTypeClean = std::decay_t<XprType>;
118+
static constexpr int Rows = XprTypeClean::Rows;
119+
static constexpr int Cols = XprTypeClean::Cols;
120+
fdapde_static_assert(Rows == 1 || Cols == 1, THIS_METHOD_IS_FOR_ROW_OR_COLUMN_VECTORS_ONLY);
121+
return diagonal_wrapper<Rows == 1 ? Cols : Rows, Cols == 1 ? Rows : Cols, XprType>(xpr);
122+
}
123+
124+
} // namespace internals
125+
91126
template <int Rows_, int Cols_, typename XprType>
92127
struct DiagonalMatrixExpr : public MatrixExpr<Rows_, Cols_, XprType> {
93128
using Base = MatrixExpr<Rows_, Cols_, XprType>;
@@ -126,7 +161,7 @@ struct DiagonalMatrixExpr : public MatrixExpr<Rows_, Cols_, XprType> {
126161
// converts the diagonal expression to a full dense matrix
127162
auto to_matrix() const { return Matrix<typename XprType::Scalar, Rows, Cols>(derived()); }
128163
// matrix inverse as 1/coeff
129-
auto inverse() const { return DiagonalMatrix<typename XprType::Scalar, Rows>(derived().cwise_inv()); }
164+
auto inverse() const { return internals::diagonal_cast(derived().diagonal().cwise_inv()); }
130165
// linear system solver Ax = b
131166
template <typename RhsXprType> constexpr auto solve(const RhsXprType& b) const {
132167
using Scalar = typename XprType::Scalar;
@@ -143,41 +178,6 @@ struct DiagonalMatrixExpr : public MatrixExpr<Rows_, Cols_, XprType> {
143178
int size_;
144179
};
145180

146-
namespace internals {
147-
148-
// class wrapping an expression of the diagonal coefficients to a square matrix. internal usage only
149-
template <int Rows_, int Cols_, typename DiagonalXprType>
150-
struct diagonal_wrapper : DiagonalMatrixExpr<Rows_, Cols_, diagonal_wrapper<Rows_, Cols_, DiagonalXprType>> {
151-
using Base = DiagonalMatrixExpr<Rows_, Cols_, diagonal_wrapper<Rows_, Cols_, DiagonalXprType>>;
152-
using DiagonalXprTypeNested = internals::ref_select_t<const DiagonalXprType>;
153-
using Scalar = typename DiagonalXprType::Scalar;
154-
static constexpr int Rows = Rows_;
155-
static constexpr int Cols = Cols_;
156-
static constexpr int NestAsRef = 0;
157-
static constexpr int ReadOnly = 1;
158-
159-
template <typename XprType>
160-
requires(std::is_constructible_v<DiagonalXprTypeNested, XprType>)
161-
constexpr diagonal_wrapper(XprType&& xpr) : Base(), xpr_(std::forward<XprType>(xpr)) { }
162-
constexpr Scalar operator()(int i, int j) const {
163-
fdapde_constexpr_assert(i >= 0 && i < Base::size_ && j >= 0 && j < Base::size_);
164-
return i == j ? xpr_[i] : Scalar(0);
165-
}
166-
private:
167-
DiagonalXprTypeNested xpr_;
168-
};
169-
170-
// helper cast function
171-
template <int Rows, int Cols, typename XprType> auto diagonal_cast(XprType&& xpr) {
172-
using XprTypeClean = std::decay_t<XprType>;
173-
static constexpr int Rows = XprTypeClean::Rows;
174-
static constexpr int Cols = XprTypeClean::Cols;
175-
fdapde_static_assert(Rows == 1 || Cols == 1, THIS_METHOD_IS_FOR_ROW_OR_COLUMN_VECTORS_ONLY);
176-
return diagonal_wrapper<Rows == 1 ? Cols : Rows, Cols == 1 ? Rows : Cols, XprType>(xpr);
177-
}
178-
179-
} // namespace internals
180-
181181
// diagonal matrix arithmetic (O(n) operations on the diagonal coefficients)
182182
template <typename LhsXprType, typename RhsXprType>
183183
constexpr auto operator+(
@@ -380,4 +380,4 @@ class DiagonalMatrixView : public DiagonalMatrixExpr<Rows_, Rows_, DiagonalMatri
380380

381381
} // namespace fdapde
382382

383-
#endif // __FDAPDE_DIAGONAL_MATRIX_H__
383+
#endif // __FDAPDE_LINALG_DIAGONAL_H__

‎fdaPDE/src/linear_algebra/evd.h‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -14,8 +14,10 @@
1414
// You should have received a copy of the GNU General Public License
1515
// along with this program. If not, see <http://www.gnu.org/licenses/>.
1616

17-
#ifndef __FDAPDE_EVD_H__
18-
#define __FDAPDE_EVD_H__
17+
#ifndef __FDAPDE_LINALG_EVD_H__
18+
#define __FDAPDE_LINALG_EVD_H__
19+
20+
#include "header_check.h"
1921

2022
namespace fdapde {
2123

@@ -181,4 +183,4 @@ template <typename Scalar_, int Size_> class EVD {
181183

182184
} // namespace fdapde
183185

184-
#endif // __FDAPDE_EVD_H__
186+
#endif // __FDAPDE_LINALG_EVD_H__

0 commit comments

Comments
 (0)