Skip to content

Commit a1ea3e2

Browse files
AlePaluAlePalu
committed
matrix type subsystems detection traits, unified Rows, Cols interface across square and non-square matrix types (enforce squareness by assert, not by API), symmetric type system, stop delegating copy constructor/assignment to base, use this-> instead of Base:: whenever possible, fixed header guards, minor changes
1 parent 7956032 commit a1ea3e2

File tree

16 files changed

+532
-339
lines changed

16 files changed

+532
-339
lines changed

‎fdaPDE/linear_algebra.h‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -64,8 +64,7 @@ struct unchecked_t { };
6464
#include "src/linear_algebra/triangular.h"
6565
#include "src/linear_algebra/bool.h"
6666
#include "src/linear_algebra/coeffwise.h"
67-
68-
// #include "src/linear_algebra/symmetric.h"
67+
#include "src/linear_algebra/symmetric.h"
6968

7069
// // #include "src/linear_algebra/skew.h"
7170
#include "src/linear_algebra/permutation.h"
@@ -78,7 +77,7 @@ struct unchecked_t { };
7877
#include "src/linear_algebra/xpr.h"
7978

8079
// // algorithms
81-
// #include "src/linear_algebra/evd.h"
80+
#include "src/linear_algebra/evd.h"
8281

8382
// #include "src/linear_algebra/qr.h"
8483

‎fdaPDE/src/linear_algebra/bool.h‎

Lines changed: 17 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -962,6 +962,17 @@ constexpr bool operator!=(const BoolMatrixExpr<LhsXprType>& op1, const BoolMatri
962962
return !(op1 == op2);
963963
}
964964

965+
// detection trait
966+
template <typename XprType> struct is_boolean_matrix {
967+
static constexpr bool value = std::is_base_of_v<BoolMatrixExpr<std::decay_t<XprType>>, XprType>;
968+
};
969+
template <typename XprType> static constexpr bool is_boolean_matrix_v = is_boolean_matrix<XprType>::value;
970+
template <typename XprType> struct is_boolean_vector {
971+
static constexpr bool value =
972+
is_boolean_matrix_v<XprType> && (std::decay_t<XprType>::Cols == 1 || std::decay_t<XprType>::Rows == 1);
973+
};
974+
template <typename XprType> static constexpr bool is_boolean_vector_v = is_boolean_vector<XprType>::value;
975+
965976
// indexes of true elements in the boolean expression
966977
template <typename XprType> std::vector<int> which(const BoolMatrixExpr<XprType>& mtx) { return mtx.which(true); }
967978

@@ -991,7 +1002,8 @@ Matrix<bool, Dynamic, Dynamic> nan_indicator(DataType&& data) {
9911002
Matrix<bool, Dynamic, Dynamic> mask(rows, cols);
9921003
for (int i = 0; i < rows; ++i) {
9931004
for (int j = 0; j < cols; ++j) {
994-
const auto val = internals::is_vector_like_v<DataType> ? internals::vector_like_access(data, i) : data(i, j);
1005+
const auto val =
1006+
internals::is_vector_like_v<DataType> ? internals::vector_like_access(data, i) : data(i, j);
9951007
if (std::isnan(val)) mask.set(i, j);
9961008
}
9971009
}
@@ -1062,19 +1074,19 @@ class MatrixView<bool, Rows_, Cols_, StorageOrder_> :
10621074
constexpr const bitpack_t* data() const { return data_; }
10631075
constexpr bitpack_t* data() { return data_; }
10641076
// modifiers
1065-
constexpr void set(int i, int j) { Base::operator()(i, j).set(); }
1077+
constexpr void set(int i, int j) { this->operator()(i, j).set(); }
10661078
constexpr void set(int i) {
10671079
fdapde_static_assert(Rows == 1 || Cols == 1, THIS_METHOD_IS_FOR_ROW_OR_COLUMN_VECTORS_ONLY);
1068-
Base::operator[](i).set();
1080+
this->operator[](i).set();
10691081
}
10701082
constexpr void set() {
10711083
for (int i = 0; i < bitpacks_ - 1; ++i) { data_[i] = ~bitpack_t(0); }
10721084
data_[bitpacks_ - 1] |= last_bitpack_mask_;
10731085
}
1074-
constexpr void clear(int i, int j) { Base::operator()(i, j).clear(); }
1086+
constexpr void clear(int i, int j) { this->operator()(i, j).clear(); }
10751087
constexpr void clear(int i) {
10761088
fdapde_static_assert(Rows == 1 || Cols == 1, THIS_METHOD_IS_FOR_ROW_OR_COLUMN_VECTORS_ONLY);
1077-
Base::operator[](i).clear();
1089+
this->operator[](i).clear();
10781090
}
10791091
constexpr void clear() {
10801092
for (int i = 0; i < bitpacks_ - 1; ++i) { data_[i] = bitpack_t(0); }

‎fdaPDE/src/linear_algebra/coeffwise.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121

2222
namespace fdapde {
2323

24+
// coefficient-wise type system
2425
template <typename XprType_> struct MatrixCoeffWiseExpr;
2526

2627
namespace internals {

‎fdaPDE/src/linear_algebra/diagonal.h‎

Lines changed: 22 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -22,8 +22,8 @@
2222
namespace fdapde {
2323

2424
// diagonal matrix type system
25-
template <typename XprType> struct Diagonal;
26-
25+
template <typename XprType> struct DiagonalMatrixExpr;
26+
2727
namespace internals {
2828

2929
struct diagonal_assignment_executor {
@@ -48,7 +48,7 @@ struct diagonal_assignment_executor {
4848

4949
} // namespace internals
5050

51-
// expression of the diagonal of a matrix
51+
// expression of the diagonal of a matrix (this node acts as a vector expression)
5252
template <typename XprType_> struct Diagonal : public MatrixExpr<Diagonal<XprType_>> {
5353
private:
5454
using Base = MatrixExpr<Diagonal<XprType_>>;
@@ -84,18 +84,14 @@ template <typename XprType_> struct Diagonal : public MatrixExpr<Diagonal<XprTyp
8484
using Base::operator=;
8585
// const access
8686
constexpr Scalar operator()(int i, int j) const {
87-
fdapde_assert(i >= 0 && i < rows() && j >= 0 && j < cols());
87+
fdapde_assert(i >= 0 && i < rows() && j >= 0 && j < 1);
8888
return xpr_(i, i);
8989
}
9090
constexpr Scalar operator[](int i) const {
9191
fdapde_assert(i >= 0 && i < rows());
9292
return xpr_(i, i);
9393
}
9494
// non-const access
95-
constexpr Scalar& operator()(int i, int j) {
96-
fdapde_assert(i >= 0 && i < rows() && j >= 0 && j < cols());
97-
return xpr_(i, i);
98-
}
9995
constexpr Scalar& operator[](int i) {
10096
fdapde_assert(i >= 0 && i < rows());
10197
return xpr_(i, i);
@@ -108,8 +104,6 @@ template <typename XprType_> struct Diagonal : public MatrixExpr<Diagonal<XprTyp
108104
XprTypeNested xpr_;
109105
};
110106

111-
template <typename XprType> struct DiagonalMatrixExpr;
112-
113107
namespace internals {
114108

115109
// class wrapping an expression of the diagonal coefficients to a square matrix. internal usage only
@@ -294,15 +288,20 @@ class DiagonalMatrix : public DiagonalMatrixExpr<DiagonalMatrix<Scalar_, Rows_>>
294288
using assignment_executor = internals::diagonal_assignment_executor;
295289

296290
constexpr DiagonalMatrix() : Base(), data_() { }
291+
// copy semantic
292+
constexpr DiagonalMatrix(const DiagonalMatrix& rhs) : Base() { clone_(rhs); }
293+
constexpr DiagonalMatrix& operator=(const DiagonalMatrix& rhs) {
294+
clone_(rhs);
295+
return *this;
296+
}
297297
constexpr explicit DiagonalMatrix(int size) : Base(size), data_() {
298298
fdapde_static_assert(
299299
Rows == Dynamic || Cols == Dynamic, THIS_METHOD_IS_FOR_DYNAMIC_SIZED_DIAGONAL_MATRICES_ONLY);
300300
data_.resize(size);
301301
}
302302
template <typename RhsXprType_> constexpr DiagonalMatrix(const MatrixExpr<RhsXprType_>& rhs) : Base() {
303303
fdapde_assert(this->rows() == rhs.rows() && this->cols() == rhs.cols());
304-
if constexpr (Rows == Dynamic || Cols == Dynamic) { resize(rhs.rows()); }
305-
assignment_executor::run(*this, rhs.derived(), [](auto&& l, const auto& r) { l = r; });
304+
clone_(rhs.derived());
306305
}
307306
constexpr explicit DiagonalMatrix(const std::vector<Scalar>& vec) : Base(vec.size()), data_() {
308307
if constexpr (Rows == Dynamic || Cols == Dynamic) { data_.resize(vec.size()); }
@@ -331,8 +330,6 @@ class DiagonalMatrix : public DiagonalMatrixExpr<DiagonalMatrix<Scalar_, Rows_>>
331330
data_[2] = z;
332331
}
333332
constexpr const StorageType& diagonal() const { return data_; }
334-
// inherit assignment from Base
335-
using Base::operator=;
336333
// modifiers
337334
void resize(int size) {
338335
fdapde_static_assert(Rows == Dynamic || Cols == Dynamic, THIS_METHOD_IS_FOR_DYNAMIC_SIZED_MATRICES_ONLY);
@@ -348,6 +345,11 @@ class DiagonalMatrix : public DiagonalMatrixExpr<DiagonalMatrix<Scalar_, Rows_>>
348345
constexpr const Scalar* data() const { return data_.data(); }
349346
constexpr Scalar* data() { return data_.data(); }
350347
private:
348+
template <typename RhsXprType> constexpr void clone_(const RhsXprType& rhs) {
349+
if constexpr (Rows == Dynamic || Cols == Dynamic) { resize(rhs.rows()); }
350+
assignment_executor::run(*this, rhs, [](auto&& l, const auto& r) { l = r; });
351+
return;
352+
}
351353
StorageType data_;
352354
};
353355

@@ -384,6 +386,12 @@ class DiagonalMatrixView : public DiagonalMatrixExpr<DiagonalMatrixView<Scalar_,
384386
StorageType data_;
385387
};
386388

389+
// detection trait
390+
template <typename XprType> struct is_diagonal_matrix {
391+
static constexpr bool value = std::is_base_of_v<DiagonalMatrixExpr<std::decay_t<XprType>>, XprType>;
392+
};
393+
template <typename XprType> static constexpr bool is_diagonal_matrix_v = is_diagonal_matrix<XprType>::value;
394+
387395
} // namespace fdapde
388396

389397
#endif // __FDAPDE_LINALG_DIAGONAL_H__

‎fdaPDE/src/linear_algebra/evd.h‎

Lines changed: 16 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -27,17 +27,18 @@ template <typename XprType_> class EVD {
2727
fdapde_static_assert(
2828
XprType::Rows == Dynamic || XprType::Cols == Dynamic || XprType::Rows == XprType::Cols,
2929
THIS_CLASS_IS_FOR_SQUARE_MATRICES_ONLY);
30-
static constexpr int Size = XprType::Rows;
30+
static constexpr int Rows = XprType::Rows;
31+
static constexpr int Cols = XprType::Cols;
3132
using Scalar = typename XprType::Scalar;
3233

3334
// tridiagonalize a symmetric matrix via householder reflectors
3435
// see "Golub, G. H., & Van Loan, C. F. (2013). Matrix computations. JHU press. Sec.8.3.1"
3536
template <typename XprType>
36-
std::tuple<Matrix<Scalar, Size, Size>, OrthogonalMatrix<Scalar, Size>>
37+
std::tuple<Matrix<Scalar, Rows, Cols>, OrthogonalMatrix<Scalar, Rows, Cols>>
3738
householder_tridiagonalize_(const XprType& m) const {
3839
const int n = m.rows();
39-
Matrix<Scalar, Size, Size> T = m;
40-
Matrix<Scalar, Size, Size> Q = Matrix<Scalar, Size, Size>::Identity(n, n);
40+
Matrix<Scalar, Rows, Cols> T = m;
41+
Matrix<Scalar, Rows, Cols> Q = Matrix<Scalar, Rows, Cols>::Identity(n, n);
4142

4243
for (int k = 0; k < n - 2; ++k) {
4344
const int m = n - k - 1;
@@ -70,32 +71,29 @@ template <typename XprType_> class EVD {
7071
Vector<Scalar, Dynamic> zQ = Q.block(0, k + 1, n, m) * u;
7172
Q.block(0, k + 1, n, m) -= zQ * (beta * u.transpose());
7273
}
73-
return std::make_pair(T, internals::orthogonal_wrapper<Size, Size, Matrix<Scalar, Size, Size>>(Q));
74+
return std::make_pair(T, internals::orthogonal_wrapper<Matrix<Scalar, Rows, Cols>>(Q));
7475
}
7576

7677
int max_iter_ = 30; // taken from LAPACK, actual number of iteration is scaled by matrix size
7778
public:
7879
constexpr EVD() = default;
79-
80-
template <int Rows, int Cols, typename XprType>
81-
constexpr explicit EVD(const SymmetricMatrixExpr<Rows, Cols, XprType>& m) {
80+
template <typename XprType__>
81+
constexpr explicit EVD(const SymmetricMatrixExpr<XprType__>& m) {
8282
if constexpr (Rows == Dynamic || Cols == Dynamic) { fdapde_assert(m.rows() == m.cols()); }
8383
compute(m);
8484
}
8585

8686
// computes the EVD of a symmetric matrix using the implicit QR-iteration with Wilkinson shift
8787
// see "Golub, G. H., & Van Loan, C. F. (2013). Matrix computations. JHU press. Ch.8.3"
88-
template <int Rows, int Cols, typename XprType>
89-
constexpr void compute(const SymmetricMatrixExpr<Rows, Cols, XprType>& mtx) {
90-
fdapde_static_assert(Rows == Cols && Rows == Size, THIS_METHOD_IS_FOR_SQUARE_MATRICES_ONLY);
88+
template <typename XprType__> constexpr void compute(const SymmetricMatrixExpr<XprType__>& mtx) {
9189
auto [T, Q_] = householder_tridiagonalize_(mtx.derived());
9290
const int n = T.rows();
9391
const int max_iter = max_iter_ * n;
9492
Matrix<Scalar, Rows, Cols> Q = Matrix<Scalar, Rows, Cols>::Identity(n, n);
9593
// extract diagonal and subdiagonal
96-
Vector<Scalar, Size> dd;
97-
Vector<Scalar, Size == Dynamic ? Dynamic : (Size - 1)> sd;
98-
if constexpr (Size == Dynamic) {
94+
Vector<Scalar, Rows> dd;
95+
Vector<Scalar, Rows == Dynamic ? Dynamic : (Rows - 1)> sd;
96+
if constexpr (Rows == Dynamic) {
9997
dd.resize(n);
10098
sd.resize(n - 1);
10199
}
@@ -177,13 +175,13 @@ template <typename XprType_> class EVD {
177175
eigenvalues_ = dd;
178176
}
179177
// observers
180-
constexpr const Vector<Scalar, Size>& eigenvalues() const { return eigenvalues_; }
178+
constexpr const Vector<Scalar, Rows>& eigenvalues() const { return eigenvalues_; }
181179
constexpr auto eigenvectors() const {
182-
return internals::orthogonal_wrapper<Size, Size, Matrix<Scalar, Size, Size>>(eigenvectors_);
180+
return internals::orthogonal_wrapper<Matrix<Scalar, Rows, Cols>>(eigenvectors_);
183181
}
184182
private:
185-
Matrix<Scalar, Size, Size> eigenvectors_;
186-
Vector<Scalar, Size> eigenvalues_;
183+
Matrix<Scalar, Rows, Cols> eigenvectors_;
184+
Vector<Scalar, Rows> eigenvalues_;
187185
};
188186

189187
} // namespace fdapde

‎fdaPDE/src/linear_algebra/matrix.h‎

Lines changed: 15 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -216,25 +216,18 @@ class Matrix : public MatrixBase<Scalar_, Rows_, Cols_, StorageOrder_, Matrix<Sc
216216
static constexpr int NestAsRef = 1;
217217

218218
constexpr Matrix() : data_() { }
219-
constexpr Matrix(const Matrix& other) : Base() {
220-
if constexpr (Rows_ == Dynamic || Cols_ == Dynamic) { resize(other.rows(), other.cols()); }
221-
using assignment = typename Base::assignment_executor;
222-
assignment::run(*this, other, [](Scalar& l, const Scalar& r) { l = r; });
223-
}
219+
// copy semantic
220+
constexpr Matrix(const Matrix& other) : Base() { clone_(other); }
224221
constexpr Matrix& operator=(const Matrix& other) {
225-
Base::operator=(other);
222+
clone_(other);
226223
return *this;
227224
}
228225
template <typename RhsXprType_> // construct from plain MatrixExpr
229226
constexpr Matrix(const MatrixExpr<RhsXprType_>& rhs) : Base(), data_() {
230-
if constexpr (Rows_ == Dynamic || Cols_ == Dynamic) { resize(rhs.rows(), rhs.cols()); }
231-
using assignment = typename Base::assignment_executor;
232-
assignment::run(*this, rhs.derived(), [](Scalar& l, const Scalar& r) { l = r; });
227+
clone_(rhs.derived());
233228
}
234229
template <typename RhsXprType_>
235230
constexpr Matrix(const MatrixCoeffWiseExpr<RhsXprType_>& rhs) : Matrix(rhs.mwise()) { }
236-
// inherit assignment from base
237-
using Base::operator=;
238231

239232
// Matrix API
240233
// value-initialized static-sized matrix, avoid 1D vectors
@@ -370,12 +363,12 @@ class Matrix : public MatrixBase<Scalar_, Rows_, Cols_, StorageOrder_, Matrix<Sc
370363
fdapde_static_assert(Rows_ == Dynamic || Cols_ == Dynamic, THIS_METHOD_IS_FOR_DYNAMIC_SIZED_MATRICES_ONLY);
371364
const int rows_ = Rows_ == Dynamic ? rows : Rows_;
372365
const int cols_ = Cols_ == Dynamic ? cols : Cols_;
373-
if (rows_ == Base::rows_ && cols_ == Base::cols_) return; // do not reallocate memory if sizes didn't changed
366+
if (rows_ == this->rows_ && cols_ == this->cols_) return; // do not reallocate memory if sizes didn't changed
374367
// update and reallocate memory
375-
Base::rows_ = rows_;
376-
Base::cols_ = cols_;
377-
Base::row_stride_ = StorageOrder_ == RowMajor ? cols_ : 1;
378-
Base::col_stride_ = StorageOrder_ == RowMajor ? 1 : rows_;
368+
this->rows_ = rows_;
369+
this->cols_ = cols_;
370+
this->row_stride_ = StorageOrder_ == RowMajor ? cols_ : 1;
371+
this->col_stride_ = StorageOrder_ == RowMajor ? 1 : rows_;
379372
data_.resize(rows * cols);
380373
return;
381374
}
@@ -395,6 +388,12 @@ class Matrix : public MatrixBase<Scalar_, Rows_, Cols_, StorageOrder_, Matrix<Sc
395388
constexpr iterator end() { return data_.end(); }
396389
constexpr const_iterator end() const { return data_.end(); }
397390
private:
391+
template <typename RhsXprType> constexpr void clone_(const RhsXprType& rhs) {
392+
if constexpr (Rows_ == Dynamic || Cols_ == Dynamic) { resize(rhs.rows(), rhs.cols()); }
393+
using assignment_executor = typename Base::assignment_executor;
394+
assignment_executor::run(*this, rhs, [](auto&& l, const auto& r) { l = r; });
395+
return;
396+
}
398397
StorageType data_;
399398
};
400399

0 commit comments

Comments
 (0)