diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index 1f7fe37e571447b104b81ff61cf38dcf2144d897..c2df50e454fa883bf402c56ac81f8617cf9b7742 100644
--- a/src/algebra/CRSMatrix.hpp
+++ b/src/algebra/CRSMatrix.hpp
@@ -27,6 +27,117 @@ class CRSMatrix
Array<const DataType> m_values;
Array<const IndexType> m_column_indices;
+ template <typename DataType2, typename BinOp>
+ CRSMatrix
+ _binOp(const CRSMatrix& A, const CRSMatrix<DataType2, IndexType>& B) const
+ {
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot add matrices of different types");
+ Assert(A.numberOfRows() == B.numberOfRows() and A.numberOfColumns() == B.numberOfColumns(),
+ "cannot apply inner binary operator on matrices of different sizes");
+
+ Array<int> non_zeros(A.m_nb_rows);
+ for (IndexType i_row = 0; i_row < A.m_nb_rows; ++i_row) {
+ const auto A_row_begin = A.m_row_map[i_row];
+ const auto A_row_end = A.m_row_map[i_row + 1];
+ const auto B_row_begin = B.m_row_map[i_row];
+ const auto B_row_end = B.m_row_map[i_row + 1];
+ IndexType i_row_A = A_row_begin;
+ IndexType i_row_B = B_row_begin;
+
+ int row_nb_columns = 0;
+
+ while (i_row_A < A_row_end or i_row_B < B_row_end) {
+ const IndexType A_column_idx = [&] {
+ if (i_row_A < A_row_end) {
+ return A.m_column_indices[i_row_A];
+ } else {
+ return std::numeric_limits<IndexType>::max();
+ }
+ }();
+
+ const IndexType B_column_idx = [&] {
+ if (i_row_B < B_row_end) {
+ return B.m_column_indices[i_row_B];
+ } else {
+ return std::numeric_limits<IndexType>::max();
+ }
+ }();
+
+ if (A_column_idx == B_column_idx) {
+ ++row_nb_columns;
+ ++i_row_A;
+ ++i_row_B;
+ } else if (A_column_idx < B_column_idx) {
+ ++row_nb_columns;
+ ++i_row_A;
+ } else {
+ Assert(B_column_idx < A_column_idx);
+ ++row_nb_columns;
+ ++i_row_B;
+ }
+ }
+ non_zeros[i_row] = row_nb_columns;
+ }
+
+ Array<IndexType> row_map(A.m_nb_rows + 1);
+ row_map[0] = 0;
+ for (IndexType i = 0; i < A.m_nb_rows; ++i) {
+ row_map[i + 1] = row_map[i] + non_zeros[i];
+ }
+
+ const IndexType nb_values = row_map[row_map.size() - 1];
+ Array<DataType> values(nb_values);
+ Array<IndexType> column_indices(nb_values);
+
+ IndexType i_value = 0;
+ for (IndexType i_row = 0; i_row < A.m_nb_rows; ++i_row) {
+ const auto A_row_begin = A.m_row_map[i_row];
+ const auto A_row_end = A.m_row_map[i_row + 1];
+ const auto B_row_begin = B.m_row_map[i_row];
+ const auto B_row_end = B.m_row_map[i_row + 1];
+ IndexType i_row_A = A_row_begin;
+ IndexType i_row_B = B_row_begin;
+
+ while (i_row_A < A_row_end or i_row_B < B_row_end) {
+ const IndexType A_column_idx = [&] {
+ if (i_row_A < A_row_end) {
+ return A.m_column_indices[i_row_A];
+ } else {
+ return std::numeric_limits<IndexType>::max();
+ }
+ }();
+
+ const IndexType B_column_idx = [&] {
+ if (i_row_B < B_row_end) {
+ return B.m_column_indices[i_row_B];
+ } else {
+ return std::numeric_limits<IndexType>::max();
+ }
+ }();
+
+ if (A_column_idx == B_column_idx) {
+ column_indices[i_value] = A_column_idx;
+ values[i_value] = BinOp()(A.m_values[i_row_A], B.m_values[i_row_B]);
+ ++i_row_A;
+ ++i_row_B;
+ } else if (A_column_idx < B_column_idx) {
+ column_indices[i_value] = A_column_idx;
+ values[i_value] = BinOp()(A.m_values[i_row_A], 0);
+ ++i_row_A;
+ } else {
+ Assert(B_column_idx < A_column_idx);
+ column_indices[i_value] = B_column_idx;
+ values[i_value] = BinOp()(0, B.m_values[i_row_B]);
+ ++i_row_B;
+ }
+ ++i_value;
+ }
+ }
+
+ return CRSMatrix(A.m_nb_rows, A.m_nb_columns, row_map, values, column_indices);
+ }
+
public:
PUGS_INLINE
bool
@@ -72,9 +183,9 @@ class CRSMatrix
operator*(const Vector<DataType2>& x) const
{
static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
- "Cannot multiply matrix and vector of different type");
+ "Cannot multiply matrix and vector of different types");
- Assert(x.size() - m_nb_columns == 0);
+ Assert(x.size() - m_nb_columns == 0, "cannot compute product: incompatible matrix and vector sizes");
Vector<MutableDataType> Ax(m_nb_rows);
@@ -92,6 +203,20 @@ class CRSMatrix
return Ax;
}
+ template <typename DataType2>
+ CRSMatrix<DataType>
+ operator+(const CRSMatrix<DataType2>& B) const
+ {
+ return this->_binOp<DataType2, std::plus<DataType>>(*this, B);
+ }
+
+ template <typename DataType2>
+ CRSMatrix<DataType>
+ operator-(const CRSMatrix<DataType2>& B) const
+ {
+ return this->_binOp<DataType2, std::minus<DataType>>(*this, B);
+ }
+
friend PUGS_INLINE std::ostream&
operator<<(std::ostream& os, const CRSMatrix& A)
{
diff --git a/src/algebra/CRSMatrixDescriptor.hpp b/src/algebra/CRSMatrixDescriptor.hpp
index 4be6a323630a9845ecb142b32f73ceca5d6dc7f6..9c2e0894152ead3e01d4e38396f581a05c283fa2 100644
--- a/src/algebra/CRSMatrixDescriptor.hpp
+++ b/src/algebra/CRSMatrixDescriptor.hpp
@@ -248,9 +248,9 @@ class CRSMatrixDescriptor
}
}
- CRSMatrixDescriptor(const IndexType nb_rows, const IndexType nb_columns, const Array<IndexType>& non_zeros)
- : m_nb_rows{nb_rows},
- m_nb_columns{nb_columns},
+ CRSMatrixDescriptor(const size_t nb_rows, const size_t nb_columns, const Array<IndexType>& non_zeros)
+ : m_nb_rows(nb_rows),
+ m_nb_columns(nb_columns),
m_row_map{this->_computeRowMap(non_zeros)},
m_values(m_row_map[m_row_map.size() - 1]),
m_column_indices(m_row_map[m_row_map.size() - 1]),
@@ -268,6 +268,10 @@ class CRSMatrixDescriptor
}
}
+ CRSMatrixDescriptor(const size_t nb_rows, const Array<IndexType>& non_zeros)
+ : CRSMatrixDescriptor(nb_rows, nb_rows, non_zeros)
+ {}
+
~CRSMatrixDescriptor() = default;
};
diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 1d96d01f992d22a8346417dce85750d9b186749b..6125e67e6c5fad94a122028ce78bd78ef7b12e6c 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -63,7 +63,7 @@ struct EigenvalueSolver::Internals
};
void
-EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues)
+EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues)
{
EPS eps;
@@ -76,7 +76,7 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Arr
PetscInt nb_eigenvalues;
EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
- eigenvalues = Array<double>(nb_eigenvalues);
+ eigenvalues = SmallArray<double>(nb_eigenvalues);
for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, nullptr, nullptr);
}
@@ -86,8 +86,8 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Arr
void
EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
- Array<double>& eigenvalues,
- std::vector<Vector<double>>& eigenvector_list)
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvector_list)
{
EPS eps;
@@ -100,11 +100,11 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
PetscInt nb_eigenvalues;
EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
- eigenvalues = Array<double>(nb_eigenvalues);
+ eigenvalues = SmallArray<double>(nb_eigenvalues);
eigenvector_list.reserve(nb_eigenvalues);
for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
Vec Vr;
- Vector<double> eigenvector{A.numberOfRows()};
+ SmallVector<double> eigenvector{A.numberOfRows()};
VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr);
EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr);
VecDestroy(&Vr);
@@ -116,8 +116,8 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
void
EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
- Array<double>& eigenvalues,
- DenseMatrix<double>& P)
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P)
{
EPS eps;
@@ -130,8 +130,8 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
PetscInt nb_eigenvalues;
EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
- eigenvalues = Array<double>(nb_eigenvalues);
- P = DenseMatrix<double>(nb_eigenvalues, nb_eigenvalues);
+ eigenvalues = SmallArray<double>(nb_eigenvalues);
+ P = SmallMatrix<double>(nb_eigenvalues, nb_eigenvalues);
Array<double> eigenvector(nb_eigenvalues);
for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index a9dd941b9bd87d3f0a79be5f36e5683ae8a1afb6..a9b462509514ace4b736424c65f7c9abe4572f26 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -3,10 +3,10 @@
#include <algebra/PETScUtils.hpp>
-#include <algebra/DenseMatrix.hpp>
-#include <algebra/Vector.hpp>
-#include <utils/Array.hpp>
+#include <algebra/SmallMatrix.hpp>
+#include <algebra/SmallVector.hpp>
#include <utils/Exceptions.hpp>
+#include <utils/SmallArray.hpp>
class EigenvalueSolver
{
@@ -14,19 +14,21 @@ class EigenvalueSolver
struct Internals;
#ifdef PUGS_HAS_SLEPC
- void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues);
+ void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues);
void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
- Array<double>& eigenvalues,
- std::vector<Vector<double>>& eigenvectors);
+ SmallArray<double>& eigenvalues,
+ std::vector<SmallVector<double>>& eigenvectors);
- void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues, DenseMatrix<double>& P);
+ void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
+ SmallArray<double>& eigenvalues,
+ SmallMatrix<double>& P);
#endif // PUGS_HAS_SLEPC
public:
template <typename MatrixType>
void
- computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] Array<double>& eigenvalues)
+ computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] SmallArray<double>& eigenvalues)
{
#ifdef PUGS_HAS_SLEPC
this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues);
@@ -38,8 +40,8 @@ class EigenvalueSolver
template <typename MatrixType>
void
computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A,
- [[maybe_unused]] Array<double>& eigenvalues,
- [[maybe_unused]] std::vector<Vector<double>>& eigenvectors)
+ [[maybe_unused]] SmallArray<double>& eigenvalues,
+ [[maybe_unused]] std::vector<SmallVector<double>>& eigenvectors)
{
#ifdef PUGS_HAS_SLEPC
this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, eigenvectors);
@@ -51,8 +53,8 @@ class EigenvalueSolver
template <typename MatrixType>
void
computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A,
- [[maybe_unused]] Array<double>& eigenvalues,
- [[maybe_unused]] DenseMatrix<double>& P)
+ [[maybe_unused]] SmallArray<double>& eigenvalues,
+ [[maybe_unused]] SmallMatrix<double>& P)
{
#ifdef PUGS_HAS_SLEPC
this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, P);
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index 3cffca7298b0a7377996c130195ac9e43a66ca44..99bff927c744c9ba216e575dcc914a2bdba1d4da 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -66,7 +66,7 @@ struct LinearSolver::Internals
case LSMethod::bicgstab2:
case LSMethod::gmres:
case LSMethod::lu:
- case LSMethod::choleski: {
+ case LSMethod::cholesky: {
break;
}
// LCOV_EXCL_START
@@ -97,7 +97,7 @@ struct LinearSolver::Internals
case LSPrecond::none:
case LSPrecond::amg:
case LSPrecond::diagonal:
- case LSPrecond::incomplete_choleski:
+ case LSPrecond::incomplete_cholesky:
case LSPrecond::incomplete_LU: {
break;
}
@@ -221,9 +221,10 @@ struct LinearSolver::Internals
direct_solver = true;
break;
}
- case LSMethod::choleski: {
+ case LSMethod::cholesky: {
KSPSetType(ksp, KSPPREONLY);
PCSetType(pc, PCCHOLESKY);
+ PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
direct_solver = true;
break;
}
@@ -248,7 +249,7 @@ struct LinearSolver::Internals
PCSetType(pc, PCILU);
break;
}
- case LSPrecond::incomplete_choleski: {
+ case LSPrecond::incomplete_cholesky: {
PCSetType(pc, PCICC);
break;
}
@@ -335,7 +336,7 @@ LinearSolver::solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>&
}
void
-LinearSolver::solveLocalSystem(const DenseMatrix<double>& A, Vector<double>& x, const Vector<const double>& b)
+LinearSolver::solveLocalSystem(const SmallMatrix<double>& A, SmallVector<double>& x, const SmallVector<const double>& b)
{
Internals::solveLocalSystem(A, x, b, m_options);
}
diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
index 1ab43a1405ee312fdb9255c6b1f4508ade4e168e..e90402c8ddc78c9426f97300ab6b10730b4c503b 100644
--- a/src/algebra/LinearSolver.hpp
+++ b/src/algebra/LinearSolver.hpp
@@ -2,8 +2,8 @@
#define LINEAR_SOLVER_HPP
#include <algebra/CRSMatrix.hpp>
-#include <algebra/DenseMatrix.hpp>
#include <algebra/LinearSolverOptions.hpp>
+#include <algebra/SmallMatrix.hpp>
#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
#include <algebra/Vector.hpp>
@@ -24,13 +24,10 @@ class LinearSolver
void
solveLocalSystem(const TinyMatrix<N>& A, TinyVector<N>& x, const TinyVector<N>& b)
{
- DenseMatrix A_dense{A};
- Vector<double> x_vector{N};
- Vector<double> b_vector{N};
- for (size_t i = 0; i < N; ++i) {
- x_vector[i] = x[i];
- b_vector[i] = b[i];
- }
+ SmallMatrix A_dense{A};
+
+ SmallVector x_vector{x};
+ SmallVector b_vector{b};
this->solveLocalSystem(A_dense, x_vector, b_vector);
@@ -40,7 +37,7 @@ class LinearSolver
}
void solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b);
- void solveLocalSystem(const DenseMatrix<double>& A, Vector<double>& x, const Vector<const double>& b);
+ void solveLocalSystem(const SmallMatrix<double>& A, SmallVector<double>& x, const SmallVector<const double>& b);
LinearSolver();
LinearSolver(const LinearSolverOptions& options);
diff --git a/src/algebra/LinearSolverOptions.hpp b/src/algebra/LinearSolverOptions.hpp
index 014ee86d64bee4d2a128e820ab8876450a195369..940ce1ed63a2dbbf046022aee838c3f9438c1906 100644
--- a/src/algebra/LinearSolverOptions.hpp
+++ b/src/algebra/LinearSolverOptions.hpp
@@ -24,7 +24,7 @@ enum class LSMethod : int8_t
bicgstab2,
gmres,
lu,
- choleski,
+ cholesky,
//
LS__end
};
@@ -35,7 +35,7 @@ enum class LSPrecond : int8_t
//
none = LS__begin,
diagonal,
- incomplete_choleski,
+ incomplete_cholesky,
incomplete_LU,
amg,
//
@@ -77,8 +77,8 @@ name(const LSMethod method)
case LSMethod::lu: {
return "LU";
}
- case LSMethod::choleski: {
- return "Choleski";
+ case LSMethod::cholesky: {
+ return "Cholesky";
}
case LSMethod::LS__end: {
}
@@ -96,8 +96,8 @@ name(const LSPrecond precond)
case LSPrecond::diagonal: {
return "diagonal";
}
- case LSPrecond::incomplete_choleski: {
- return "ICholeski";
+ case LSPrecond::incomplete_cholesky: {
+ return "ICholesky";
}
case LSPrecond::incomplete_LU: {
return "ILU";
diff --git a/src/algebra/PETScUtils.cpp b/src/algebra/PETScUtils.cpp
index 63c5a68ac5f3a248ba7e5bc51bf381f5418dbe4a..29904103577a02f454b6c1fbb91f90fb637b8bec 100644
--- a/src/algebra/PETScUtils.cpp
+++ b/src/algebra/PETScUtils.cpp
@@ -12,24 +12,25 @@ PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const size_t nb_rows, const size_
MatSetUp(m_petscMat);
{
- Array<PetscInt> row_indices(nb_rows);
+ SmallArray<PetscInt> row_indices(nb_rows);
for (size_t i = 0; i < nb_rows; ++i) {
row_indices[i] = i;
}
- m_row_indices = row_indices;
+ m_small_row_indices = row_indices;
}
if (nb_rows == nb_columns) {
- m_column_indices = m_row_indices;
+ m_small_column_indices = m_small_row_indices;
} else {
- Array<PetscInt> column_indices(nb_columns);
+ SmallArray<PetscInt> column_indices(nb_columns);
for (size_t i = 0; i < nb_columns; ++i) {
column_indices[i] = i;
}
- m_column_indices = column_indices;
+ m_small_column_indices = column_indices;
}
- MatSetValuesBlocked(m_petscMat, nb_rows, &(m_row_indices[0]), nb_columns, &(m_column_indices[0]), A, INSERT_VALUES);
+ MatSetValuesBlocked(m_petscMat, nb_rows, &(m_small_row_indices[0]), nb_columns, &(m_small_column_indices[0]), A,
+ INSERT_VALUES);
MatAssemblyBegin(m_petscMat, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(m_petscMat, MAT_FINAL_ASSEMBLY);
@@ -42,9 +43,10 @@ PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const CRSMatrix<double, int>& A)
"index type conversion is not implemented yet");
static_assert(std::is_same_v<PetscScalar, typename std::decay_t<decltype(A)>::data_type>,
"value type conversion is not implemented yet");
- m_values = copy(A.values());
- m_row_indices = copy(A.rowMap());
- m_column_indices = copy(A.columnIndices());
+
+ m_values = A.values();
+ m_row_indices = A.rowMap();
+ m_column_indices = A.columnIndices();
MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.numberOfRows(), A.numberOfColumns(),
const_cast<PetscInt*>(&m_row_indices[0]), // These ugly const_cast
diff --git a/src/algebra/PETScUtils.hpp b/src/algebra/PETScUtils.hpp
index 3e7cff5b87b843f84575ed9b2af4a425e6ed4f5a..c4db9ee5910ef8e58333dbfb6dc20c8a1987bb46 100644
--- a/src/algebra/PETScUtils.hpp
+++ b/src/algebra/PETScUtils.hpp
@@ -6,7 +6,7 @@
#ifdef PUGS_HAS_PETSC
#include <algebra/CRSMatrix.hpp>
-#include <algebra/DenseMatrix.hpp>
+#include <algebra/SmallMatrix.hpp>
#include <algebra/TinyMatrix.hpp>
#include <petsc.h>
@@ -18,6 +18,8 @@ class PETScAijMatrixEmbedder
Array<const PetscScalar> m_values;
Array<const PetscInt> m_row_indices;
Array<const PetscInt> m_column_indices;
+ SmallArray<const PetscInt> m_small_row_indices;
+ SmallArray<const PetscInt> m_small_column_indices;
const size_t m_nb_rows;
const size_t m_nb_columns;
@@ -54,7 +56,7 @@ class PETScAijMatrixEmbedder
PETScAijMatrixEmbedder(const TinyMatrix<N>& A) : PETScAijMatrixEmbedder{N, N, &A(0, 0)}
{}
- PETScAijMatrixEmbedder(const DenseMatrix<double>& A)
+ PETScAijMatrixEmbedder(const SmallMatrix<double>& A)
: PETScAijMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), &A(0, 0)}
{}
diff --git a/src/algebra/DenseMatrix.hpp b/src/algebra/SmallMatrix.hpp
similarity index 62%
rename from src/algebra/DenseMatrix.hpp
rename to src/algebra/SmallMatrix.hpp
index 952342fedbd04246879356d67425275d471af2c9..d0d6ee329f02c41626900384c794d1cfa490089c 100644
--- a/src/algebra/DenseMatrix.hpp
+++ b/src/algebra/SmallMatrix.hpp
@@ -1,18 +1,18 @@
-#ifndef DENSE_MATRIX_HPP
-#define DENSE_MATRIX_HPP
+#ifndef SMALL_MATRIX_HPP
+#define SMALL_MATRIX_HPP
+#include <algebra/SmallVector.hpp>
#include <algebra/TinyMatrix.hpp>
-#include <algebra/Vector.hpp>
-#include <utils/Array.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
#include <utils/PugsUtils.hpp>
+#include <utils/SmallArray.hpp>
#include <utils/Types.hpp>
#include <iostream>
template <typename DataType>
-class DenseMatrix // LCOV_EXCL_LINE
+class [[nodiscard]] SmallMatrix // LCOV_EXCL_LINE
{
public:
using data_type = DataType;
@@ -21,32 +21,29 @@ class DenseMatrix // LCOV_EXCL_LINE
private:
size_t m_nb_rows;
size_t m_nb_columns;
- Array<DataType> m_values;
+ SmallArray<DataType> m_values;
static_assert(std::is_same_v<typename decltype(m_values)::index_type, index_type>);
- static_assert(std::is_arithmetic_v<DataType>, "Dense matrices expect arithmetic data");
+ static_assert(std::is_arithmetic_v<DataType>, "Small matrices expect arithmetic data");
// Allows const version to access our data
- friend DenseMatrix<std::add_const_t<DataType>>;
+ friend SmallMatrix<std::add_const_t<DataType>>;
public:
PUGS_INLINE
- bool
- isSquare() const noexcept
+ bool isSquare() const noexcept
{
return m_nb_rows == m_nb_columns;
}
- friend PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
- copy(const DenseMatrix& A) noexcept
+ friend PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> copy(const SmallMatrix& A) noexcept
{
return {A.m_nb_rows, A.m_nb_columns, copy(A.m_values)};
}
- friend PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
- transpose(const DenseMatrix& A)
+ friend PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> transpose(const SmallMatrix& A)
{
- DenseMatrix<std::remove_const_t<DataType>> A_transpose{A.m_nb_columns, A.m_nb_rows};
+ SmallMatrix<std::remove_const_t<DataType>> A_transpose{A.m_nb_columns, A.m_nb_rows};
for (size_t i = 0; i < A.m_nb_rows; ++i) {
for (size_t j = 0; j < A.m_nb_columns; ++j) {
A_transpose(j, i) = A(i, j);
@@ -55,22 +52,20 @@ class DenseMatrix // LCOV_EXCL_LINE
return A_transpose;
}
- friend PUGS_INLINE DenseMatrix
- operator*(const DataType& a, const DenseMatrix& A)
+ friend PUGS_INLINE SmallMatrix operator*(const DataType& a, const SmallMatrix& A)
{
- DenseMatrix<std::remove_const_t<DataType>> aA = copy(A);
+ SmallMatrix<std::remove_const_t<DataType>> aA = copy(A);
return aA *= a;
}
template <typename DataType2>
- PUGS_INLINE Vector<std::remove_const_t<DataType>>
- operator*(const Vector<DataType2>& x) const
+ PUGS_INLINE SmallVector<std::remove_const_t<DataType>> operator*(const SmallVector<DataType2>& x) const
{
static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
"incompatible data types");
Assert(m_nb_columns == x.size(), "cannot compute matrix-vector product: incompatible sizes");
- const DenseMatrix& A = *this;
- Vector<std::remove_const_t<DataType>> Ax{m_nb_rows};
+ const SmallMatrix& A = *this;
+ SmallVector<std::remove_const_t<DataType>> Ax{m_nb_rows};
for (size_t i = 0; i < m_nb_rows; ++i) {
std::remove_const_t<DataType> Axi = A(i, 0) * x[0];
for (size_t j = 1; j < m_nb_columns; ++j) {
@@ -82,14 +77,13 @@ class DenseMatrix // LCOV_EXCL_LINE
}
template <typename DataType2>
- PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
- operator*(const DenseMatrix<DataType2>& B) const
+ PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> operator*(const SmallMatrix<DataType2>& B) const
{
static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
"incompatible data types");
Assert(m_nb_columns == B.numberOfRows(), "cannot compute matrix product: incompatible sizes");
- const DenseMatrix& A = *this;
- DenseMatrix<std::remove_const_t<DataType>> AB{m_nb_rows, B.numberOfColumns()};
+ const SmallMatrix& A = *this;
+ SmallMatrix<std::remove_const_t<DataType>> AB{m_nb_rows, B.numberOfColumns()};
for (size_t i = 0; i < m_nb_rows; ++i) {
for (size_t j = 0; j < B.numberOfColumns(); ++j) {
@@ -104,16 +98,14 @@ class DenseMatrix // LCOV_EXCL_LINE
}
template <typename DataType2>
- PUGS_INLINE DenseMatrix&
- operator/=(const DataType2& a)
+ PUGS_INLINE SmallMatrix& operator/=(const DataType2& a)
{
const auto inv_a = 1. / a;
return (*this) *= inv_a;
}
template <typename DataType2>
- PUGS_INLINE DenseMatrix&
- operator*=(const DataType2& a)
+ PUGS_INLINE SmallMatrix& operator*=(const DataType2& a)
{
parallel_for(
m_values.size(), PUGS_LAMBDA(index_type i) { m_values[i] *= a; });
@@ -121,8 +113,7 @@ class DenseMatrix // LCOV_EXCL_LINE
}
template <typename DataType2>
- PUGS_INLINE DenseMatrix&
- operator-=(const DenseMatrix<DataType2>& B)
+ PUGS_INLINE SmallMatrix& operator-=(const SmallMatrix<DataType2>& B)
{
static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
"incompatible data types");
@@ -135,8 +126,7 @@ class DenseMatrix // LCOV_EXCL_LINE
}
template <typename DataType2>
- PUGS_INLINE DenseMatrix&
- operator+=(const DenseMatrix<DataType2>& B)
+ PUGS_INLINE SmallMatrix& operator+=(const SmallMatrix<DataType2>& B)
{
static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
"incompatible data types");
@@ -148,15 +138,14 @@ class DenseMatrix // LCOV_EXCL_LINE
}
template <typename DataType2>
- PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
- operator+(const DenseMatrix<DataType2>& B) const
+ PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> operator+(const SmallMatrix<DataType2>& B) const
{
static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
"incompatible data types");
Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()),
"cannot compute matrix sum: incompatible sizes");
- DenseMatrix<std::remove_const_t<DataType>> sum{B.numberOfRows(), B.numberOfColumns()};
+ SmallMatrix<std::remove_const_t<DataType>> sum{B.numberOfRows(), B.numberOfColumns()};
parallel_for(
m_values.size(), PUGS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] + B.m_values[i]; });
@@ -165,15 +154,14 @@ class DenseMatrix // LCOV_EXCL_LINE
}
template <typename DataType2>
- PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
- operator-(const DenseMatrix<DataType2>& B) const
+ PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> operator-(const SmallMatrix<DataType2>& B) const
{
static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
"incompatible data types");
Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()),
"cannot compute matrix difference: incompatible sizes");
- DenseMatrix<std::remove_const_t<DataType>> difference{B.numberOfRows(), B.numberOfColumns()};
+ SmallMatrix<std::remove_const_t<DataType>> difference{B.numberOfRows(), B.numberOfColumns()};
parallel_for(
m_values.size(), PUGS_LAMBDA(index_type i) { difference.m_values[i] = m_values[i] - B.m_values[i]; });
@@ -182,40 +170,36 @@ class DenseMatrix // LCOV_EXCL_LINE
}
PUGS_INLINE
- DataType&
- operator()(index_type i, index_type j) const noexcept(NO_ASSERT)
+ DataType& operator()(index_type i, index_type j) const noexcept(NO_ASSERT)
{
Assert(i < m_nb_rows and j < m_nb_columns, "cannot access element: invalid indices");
return m_values[i * m_nb_columns + j];
}
PUGS_INLINE
- size_t
- numberOfRows() const noexcept
+ size_t numberOfRows() const noexcept
{
return m_nb_rows;
}
PUGS_INLINE
- size_t
- numberOfColumns() const noexcept
+ size_t numberOfColumns() const noexcept
{
return m_nb_columns;
}
- PUGS_INLINE void
- fill(const DataType& value) noexcept
+ PUGS_INLINE void fill(const DataType& value) noexcept
{
m_values.fill(value);
}
- PUGS_INLINE DenseMatrix& operator=(ZeroType) noexcept
+ PUGS_INLINE SmallMatrix& operator=(ZeroType) noexcept
{
m_values.fill(0);
return *this;
}
- PUGS_INLINE DenseMatrix& operator=(IdentityType) noexcept(NO_ASSERT)
+ PUGS_INLINE SmallMatrix& operator=(IdentityType) noexcept(NO_ASSERT)
{
Assert(m_nb_rows == m_nb_columns, "identity must be a square matrix");
@@ -226,15 +210,14 @@ class DenseMatrix // LCOV_EXCL_LINE
}
template <typename DataType2>
- PUGS_INLINE DenseMatrix&
- operator=(const DenseMatrix<DataType2>& A) noexcept
+ PUGS_INLINE SmallMatrix& operator=(const SmallMatrix<DataType2>& A) noexcept
{
// ensures that DataType is the same as source DataType2
static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
- "Cannot assign DenseMatrix of different type");
+ "Cannot assign SmallMatrix of different type");
// ensures that const is not lost through copy
static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
- "Cannot assign DenseMatrix of const to DenseMatrix of non-const");
+ "Cannot assign SmallMatrix of const to SmallMatrix of non-const");
m_nb_rows = A.m_nb_rows;
m_nb_columns = A.m_nb_columns;
@@ -243,13 +226,12 @@ class DenseMatrix // LCOV_EXCL_LINE
}
PUGS_INLINE
- DenseMatrix& operator=(const DenseMatrix&) = default;
+ SmallMatrix& operator=(const SmallMatrix&) = default;
PUGS_INLINE
- DenseMatrix& operator=(DenseMatrix&&) = default;
+ SmallMatrix& operator=(SmallMatrix&&) = default;
- friend std::ostream&
- operator<<(std::ostream& os, const DenseMatrix& A)
+ friend std::ostream& operator<<(std::ostream& os, const SmallMatrix& A)
{
for (size_t i = 0; i < A.numberOfRows(); ++i) {
os << i << '|';
@@ -262,51 +244,42 @@ class DenseMatrix // LCOV_EXCL_LINE
}
template <typename DataType2>
- DenseMatrix(const DenseMatrix<DataType2>& A)
+ SmallMatrix(const SmallMatrix<DataType2>& A)
{
// ensures that DataType is the same as source DataType2
static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
- "Cannot assign DenseMatrix of different type");
+ "Cannot assign SmallMatrix of different type");
// ensures that const is not lost through copy
static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
- "Cannot assign DenseMatrix of const to DenseMatrix of non-const");
+ "Cannot assign SmallMatrix of const to SmallMatrix of non-const");
this->operator=(A);
}
- DenseMatrix(const DenseMatrix&) = default;
+ SmallMatrix(const SmallMatrix&) = default;
- DenseMatrix(DenseMatrix&&) = default;
+ SmallMatrix(SmallMatrix &&) = default;
- explicit DenseMatrix(size_t nb_rows, size_t nb_columns) noexcept
+ explicit SmallMatrix(size_t nb_rows, size_t nb_columns) noexcept
: m_nb_rows{nb_rows}, m_nb_columns{nb_columns}, m_values{nb_rows * nb_columns}
{}
- explicit DenseMatrix(size_t nb_rows) noexcept : m_nb_rows{nb_rows}, m_nb_columns{nb_rows}, m_values{nb_rows * nb_rows}
+ explicit SmallMatrix(size_t nb_rows) noexcept : m_nb_rows{nb_rows}, m_nb_columns{nb_rows}, m_values{nb_rows * nb_rows}
{}
- template <size_t N>
- explicit DenseMatrix(const TinyMatrix<N, DataType>& M) noexcept : m_nb_rows{N}, m_nb_columns{N}, m_values{N * N}
+ template <size_t M, size_t N>
+ explicit SmallMatrix(const TinyMatrix<M, N, DataType>& A) noexcept : m_nb_rows{M}, m_nb_columns{N}, m_values{M * N}
{
- parallel_for(
- N, PUGS_LAMBDA(const index_type i) {
- for (size_t j = 0; j < N; ++j) {
- m_values[i * N + j] = M(i, j);
- }
- });
+ for (size_t i = 0; i < M; ++i) {
+ for (size_t j = 0; j < N; ++j) {
+ m_values[i * N + j] = A(i, j);
+ }
+ }
}
- DenseMatrix() noexcept : m_nb_rows{0}, m_nb_columns{0} {}
+ SmallMatrix() noexcept : m_nb_rows{0}, m_nb_columns{0} {}
- private:
- DenseMatrix(size_t nb_rows, size_t nb_columns, const Array<DataType> values) noexcept(NO_ASSERT)
- : m_nb_rows{nb_rows}, m_nb_columns{nb_columns}, m_values{values}
- {
- Assert(m_values.size() == m_nb_rows * m_nb_columns, "cannot create matrix: incorrect number of values");
- }
-
- public:
- ~DenseMatrix() = default;
+ ~SmallMatrix() = default;
};
-#endif // DENSE_MATRIX_HPP
+#endif // SMALL_MATRIX_HPP
diff --git a/src/algebra/SmallVector.hpp b/src/algebra/SmallVector.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d20bb48f4cdb700f61a0881898fb5273217693aa
--- /dev/null
+++ b/src/algebra/SmallVector.hpp
@@ -0,0 +1,216 @@
+#ifndef SMALL_VECTOR_HPP
+#define SMALL_VECTOR_HPP
+
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+#include <utils/PugsUtils.hpp>
+#include <utils/SmallArray.hpp>
+#include <utils/Types.hpp>
+
+template <typename DataType>
+class SmallVector // LCOV_EXCL_LINE
+{
+ public:
+ using data_type = DataType;
+ using index_type = size_t;
+
+ private:
+ SmallArray<DataType> m_values;
+
+ static_assert(std::is_same_v<typename decltype(m_values)::index_type, index_type>);
+ static_assert(std::is_arithmetic_v<DataType>, "SmallVector is only defined for arithmetic values");
+
+ // Allows const version to access our data
+ friend SmallVector<std::add_const_t<DataType>>;
+
+ public:
+ friend std::ostream&
+ operator<<(std::ostream& os, const SmallVector& x)
+ {
+ return os << x.m_values;
+ }
+
+ friend SmallVector<std::remove_const_t<DataType>>
+ copy(const SmallVector& x)
+ {
+ auto values = copy(x.m_values);
+
+ SmallVector<std::remove_const_t<DataType>> x_copy{0};
+ x_copy.m_values = values;
+ return x_copy;
+ }
+
+ friend SmallVector operator*(const DataType& a, const SmallVector& x)
+ {
+ SmallVector<std::remove_const_t<DataType>> y = copy(x);
+ return y *= a;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE friend auto
+ dot(const SmallVector& x, const SmallVector<DataType2>& y)
+ {
+ Assert(x.size() == y.size(), "cannot compute dot product: incompatible vector sizes");
+ // Quite ugly, TODO: fix me in C++20
+ auto promoted = [] {
+ DataType a{0};
+ DataType2 b{0};
+ return a * b;
+ }();
+
+ decltype(promoted) sum = 0;
+
+ // Does not use parallel_for to preserve sum order
+ for (index_type i = 0; i < x.size(); ++i) {
+ sum += x[i] * y[i];
+ }
+
+ return sum;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE SmallVector&
+ operator/=(const DataType2& a)
+ {
+ const auto inv_a = 1. / a;
+ return (*this) *= inv_a;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE SmallVector&
+ operator*=(const DataType2& a)
+ {
+ parallel_for(
+ this->size(), PUGS_LAMBDA(index_type i) { m_values[i] *= a; });
+ return *this;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE SmallVector&
+ operator-=(const SmallVector<DataType2>& y)
+ {
+ Assert(this->size() == y.size(), "cannot substract vector: incompatible sizes");
+
+ parallel_for(
+ this->size(), PUGS_LAMBDA(index_type i) { m_values[i] -= y[i]; });
+
+ return *this;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE SmallVector&
+ operator+=(const SmallVector<DataType2>& y)
+ {
+ Assert(this->size() == y.size(), "cannot add vector: incompatible sizes");
+
+ parallel_for(
+ this->size(), PUGS_LAMBDA(index_type i) { m_values[i] += y[i]; });
+
+ return *this;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE SmallVector<std::remove_const_t<DataType>>
+ operator+(const SmallVector<DataType2>& y) const
+ {
+ Assert(this->size() == y.size(), "cannot compute vector sum: incompatible sizes");
+ SmallVector<std::remove_const_t<DataType>> sum{y.size()};
+
+ parallel_for(
+ this->size(), PUGS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] + y[i]; });
+
+ return sum;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE SmallVector<std::remove_const_t<DataType>>
+ operator-(const SmallVector<DataType2>& y) const
+ {
+ Assert(this->size() == y.size(), "cannot compute vector difference: incompatible sizes");
+ SmallVector<std::remove_const_t<DataType>> sum{y.size()};
+
+ parallel_for(
+ this->size(), PUGS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] - y[i]; });
+
+ return sum;
+ }
+
+ PUGS_INLINE
+ DataType& operator[](index_type i) const noexcept(NO_ASSERT)
+ {
+ return m_values[i];
+ }
+
+ PUGS_INLINE
+ size_t
+ size() const noexcept
+ {
+ return m_values.size();
+ }
+
+ PUGS_INLINE SmallVector&
+ fill(const DataType& value) noexcept
+ {
+ m_values.fill(value);
+ return *this;
+ }
+
+ PUGS_INLINE SmallVector& operator=(const DataType&) = delete;
+
+ PUGS_INLINE SmallVector&
+ operator=(const ZeroType&) noexcept
+ {
+ m_values.fill(0);
+ return *this;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE SmallVector&
+ operator=(const SmallVector<DataType2>& x)
+ {
+ // ensures that DataType is the same as source DataType2
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot assign SmallVector of different type");
+ // ensures that const is not lost through copy. If this point is
+ // reached data types only differ through their constness
+ static_assert((not std::is_const_v<DataType2>), "Cannot assign SmallVector of const to SmallVector of non-const");
+
+ m_values = x.m_values;
+ return *this;
+ }
+
+ PUGS_INLINE
+ SmallVector& operator=(const SmallVector&) = default;
+
+ PUGS_INLINE
+ SmallVector& operator=(SmallVector&& x) = default;
+
+ template <typename DataType2>
+ SmallVector(const SmallVector<DataType2>& x)
+ {
+ // ensures that DataType is the same as source DataType2
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot assign SmallVector of different type");
+ // ensures that const is not lost through copy
+ static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
+ "Cannot assign SmallVector of const to SmallVector of non-const");
+
+ m_values = x.m_values;
+ }
+
+ template <size_t N>
+ explicit SmallVector(const TinyVector<N, DataType>& x) : m_values{N}
+ {
+ std::copy(&x[0], &x[0] + N, &m_values[0]);
+ }
+
+ SmallVector(const SmallVector&) = default;
+ SmallVector(SmallVector&&) = default;
+
+ SmallVector(size_t size) : m_values{size} {}
+
+ SmallVector() = default;
+ ~SmallVector() = default;
+};
+
+#endif // SMALL_VECTOR_HPP
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 98030dd447ee6cac14b6de0f3c45d55b68cdf7ba..0d3b9c0e5fef46f85fa576b68f940708c8ed529d 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -10,15 +10,14 @@
#include <iostream>
-template <size_t N, typename T = double>
+template <size_t M, size_t N = M, typename T = double>
class [[nodiscard]] TinyMatrix
{
public:
using data_type = T;
private:
- T m_values[N * N];
- static_assert((N > 0), "TinyMatrix size must be strictly positive");
+ T m_values[M * N];
PUGS_FORCEINLINE
constexpr size_t _index(size_t i, size_t j) const noexcept // LCOV_EXCL_LINE (due to forced inline)
@@ -29,7 +28,7 @@ class [[nodiscard]] TinyMatrix
template <typename... Args>
PUGS_FORCEINLINE constexpr void _unpackVariadicInput(const T& t, Args&&... args) noexcept
{
- m_values[N * N - 1 - sizeof...(args)] = t;
+ m_values[M * N - 1 - sizeof...(args)] = t;
if constexpr (sizeof...(args) > 0) {
this->_unpackVariadicInput(std::forward<Args>(args)...);
}
@@ -39,13 +38,25 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr bool isSquare() const noexcept
{
- return true;
+ return M == N;
+ }
+
+ PUGS_INLINE
+ constexpr friend TinyMatrix<N, M, T> transpose(const TinyMatrix& A)
+ {
+ TinyMatrix<N, M, T> tA;
+ for (size_t i = 0; i < M; ++i) {
+ for (size_t j = 0; j < N; ++j) {
+ tA(j, i) = A(i, j);
+ }
+ }
+ return tA;
}
PUGS_INLINE
constexpr size_t dimension() const
{
- return N * N;
+ return M * N;
}
PUGS_INLINE
@@ -57,7 +68,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr size_t numberOfRows() const
{
- return N;
+ return M;
}
PUGS_INLINE
@@ -70,7 +81,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix operator-() const
{
TinyMatrix opposite;
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
opposite.m_values[i] = -m_values[i];
}
return opposite;
@@ -92,19 +103,19 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr TinyMatrix& operator*=(const T& t)
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] *= t;
}
return *this;
}
- PUGS_INLINE
- constexpr TinyMatrix operator*(const TinyMatrix& B) const
+ template <size_t P>
+ PUGS_INLINE constexpr TinyMatrix<M, P, T> operator*(const TinyMatrix<N, P, T>& B) const
{
const TinyMatrix& A = *this;
- TinyMatrix AB;
- for (size_t i = 0; i < N; ++i) {
- for (size_t j = 0; j < N; ++j) {
+ TinyMatrix<M, P, T> AB;
+ for (size_t i = 0; i < M; ++i) {
+ for (size_t j = 0; j < P; ++j) {
T sum = A(i, 0) * B(0, j);
for (size_t k = 1; k < N; ++k) {
sum += A(i, k) * B(k, j);
@@ -116,11 +127,11 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyVector<N, T> operator*(const TinyVector<N, T>& x) const
+ constexpr TinyVector<M, T> operator*(const TinyVector<N, T>& x) const
{
const TinyMatrix& A = *this;
- TinyVector<N, T> Ax;
- for (size_t i = 0; i < N; ++i) {
+ TinyVector<M, T> Ax;
+ for (size_t i = 0; i < M; ++i) {
T sum = A(i, 0) * x[0];
for (size_t j = 1; j < N; ++j) {
sum += A(i, j) * x[j];
@@ -134,7 +145,7 @@ class [[nodiscard]] TinyMatrix
constexpr friend std::ostream& operator<<(std::ostream& os, const TinyMatrix& A)
{
os << '[';
- for (size_t i = 0; i < N; ++i) {
+ for (size_t i = 0; i < M; ++i) {
os << '(' << NaNHelper(A(i, 0));
for (size_t j = 1; j < N; ++j) {
os << ',' << NaNHelper(A(i, j));
@@ -149,7 +160,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr bool operator==(const TinyMatrix& A) const
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
if (m_values[i] != A.m_values[i])
return false;
}
@@ -166,7 +177,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix operator+(const TinyMatrix& A) const
{
TinyMatrix sum;
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
sum.m_values[i] = m_values[i] + A.m_values[i];
}
return sum;
@@ -175,9 +186,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr TinyMatrix operator+(TinyMatrix&& A) const
{
- for (size_t i = 0; i < N * N; ++i) {
- A.m_values[i] += m_values[i];
- }
+ A += *this;
return std::move(A);
}
@@ -185,7 +194,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix operator-(const TinyMatrix& A) const
{
TinyMatrix difference;
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
difference.m_values[i] = m_values[i] - A.m_values[i];
}
return difference;
@@ -194,7 +203,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr TinyMatrix operator-(TinyMatrix&& A) const
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
A.m_values[i] = m_values[i] - A.m_values[i];
}
return std::move(A);
@@ -203,7 +212,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr TinyMatrix& operator+=(const TinyMatrix& A)
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] += A.m_values[i];
}
return *this;
@@ -212,7 +221,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr void operator+=(const volatile TinyMatrix& A) volatile
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] += A.m_values[i];
}
}
@@ -220,7 +229,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr TinyMatrix& operator-=(const TinyMatrix& A)
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] -= A.m_values[i];
}
return *this;
@@ -229,14 +238,14 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr T& operator()(size_t i, size_t j) noexcept(NO_ASSERT)
{
- Assert((i < N) and (j < N));
+ Assert((i < M) and (j < N));
return m_values[_index(i, j)];
}
PUGS_INLINE
constexpr const T& operator()(size_t i, size_t j) const noexcept(NO_ASSERT)
{
- Assert((i < N) and (j < N));
+ Assert((i < M) and (j < N));
return m_values[_index(i, j)];
}
@@ -244,7 +253,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix& operator=(ZeroType) noexcept
{
static_assert(std::is_arithmetic<T>(), "Cannot assign 'zero' value for non-arithmetic types");
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] = 0;
}
return *this;
@@ -254,7 +263,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix& operator=(IdentityType) noexcept
{
static_assert(std::is_arithmetic<T>(), "Cannot assign 'identity' value for non-arithmetic types");
- for (size_t i = 0; i < N; ++i) {
+ for (size_t i = 0; i < M; ++i) {
for (size_t j = 0; j < N; ++j) {
m_values[_index(i, j)] = (i == j) ? 1 : 0;
}
@@ -271,7 +280,7 @@ class [[nodiscard]] TinyMatrix
template <typename... Args>
PUGS_INLINE constexpr TinyMatrix(const T& t, Args&&... args) noexcept
{
- static_assert(sizeof...(args) == N * N - 1, "wrong number of parameters");
+ static_assert(sizeof...(args) == M * N - 1, "wrong number of parameters");
this->_unpackVariadicInput(t, std::forward<Args>(args)...);
}
@@ -281,7 +290,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix() noexcept
{
#ifndef NDEBUG
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] = invalid_data_v<T>;
}
#endif // NDEBUG
@@ -292,7 +301,7 @@ class [[nodiscard]] TinyMatrix
{
static_assert(std::is_arithmetic<T>(), "Cannot construct from 'zero' value "
"for non-arithmetic types");
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] = 0;
}
}
@@ -302,7 +311,7 @@ class [[nodiscard]] TinyMatrix
{
static_assert(std::is_arithmetic<T>(), "Cannot construct from 'identity' "
"value for non-arithmetic types");
- for (size_t i = 0; i < N; ++i) {
+ for (size_t i = 0; i < M; ++i) {
for (size_t j = 0; j < N; ++j) {
m_values[_index(i, j)] = (i == j) ? 1 : 0;
}
@@ -319,12 +328,12 @@ class [[nodiscard]] TinyMatrix
~TinyMatrix() = default;
};
-template <size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<N, T>
-tensorProduct(const TinyVector<N, T>& x, const TinyVector<N, T>& y)
+template <size_t M, size_t N, typename T>
+PUGS_INLINE constexpr TinyMatrix<M, N, T>
+tensorProduct(const TinyVector<M, T>& x, const TinyVector<N, T>& y)
{
- TinyMatrix<N, T> A;
- for (size_t i = 0; i < N; ++i) {
+ TinyMatrix<M, N, T> A;
+ for (size_t i = 0; i < M; ++i) {
for (size_t j = 0; j < N; ++j) {
A(i, j) = x[i] * y[j];
}
@@ -334,16 +343,17 @@ tensorProduct(const TinyVector<N, T>& x, const TinyVector<N, T>& y)
template <size_t N, typename T>
PUGS_INLINE constexpr T
-det(const TinyMatrix<N, T>& A)
+det(const TinyMatrix<N, N, T>& A)
{
- static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
- static_assert(std::is_floating_point<T>::value, "determinent for arbitrary dimension N is defined for floating "
+ static_assert(std::is_arithmetic<T>::value, "determinant is not defined for non-arithmetic types");
+ static_assert(std::is_floating_point<T>::value, "determinant for arbitrary dimension N is defined for floating "
"point types only");
- TinyMatrix<N, T> M = A;
+ TinyMatrix B = A;
TinyVector<N, size_t> index;
- for (size_t i = 0; i < N; ++i)
+ for (size_t i = 0; i < N; ++i) {
index[i] = i;
+ }
T determinent = 1;
for (size_t i = 0; i < N; ++i) {
@@ -351,7 +361,7 @@ det(const TinyMatrix<N, T>& A)
size_t l = j;
const size_t J = index[j];
for (size_t k = j + 1; k < N; ++k) {
- if (std::abs(M(index[k], i)) > std::abs(M(J, i))) {
+ if (std::abs(B(index[k], i)) > std::abs(B(J, i))) {
l = k;
}
}
@@ -361,14 +371,14 @@ det(const TinyMatrix<N, T>& A)
}
}
const size_t I = index[i];
- const T Mii = M(I, i);
- if (Mii != 0) {
- const T inv_Mii = 1. / M(I, i);
+ const T Bii = B(I, i);
+ if (Bii != 0) {
+ const T inv_Bii = 1. / B(I, i);
for (size_t k = i + 1; k < N; ++k) {
const size_t K = index[k];
- const T factor = M(K, i) * inv_Mii;
+ const T factor = B(K, i) * inv_Bii;
for (size_t l = i + 1; l < N; ++l) {
- M(K, l) -= factor * M(I, l);
+ B(K, l) -= factor * B(I, l);
}
}
} else {
@@ -377,14 +387,14 @@ det(const TinyMatrix<N, T>& A)
}
for (size_t i = 0; i < N; ++i) {
- determinent *= M(index[i], i);
+ determinent *= B(index[i], i);
}
return determinent;
}
template <typename T>
PUGS_INLINE constexpr T
-det(const TinyMatrix<1, T>& A)
+det(const TinyMatrix<1, 1, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
return A(0, 0);
@@ -392,7 +402,7 @@ det(const TinyMatrix<1, T>& A)
template <typename T>
PUGS_INLINE constexpr T
-det(const TinyMatrix<2, T>& A)
+det(const TinyMatrix<2, 2, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
return A(0, 0) * A(1, 1) - A(1, 0) * A(0, 1);
@@ -400,56 +410,56 @@ det(const TinyMatrix<2, T>& A)
template <typename T>
PUGS_INLINE constexpr T
-det(const TinyMatrix<3, T>& A)
+det(const TinyMatrix<3, 3, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
return A(0, 0) * (A(1, 1) * A(2, 2) - A(2, 1) * A(1, 2)) - A(1, 0) * (A(0, 1) * A(2, 2) - A(2, 1) * A(0, 2)) +
A(2, 0) * (A(0, 1) * A(1, 2) - A(1, 1) * A(0, 2));
}
-template <size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<N - 1, T>
-getMinor(const TinyMatrix<N, T>& A, size_t I, size_t J)
+template <size_t M, size_t N, typename T>
+PUGS_INLINE constexpr TinyMatrix<M - 1, N - 1, T>
+getMinor(const TinyMatrix<M, N, T>& A, size_t I, size_t J)
{
- static_assert(N >= 2, "minor calculation requires at least 2x2 matrices");
- Assert((I < N) and (J < N));
- TinyMatrix<N - 1, T> M;
+ static_assert(M >= 2 and N >= 2, "minor calculation requires at least 2x2 matrices");
+ Assert((I < M) and (J < N));
+ TinyMatrix<M - 1, N - 1, T> minor;
for (size_t i = 0; i < I; ++i) {
for (size_t j = 0; j < J; ++j) {
- M(i, j) = A(i, j);
+ minor(i, j) = A(i, j);
}
for (size_t j = J + 1; j < N; ++j) {
- M(i, j - 1) = A(i, j);
+ minor(i, j - 1) = A(i, j);
}
}
- for (size_t i = I + 1; i < N; ++i) {
+ for (size_t i = I + 1; i < M; ++i) {
for (size_t j = 0; j < J; ++j) {
- M(i - 1, j) = A(i, j);
+ minor(i - 1, j) = A(i, j);
}
for (size_t j = J + 1; j < N; ++j) {
- M(i - 1, j - 1) = A(i, j);
+ minor(i - 1, j - 1) = A(i, j);
}
}
- return M;
+ return minor;
}
template <size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<N, T> inverse(const TinyMatrix<N, T>& A);
+PUGS_INLINE constexpr TinyMatrix<N, N, T> inverse(const TinyMatrix<N, N, T>& A);
template <typename T>
-PUGS_INLINE constexpr TinyMatrix<1, T>
-inverse(const TinyMatrix<1, T>& A)
+PUGS_INLINE constexpr TinyMatrix<1, 1, T>
+inverse(const TinyMatrix<1, 1, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
- TinyMatrix<1, T> A_1(1. / A(0, 0));
+ TinyMatrix<1, 1, T> A_1(1. / A(0, 0));
return A_1;
}
template <size_t N, typename T>
PUGS_INLINE constexpr T
-cofactor(const TinyMatrix<N, T>& A, size_t i, size_t j)
+cofactor(const TinyMatrix<N, N, T>& A, size_t i, size_t j)
{
static_assert(std::is_arithmetic<T>::value, "cofactor is not defined for non-arithmetic types");
const T sign = ((i + j) % 2) ? -1 : 1;
@@ -458,8 +468,8 @@ cofactor(const TinyMatrix<N, T>& A, size_t i, size_t j)
}
template <typename T>
-PUGS_INLINE constexpr TinyMatrix<2, T>
-inverse(const TinyMatrix<2, T>& A)
+PUGS_INLINE constexpr TinyMatrix<2, 2, T>
+inverse(const TinyMatrix<2, 2, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
@@ -467,22 +477,22 @@ inverse(const TinyMatrix<2, T>& A)
const T determinent = det(A);
const T inv_determinent = 1. / determinent;
- TinyMatrix<2, T> A_cofactors_T(A(1, 1), -A(0, 1), -A(1, 0), A(0, 0));
+ TinyMatrix<2, 2, T> A_cofactors_T(A(1, 1), -A(0, 1), -A(1, 0), A(0, 0));
return A_cofactors_T *= inv_determinent;
}
template <typename T>
-PUGS_INLINE constexpr TinyMatrix<3, T>
-inverse(const TinyMatrix<3, T>& A)
+PUGS_INLINE constexpr TinyMatrix<3, 3, T>
+inverse(const TinyMatrix<3, 3, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
const T determinent = det(A);
- TinyMatrix<3, T> A_cofactors_T(cofactor(A, 0, 0), cofactor(A, 1, 0), cofactor(A, 2, 0), cofactor(A, 0, 1),
- cofactor(A, 1, 1), cofactor(A, 2, 1), cofactor(A, 0, 2), cofactor(A, 1, 2),
- cofactor(A, 2, 2));
+ TinyMatrix<3, 3, T> A_cofactors_T(cofactor(A, 0, 0), cofactor(A, 1, 0), cofactor(A, 2, 0), cofactor(A, 0, 1),
+ cofactor(A, 1, 1), cofactor(A, 2, 1), cofactor(A, 0, 2), cofactor(A, 1, 2),
+ cofactor(A, 2, 2));
return A_cofactors_T *= 1. / determinent;
}
diff --git a/src/language/PEGGrammar.hpp b/src/language/PEGGrammar.hpp
index a3ad0342839eb5cc5ee6693e4a49e22d99cf41d2..723d9a710f277aff551508ef849632c24e727575 100644
--- a/src/language/PEGGrammar.hpp
+++ b/src/language/PEGGrammar.hpp
@@ -126,11 +126,7 @@ struct BREAK : seq < break_kw, ignored > {};
struct continue_kw : TAO_PEGTL_KEYWORD("continue") {};
struct CONTINUE : seq < continue_kw, ignored > {};
-struct cout_kw : TAO_PEGTL_KEYWORD("cout") {};
-struct cerr_kw : TAO_PEGTL_KEYWORD("cerr") {};
-struct clog_kw : TAO_PEGTL_KEYWORD("clog") {};
-
-struct keyword : sor < basic_type, import_kw, true_kw, false_kw, let_kw, do_kw, while_kw, for_kw, if_kw, else_kw, and_kw, or_kw, xor_kw, not_kw, break_kw, continue_kw, cout_kw, cerr_kw, clog_kw > {};
+struct keyword : sor < basic_type, import_kw, true_kw, false_kw, let_kw, do_kw, while_kw, for_kw, if_kw, else_kw, and_kw, or_kw, xor_kw, not_kw, break_kw, continue_kw> {};
struct identifier_minus_keyword : minus< identifier, keyword > {};
@@ -218,7 +214,9 @@ struct product : list_must< unary_expression, sor< multiply_op, divide_op > > {}
struct sum : list_must< product, sor< plus_op, minus_op > > {};
-struct compare : list_must<sum, sor< lesser_or_eq_op, greater_or_eq_op, lesser_op, greater_op > >{};
+struct shift : list_must<sum, sor< shift_left_op, shift_right_op > >{};
+
+struct compare : list_must<shift, sor< lesser_or_eq_op, greater_or_eq_op, lesser_op, greater_op > >{};
struct equality : list_must< compare, sor< eqeq_op, not_eq_op > >{};
@@ -289,15 +287,11 @@ struct for_statement_block;
struct for_statement : if_must< FOR, open_parent, for_init, SEMICOL, for_test, SEMICOL, for_post, close_parent, for_statement_block >{};
-struct ostream_object : seq< sor< cout_kw, cerr_kw, clog_kw >, ignored >{};
-struct ostream_statement : seq< ostream_object, star< if_must< shift_left_op, expression, ignored > > >{};
-
struct instruction
: sor<if_statement,
if_must<do_while_statement, semicol>,
while_statement,
for_statement,
- if_must< ostream_statement, semicol >,
if_must< BREAK, semicol >,
if_must< CONTINUE, semicol >,
block,
diff --git a/src/language/ast/ASTBuilder.cpp b/src/language/ast/ASTBuilder.cpp
index 0c433bbfd21704702ed87012d5956519a0f93bde..a4dc8c822bab53a67be975b6444eb9dffa9bb6be 100644
--- a/src/language/ast/ASTBuilder.cpp
+++ b/src/language/ast/ASTBuilder.cpp
@@ -214,20 +214,6 @@ struct ASTBuilder::simplify_for_post : TAO_PEGTL_NAMESPACE::parse_tree::apply<AS
}
};
-struct ASTBuilder::simplify_stream_statement
- : TAO_PEGTL_NAMESPACE::parse_tree::apply<ASTBuilder::simplify_stream_statement>
-{
- template <typename... States>
- static void
- transform(std::unique_ptr<ASTNode>& n, States&&...)
- {
- for (size_t i = 1; i < n->children.size(); ++i) {
- n->children[0]->children.emplace_back(std::move(n->children[i]));
- }
- n = std::move(n->children[0]);
- }
-};
-
template <typename Rule>
using selector = TAO_PEGTL_NAMESPACE::parse_tree::selector<
Rule,
@@ -248,9 +234,6 @@ using selector = TAO_PEGTL_NAMESPACE::parse_tree::selector<
language::vector_type,
language::matrix_type,
language::string_type,
- language::cout_kw,
- language::cerr_kw,
- language::clog_kw,
language::var_declaration,
language::fct_declaration,
language::type_mapping,
@@ -269,6 +252,7 @@ using selector = TAO_PEGTL_NAMESPACE::parse_tree::selector<
language::equality,
language::compare,
language::sum,
+ language::shift,
language::product,
language::affectation,
language::expression>,
@@ -282,6 +266,8 @@ using selector = TAO_PEGTL_NAMESPACE::parse_tree::selector<
language::name_subscript_expression>,
TAO_PEGTL_NAMESPACE::parse_tree::remove_content::on<language::plus_op,
language::minus_op,
+ language::shift_left_op,
+ language::shift_right_op,
language::multiply_op,
language::divide_op,
language::lesser_op,
@@ -308,8 +294,7 @@ using selector = TAO_PEGTL_NAMESPACE::parse_tree::selector<
ASTBuilder::simplify_statement_block::on<language::statement_block>,
ASTBuilder::simplify_for_init::on<language::for_init>,
ASTBuilder::simplify_for_test::on<language::for_test>,
- ASTBuilder::simplify_for_post::on<language::for_post>,
- ASTBuilder::simplify_stream_statement::on<language::ostream_statement>>;
+ ASTBuilder::simplify_for_post::on<language::for_post>>;
template <typename InputT>
std::unique_ptr<ASTNode>
diff --git a/src/language/ast/ASTNodeBinaryOperatorExpressionBuilder.cpp b/src/language/ast/ASTNodeBinaryOperatorExpressionBuilder.cpp
index b6904220a9c42f2fcbfd49e8ce5881fbd9b7cf3a..c96b8799e3edbd037f70f8e68453bf700ced0fa2 100644
--- a/src/language/ast/ASTNodeBinaryOperatorExpressionBuilder.cpp
+++ b/src/language/ast/ASTNodeBinaryOperatorExpressionBuilder.cpp
@@ -22,6 +22,11 @@ ASTNodeBinaryOperatorExpressionBuilder::ASTNodeBinaryOperatorExpressionBuilder(A
} else if (n.is_type<language::minus_op>()) {
return binaryOperatorMangler<language::minus_op>(lhs_data_type, rhs_data_type);
+ } else if (n.is_type<language::shift_left_op>()) {
+ return binaryOperatorMangler<language::shift_left_op>(lhs_data_type, rhs_data_type);
+ } else if (n.is_type<language::shift_right_op>()) {
+ return binaryOperatorMangler<language::shift_right_op>(lhs_data_type, rhs_data_type);
+
} else if (n.is_type<language::or_op>()) {
return binaryOperatorMangler<language::or_op>(lhs_data_type, rhs_data_type);
} else if (n.is_type<language::and_op>()) {
diff --git a/src/language/ast/ASTNodeDataTypeBuilder.cpp b/src/language/ast/ASTNodeDataTypeBuilder.cpp
index eb307e2c16893be52128166c973f4992e86f3cb4..cb9e479df5f2e1f64edb4d5cb732acebcb505af6 100644
--- a/src/language/ast/ASTNodeDataTypeBuilder.cpp
+++ b/src/language/ast/ASTNodeDataTypeBuilder.cpp
@@ -161,8 +161,6 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
} else if (n.is_type<language::literal>()) {
n.m_data_type = ASTNodeDataType::build<ASTNodeDataType::string_t>();
- } else if (n.is_type<language::cout_kw>() or n.is_type<language::cerr_kw>() or n.is_type<language::clog_kw>()) {
- n.m_data_type = ASTNodeDataType::build<ASTNodeDataType::void_t>();
} else if (n.is_type<language::var_declaration>()) {
auto& name_node = *(n.children[0]);
auto& type_node = *(n.children[1]);
@@ -460,6 +458,7 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
}
} else if (n.is_type<language::plus_op>() or n.is_type<language::minus_op>() or
n.is_type<language::multiply_op>() or n.is_type<language::divide_op>() or
+ n.is_type<language::shift_left_op>() or n.is_type<language::shift_right_op>() or
n.is_type<language::lesser_op>() or n.is_type<language::lesser_or_eq_op>() or
n.is_type<language::greater_op>() or n.is_type<language::greater_or_eq_op>() or
n.is_type<language::eqeq_op>() or n.is_type<language::not_eq_op>() or n.is_type<language::and_op>() or
@@ -482,6 +481,12 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
return operator_repository.getBinaryOperatorValueType(
binaryOperatorMangler<language::divide_op>(type_0, type_1));
+ } else if (n.is_type<language::shift_left_op>()) {
+ return operator_repository.getBinaryOperatorValueType(
+ binaryOperatorMangler<language::shift_left_op>(type_0, type_1));
+ } else if (n.is_type<language::shift_right_op>()) {
+ return operator_repository.getBinaryOperatorValueType(
+ binaryOperatorMangler<language::shift_right_op>(type_0, type_1));
} else if (n.is_type<language::lesser_op>()) {
return operator_repository.getBinaryOperatorValueType(
binaryOperatorMangler<language::lesser_op>(type_0, type_1));
diff --git a/src/language/ast/ASTNodeExpressionBuilder.cpp b/src/language/ast/ASTNodeExpressionBuilder.cpp
index d9b3b4e439183e1f3a205a7c94e47252ac7a42f2..bbbcda43e788356d6bf4ec338643a79049069fc0 100644
--- a/src/language/ast/ASTNodeExpressionBuilder.cpp
+++ b/src/language/ast/ASTNodeExpressionBuilder.cpp
@@ -18,7 +18,6 @@
#include <language/node_processor/IfProcessor.hpp>
#include <language/node_processor/LocalNameProcessor.hpp>
#include <language/node_processor/NameProcessor.hpp>
-#include <language/node_processor/OStreamProcessor.hpp>
#include <language/node_processor/TupleToTinyVectorProcessor.hpp>
#include <language/node_processor/TupleToVectorProcessor.hpp>
#include <language/node_processor/ValueProcessor.hpp>
@@ -105,17 +104,12 @@ ASTNodeExpressionBuilder::_buildExpression(ASTNode& n)
} else if (n.is_type<language::multiply_op>() or n.is_type<language::divide_op>() or n.is_type<language::plus_op>() or
n.is_type<language::minus_op>() or n.is_type<language::or_op>() or n.is_type<language::and_op>() or
n.is_type<language::xor_op>() or n.is_type<language::greater_op>() or
+ n.is_type<language::shift_left_op>() or n.is_type<language::shift_right_op>() or
n.is_type<language::greater_or_eq_op>() or n.is_type<language::lesser_op>() or
n.is_type<language::lesser_or_eq_op>() or n.is_type<language::eqeq_op>() or
n.is_type<language::not_eq_op>()) {
ASTNodeBinaryOperatorExpressionBuilder{n};
- } else if (n.is_type<language::cout_kw>()) {
- n.m_node_processor = std::make_unique<OStreamProcessor>(n, std::cout);
- } else if (n.is_type<language::cerr_kw>()) {
- n.m_node_processor = std::make_unique<OStreamProcessor>(n, std::cerr);
- } else if (n.is_type<language::clog_kw>()) {
- n.m_node_processor = std::make_unique<OStreamProcessor>(n, std::clog);
} else if (n.is_type<language::if_statement>()) {
n.m_node_processor = std::make_unique<IfProcessor>(n);
} else if (n.is_type<language::statement_block>()) {
diff --git a/src/language/modules/BuiltinModule.cpp b/src/language/modules/BuiltinModule.cpp
index 98d8f71c6c91bd83ac53ff8d45b98f419c283fe9..cc76f78f7369f057c476a1ae15b953798d88903b 100644
--- a/src/language/modules/BuiltinModule.cpp
+++ b/src/language/modules/BuiltinModule.cpp
@@ -2,6 +2,7 @@
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/TypeDescriptor.hpp>
+#include <language/utils/ValueDescriptor.hpp>
#include <utils/Exceptions.hpp>
#include <memory>
@@ -44,6 +45,17 @@ BuiltinModule::_addBuiltinFunction(const std::string& name,
}
}
+void
+BuiltinModule::_addNameValue(const std::string& name, const ASTNodeDataType& type, const DataVariant& data)
+{
+ std::shared_ptr value_descriptor = std::make_shared<ValueDescriptor>(type, data);
+
+ auto [i_type, success] = m_name_value_map.insert(std::make_pair(name, value_descriptor));
+ if (not success) {
+ throw NormalError("variable '" + name + "' cannot be added!\n");
+ }
+}
+
void
BuiltinModule::_addTypeDescriptor(const ASTNodeDataType& ast_node_data_type)
{
diff --git a/src/language/modules/BuiltinModule.hpp b/src/language/modules/BuiltinModule.hpp
index 2be1d162578fb0d111e6eac7ae93e10dae981bce..3dd9ad1dbb9895f3f79ba04ec2702cc7ee0bd15f 100644
--- a/src/language/modules/BuiltinModule.hpp
+++ b/src/language/modules/BuiltinModule.hpp
@@ -6,18 +6,22 @@
class IBuiltinFunctionEmbedder;
class TypeDescriptor;
+class ValueDescriptor;
class BuiltinModule : public IModule
{
protected:
NameBuiltinFunctionMap m_name_builtin_function_map;
NameTypeMap m_name_type_map;
+ NameValueMap m_name_value_map;
void _addBuiltinFunction(const std::string& name,
std::shared_ptr<IBuiltinFunctionEmbedder> builtin_function_embedder);
void _addTypeDescriptor(const ASTNodeDataType& type);
+ void _addNameValue(const std::string& name, const ASTNodeDataType& type, const DataVariant& data);
+
const bool m_is_mandatory;
public:
@@ -39,6 +43,12 @@ class BuiltinModule : public IModule
return m_name_type_map;
}
+ const NameValueMap&
+ getNameValueMap() const final
+ {
+ return m_name_value_map;
+ }
+
BuiltinModule(bool is_mandatory = false);
~BuiltinModule() = default;
diff --git a/src/language/modules/CoreModule.cpp b/src/language/modules/CoreModule.cpp
index 1cc699455e23c123fbd1fee059f5296917a46b9b..9b70144e8559ccf97515bed9b71b174afb1a9311 100644
--- a/src/language/modules/CoreModule.cpp
+++ b/src/language/modules/CoreModule.cpp
@@ -22,6 +22,8 @@
#include <language/utils/IncDecOperatorRegisterForN.hpp>
#include <language/utils/IncDecOperatorRegisterForR.hpp>
#include <language/utils/IncDecOperatorRegisterForZ.hpp>
+#include <language/utils/OFStream.hpp>
+#include <language/utils/OStream.hpp>
#include <language/utils/UnaryOperatorRegisterForB.hpp>
#include <language/utils/UnaryOperatorRegisterForN.hpp>
#include <language/utils/UnaryOperatorRegisterForR.hpp>
@@ -83,6 +85,25 @@ CoreModule::CoreModule() : BuiltinModule(true)
[]() { RandomEngine::instance().resetRandomSeed(); }
));
+
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const OStream>>);
+
+ this
+ ->_addBuiltinFunction("ofstream",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const OStream>(const std::string&)>>(
+
+ [](const std::string& filename) { return std::make_shared<const OFStream>(filename); }
+
+ ));
+
+ this->_addNameValue("cout", ast_node_data_type_from<std::shared_ptr<const OStream>>,
+ EmbeddedData{std::make_shared<DataHandler<const OStream>>(std::make_shared<OStream>(std::cout))});
+
+ this->_addNameValue("cerr", ast_node_data_type_from<std::shared_ptr<const OStream>>,
+ EmbeddedData{std::make_shared<DataHandler<const OStream>>(std::make_shared<OStream>(std::cerr))});
+
+ this->_addNameValue("clog", ast_node_data_type_from<std::shared_ptr<const OStream>>,
+ EmbeddedData{std::make_shared<DataHandler<const OStream>>(std::make_shared<OStream>(std::clog))});
}
void
diff --git a/src/language/modules/IModule.hpp b/src/language/modules/IModule.hpp
index ceb3fd3dffc7d0af7f6423e3a21e7e8aa9686ffc..04aefbe91b57ab1925f55bae62ad233c1ea135cf 100644
--- a/src/language/modules/IModule.hpp
+++ b/src/language/modules/IModule.hpp
@@ -1,6 +1,8 @@
#ifndef IMODULE_HPP
#define IMODULE_HPP
+#include <language/utils/DataVariant.hpp>
+
#include <memory>
#include <string>
#include <string_view>
@@ -8,12 +10,14 @@
class IBuiltinFunctionEmbedder;
class TypeDescriptor;
+class ValueDescriptor;
class IModule
{
public:
using NameBuiltinFunctionMap = std::unordered_map<std::string, std::shared_ptr<IBuiltinFunctionEmbedder>>;
using NameTypeMap = std::unordered_map<std::string, std::shared_ptr<TypeDescriptor>>;
+ using NameValueMap = std::unordered_map<std::string, std::shared_ptr<ValueDescriptor>>;
IModule() = default;
IModule(IModule&&) = default;
@@ -25,6 +29,8 @@ class IModule
virtual const NameTypeMap& getNameTypeMap() const = 0;
+ virtual const NameValueMap& getNameValueMap() const = 0;
+
virtual void registerOperators() const = 0;
virtual std::string_view name() const = 0;
diff --git a/src/language/modules/ModuleRepository.cpp b/src/language/modules/ModuleRepository.cpp
index e516ab800833f5503de27bcfdd73e8b8381e06e6..71a1368b9dc374b19562bae27cc14a5d2b54c233 100644
--- a/src/language/modules/ModuleRepository.cpp
+++ b/src/language/modules/ModuleRepository.cpp
@@ -13,6 +13,8 @@
#include <language/utils/ParseError.hpp>
#include <language/utils/SymbolTable.hpp>
#include <language/utils/TypeDescriptor.hpp>
+#include <language/utils/ValueDescriptor.hpp>
+
#include <utils/PugsAssert.hpp>
#include <algorithm>
@@ -85,6 +87,28 @@ ModuleRepository::_populateEmbedderTableT(const ASTNode& module_node,
}
}
+void
+ModuleRepository::_populateSymbolTable(const ASTNode& module_node,
+ const std::string& module_name,
+ const IModule::NameValueMap& name_value_descriptor_map,
+ SymbolTable& symbol_table)
+{
+ for (auto [symbol_name, value_descriptor] : name_value_descriptor_map) {
+ auto [i_symbol, success] = symbol_table.add(symbol_name, module_node.begin());
+
+ if (not success) {
+ std::ostringstream error_message;
+ error_message << "importing module '" << module_name << "', cannot add symbol '" << symbol_name
+ << "', it is already defined!";
+ throw ParseError(error_message.str(), module_node.begin());
+ }
+
+ i_symbol->attributes().setDataType(value_descriptor->type());
+ i_symbol->attributes().setIsInitialized();
+ i_symbol->attributes().value() = value_descriptor->value();
+ }
+}
+
void
ModuleRepository::populateSymbolTable(const ASTNode& module_name_node, SymbolTable& symbol_table)
{
@@ -110,6 +134,8 @@ ModuleRepository::populateSymbolTable(const ASTNode& module_name_node, SymbolTab
ASTNodeDataType::build<ASTNodeDataType::type_name_id_t>(), symbol_table,
symbol_table.typeEmbedderTable());
+ this->_populateSymbolTable(module_name_node, module_name, populating_module.getNameValueMap(), symbol_table);
+
for (auto [symbol_name, embedded] : populating_module.getNameTypeMap()) {
BasicAffectationRegisterFor<EmbeddedData>(ASTNodeDataType::build<ASTNodeDataType::type_id_t>(symbol_name));
}
@@ -132,6 +158,12 @@ ModuleRepository::populateMandatorySymbolTable(const ASTNode& root_node, SymbolT
ASTNodeDataType::build<ASTNodeDataType::type_name_id_t>(), symbol_table,
symbol_table.typeEmbedderTable());
+ this->_populateSymbolTable(root_node, module_name, i_module->getNameValueMap(), symbol_table);
+
+ for (auto [symbol_name, embedded] : i_module->getNameTypeMap()) {
+ BasicAffectationRegisterFor<EmbeddedData>(ASTNodeDataType::build<ASTNodeDataType::type_id_t>(symbol_name));
+ }
+
i_module->registerOperators();
}
}
diff --git a/src/language/modules/ModuleRepository.hpp b/src/language/modules/ModuleRepository.hpp
index c4c9870fcedf550c82f94374f759f94d4c9cc878..c7289c8d1c4b9c7f795925bfdbf65e3f526082e7 100644
--- a/src/language/modules/ModuleRepository.hpp
+++ b/src/language/modules/ModuleRepository.hpp
@@ -26,6 +26,11 @@ class ModuleRepository
SymbolTable& symbol_table,
EmbedderTableT& embedder_table);
+ void _populateSymbolTable(const ASTNode& module_node,
+ const std::string& module_name,
+ const IModule::NameValueMap& name_value_map,
+ SymbolTable& symbol_table);
+
public:
void populateSymbolTable(const ASTNode& module_name_node, SymbolTable& symbol_table);
void populateMandatorySymbolTable(const ASTNode& root_node, SymbolTable& symbol_table);
diff --git a/src/language/node_processor/BinaryExpressionProcessor.hpp b/src/language/node_processor/BinaryExpressionProcessor.hpp
index 2c5d7460f26577633abed7281425d234b4d95432..7f278d4c7d4a007ad5db4c9f5d444665ccca56d6 100644
--- a/src/language/node_processor/BinaryExpressionProcessor.hpp
+++ b/src/language/node_processor/BinaryExpressionProcessor.hpp
@@ -82,6 +82,34 @@ struct BinOp<language::lesser_op>
}
};
+class OStream;
+
+template <>
+struct BinOp<language::shift_left_op>
+{
+ template <typename A, typename B>
+ PUGS_INLINE auto
+ eval(const A& a, const B& b) -> decltype(a << b)
+ {
+ if constexpr (std::is_same_v<A, std::shared_ptr<const OStream>> and std::is_same_v<B, bool>) {
+ return a << std::boolalpha << b;
+ } else {
+ return a << b;
+ }
+ }
+};
+
+template <>
+struct BinOp<language::shift_right_op>
+{
+ template <typename A, typename B>
+ PUGS_INLINE auto
+ eval(const A& a, const B& b) -> decltype(a >> b)
+ {
+ return a >> b;
+ }
+};
+
template <>
struct BinOp<language::lesser_or_eq_op>
{
@@ -286,7 +314,8 @@ struct BinaryExpressionProcessor<BinaryOpT, std::shared_ptr<ValueT>, std::shared
PUGS_INLINE DataVariant
_eval(const DataVariant& a, const DataVariant& b)
{
- if constexpr ((std::is_arithmetic_v<B_DataT>) or (is_tiny_matrix_v<B_DataT>) or (is_tiny_vector_v<B_DataT>)) {
+ if constexpr ((std::is_arithmetic_v<B_DataT>) or (is_tiny_matrix_v<B_DataT>) or (is_tiny_vector_v<B_DataT>) or
+ (std::is_same_v<std::string, B_DataT>)) {
const auto& embedded_a = std::get<EmbeddedData>(a);
const auto& b_value = std::get<B_DataT>(b);
diff --git a/src/language/node_processor/OStreamProcessor.hpp b/src/language/node_processor/OStreamProcessor.hpp
deleted file mode 100644
index e67dd7b6d9fcfa2440a59367efecd42ffd35a8ca..0000000000000000000000000000000000000000
--- a/src/language/node_processor/OStreamProcessor.hpp
+++ /dev/null
@@ -1,38 +0,0 @@
-#ifndef OSTREAM_PROCESSOR_HPP
-#define OSTREAM_PROCESSOR_HPP
-
-#include <language/ast/ASTNode.hpp>
-#include <language/node_processor/INodeProcessor.hpp>
-#include <language/utils/ParseError.hpp>
-
-class OStreamProcessor final : public INodeProcessor
-{
- private:
- ASTNode& m_node;
- std::ostream& m_os;
-
- public:
- DataVariant
- execute(ExecutionPolicy& exec_policy)
- {
- for (size_t i = 0; i < m_node.children.size(); ++i) {
- m_os << m_node.children[i]->execute(exec_policy);
- }
-
- return {};
- }
-
- OStreamProcessor(ASTNode& node, std::ostream& os) : m_node{node}, m_os(os)
- {
- for (auto& child : m_node.children) {
- if ((child->m_data_type == ASTNodeDataType::type_name_id_t) or
- (child->m_data_type == ASTNodeDataType::function_t) or
- (child->m_data_type == ASTNodeDataType::builtin_function_t)) {
- throw ParseError("invalid argument, cannot print a '" + dataTypeName(child->m_data_type) + "'",
- std::vector{child->begin()});
- }
- }
- }
-};
-
-#endif // OSTREAM_PROCESSOR_HPP
diff --git a/src/language/utils/BinaryOperatorMangler.hpp b/src/language/utils/BinaryOperatorMangler.hpp
index f3648f3fb860a6207b70ad712aee9935b669e493..d7f0c6c79358ffb555f20b1f55036cde2ce47504 100644
--- a/src/language/utils/BinaryOperatorMangler.hpp
+++ b/src/language/utils/BinaryOperatorMangler.hpp
@@ -13,6 +13,9 @@ struct divide_op;
struct plus_op;
struct minus_op;
+struct shift_left_op;
+struct shift_right_op;
+
struct or_op;
struct and_op;
struct xor_op;
@@ -56,6 +59,10 @@ binaryOperatorMangler(const ASTNodeDataType& lhs_data_type, const ASTNodeDataTyp
return "==";
} else if constexpr (std::is_same_v<BinaryOperatorT, language::not_eq_op>) {
return "!=";
+ } else if constexpr (std::is_same_v<BinaryOperatorT, language::shift_left_op>) {
+ return "<<";
+ } else if constexpr (std::is_same_v<BinaryOperatorT, language::shift_right_op>) {
+ return ">>";
} else {
static_assert(std::is_same_v<language::multiply_op, BinaryOperatorT>, "undefined binary operator");
}
diff --git a/src/language/utils/BinaryOperatorRegisterForB.cpp b/src/language/utils/BinaryOperatorRegisterForB.cpp
index cba803aaf3e3d7d4c72be7b69ca1a835203a49b2..dc228471d52330e486887e523c83dfa2d1917cef 100644
--- a/src/language/utils/BinaryOperatorRegisterForB.cpp
+++ b/src/language/utils/BinaryOperatorRegisterForB.cpp
@@ -1,6 +1,18 @@
#include <language/utils/BinaryOperatorRegisterForB.hpp>
#include <language/utils/BasicBinaryOperatorRegisterComparisonOf.hpp>
+#include <language/utils/DataHandler.hpp>
+#include <language/utils/OStream.hpp>
+
+void
+BinaryOperatorRegisterForB::_register_ostream()
+{
+ OperatorRepository& repository = OperatorRepository::instance();
+
+ repository.addBinaryOperator<language::shift_left_op>(
+ std::make_shared<BinaryOperatorProcessorBuilder<language::shift_left_op, std::shared_ptr<const OStream>,
+ std::shared_ptr<const OStream>, bool>>());
+}
void
BinaryOperatorRegisterForB::_register_comparisons()
@@ -61,6 +73,7 @@ BinaryOperatorRegisterForB::_register_divide()
BinaryOperatorRegisterForB::BinaryOperatorRegisterForB()
{
+ this->_register_ostream();
this->_register_comparisons();
this->_register_logical_operators();
this->_register_plus();
diff --git a/src/language/utils/BinaryOperatorRegisterForB.hpp b/src/language/utils/BinaryOperatorRegisterForB.hpp
index 651411bf7f2cd8e8386d26c071220d36cb79a691..545b49c5010e063654fea2b8570e4fd6e490f0e6 100644
--- a/src/language/utils/BinaryOperatorRegisterForB.hpp
+++ b/src/language/utils/BinaryOperatorRegisterForB.hpp
@@ -4,6 +4,7 @@
class BinaryOperatorRegisterForB
{
private:
+ void _register_ostream();
void _register_comparisons();
void _register_logical_operators();
void _register_plus();
diff --git a/src/language/utils/BinaryOperatorRegisterForN.cpp b/src/language/utils/BinaryOperatorRegisterForN.cpp
index 66d589405545924af9e2de137f7738c9b40b8430..e32e2b53a5fb70af8b84069790dd8435a362e051 100644
--- a/src/language/utils/BinaryOperatorRegisterForN.cpp
+++ b/src/language/utils/BinaryOperatorRegisterForN.cpp
@@ -1,6 +1,18 @@
#include <language/utils/BinaryOperatorRegisterForN.hpp>
#include <language/utils/BasicBinaryOperatorRegisterComparisonOf.hpp>
+#include <language/utils/DataHandler.hpp>
+#include <language/utils/OStream.hpp>
+
+void
+BinaryOperatorRegisterForN::_register_ostream()
+{
+ OperatorRepository& repository = OperatorRepository::instance();
+
+ repository.addBinaryOperator<language::shift_left_op>(
+ std::make_shared<BinaryOperatorProcessorBuilder<language::shift_left_op, std::shared_ptr<const OStream>,
+ std::shared_ptr<const OStream>, uint64_t>>());
+}
void
BinaryOperatorRegisterForN::_register_comparisons()
@@ -42,6 +54,7 @@ BinaryOperatorRegisterForN::_register_minus()
BinaryOperatorRegisterForN::BinaryOperatorRegisterForN()
{
+ this->_register_ostream();
this->_register_comparisons();
this->_register_arithmetic<language::plus_op>();
this->_register_minus();
diff --git a/src/language/utils/BinaryOperatorRegisterForN.hpp b/src/language/utils/BinaryOperatorRegisterForN.hpp
index 80dc52b013f36e05e08f7bf561f2ac3b41ccbfbf..997b6be7ab35284723ee64b6630ecb514f773bcd 100644
--- a/src/language/utils/BinaryOperatorRegisterForN.hpp
+++ b/src/language/utils/BinaryOperatorRegisterForN.hpp
@@ -4,6 +4,7 @@
class BinaryOperatorRegisterForN
{
private:
+ void _register_ostream();
template <typename OperatorT>
void _register_arithmetic();
void _register_comparisons();
diff --git a/src/language/utils/BinaryOperatorRegisterForR.cpp b/src/language/utils/BinaryOperatorRegisterForR.cpp
index 4f69988925e373fbd7369bff3d0bb87860b1aa3d..cc8ba8a79dc59c447317200a308103b8d7dfc1d8 100644
--- a/src/language/utils/BinaryOperatorRegisterForR.cpp
+++ b/src/language/utils/BinaryOperatorRegisterForR.cpp
@@ -1,6 +1,18 @@
#include <language/utils/BinaryOperatorRegisterForR.hpp>
#include <language/utils/BasicBinaryOperatorRegisterComparisonOf.hpp>
+#include <language/utils/DataHandler.hpp>
+#include <language/utils/OStream.hpp>
+
+void
+BinaryOperatorRegisterForR::_register_ostream()
+{
+ OperatorRepository& repository = OperatorRepository::instance();
+
+ repository.addBinaryOperator<language::shift_left_op>(
+ std::make_shared<BinaryOperatorProcessorBuilder<language::shift_left_op, std::shared_ptr<const OStream>,
+ std::shared_ptr<const OStream>, double>>());
+}
void
BinaryOperatorRegisterForR::_register_comparisons()
@@ -44,6 +56,7 @@ BinaryOperatorRegisterForR::_register_arithmetic()
BinaryOperatorRegisterForR::BinaryOperatorRegisterForR()
{
+ this->_register_ostream();
this->_register_comparisons();
this->_register_arithmetic<language::plus_op>();
this->_register_arithmetic<language::minus_op>();
diff --git a/src/language/utils/BinaryOperatorRegisterForR.hpp b/src/language/utils/BinaryOperatorRegisterForR.hpp
index cbde82fb3d36e565ef84d99b95b58c6c0556eb01..c14fa410f9a4dac520ae4482d3a62114f4840029 100644
--- a/src/language/utils/BinaryOperatorRegisterForR.hpp
+++ b/src/language/utils/BinaryOperatorRegisterForR.hpp
@@ -4,6 +4,7 @@
class BinaryOperatorRegisterForR
{
private:
+ void _register_ostream();
template <typename OperatorT>
void _register_arithmetic();
void _register_comparisons();
diff --git a/src/language/utils/BinaryOperatorRegisterForRn.cpp b/src/language/utils/BinaryOperatorRegisterForRn.cpp
index 41e5b031f7e3db02d8e59a9962e028c046f813a8..6e900d036871ed79eb13dfa737445becd50a1a9a 100644
--- a/src/language/utils/BinaryOperatorRegisterForRn.cpp
+++ b/src/language/utils/BinaryOperatorRegisterForRn.cpp
@@ -1,8 +1,23 @@
#include <language/utils/BinaryOperatorRegisterForRn.hpp>
#include <language/utils/BinaryOperatorProcessorBuilder.hpp>
+#include <language/utils/DataHandler.hpp>
+#include <language/utils/OStream.hpp>
#include <language/utils/OperatorRepository.hpp>
+template <size_t Dimension>
+void
+BinaryOperatorRegisterForRn<Dimension>::_register_ostream()
+{
+ OperatorRepository& repository = OperatorRepository::instance();
+
+ using Rn = TinyVector<Dimension>;
+
+ repository.addBinaryOperator<language::shift_left_op>(
+ std::make_shared<BinaryOperatorProcessorBuilder<language::shift_left_op, std::shared_ptr<const OStream>,
+ std::shared_ptr<const OStream>, Rn>>());
+}
+
template <size_t Dimension>
void
BinaryOperatorRegisterForRn<Dimension>::_register_comparisons()
@@ -54,6 +69,8 @@ BinaryOperatorRegisterForRn<Dimension>::_register_arithmetic()
template <size_t Dimension>
BinaryOperatorRegisterForRn<Dimension>::BinaryOperatorRegisterForRn()
{
+ this->_register_ostream();
+
this->_register_comparisons();
this->_register_product_by_a_scalar();
diff --git a/src/language/utils/BinaryOperatorRegisterForRn.hpp b/src/language/utils/BinaryOperatorRegisterForRn.hpp
index b9f27d4165d236b913e3393e342a144019fe6015..fc83d3a670a18b49824d7f59360c41dac9d39b6d 100644
--- a/src/language/utils/BinaryOperatorRegisterForRn.hpp
+++ b/src/language/utils/BinaryOperatorRegisterForRn.hpp
@@ -7,10 +7,9 @@ template <size_t Dimension>
class BinaryOperatorRegisterForRn
{
private:
+ void _register_ostream();
void _register_comparisons();
-
void _register_product_by_a_scalar();
-
template <typename OperatorT>
void _register_arithmetic();
diff --git a/src/language/utils/BinaryOperatorRegisterForRnxn.cpp b/src/language/utils/BinaryOperatorRegisterForRnxn.cpp
index 3d8850e79a4129c9a2098baa14b869f48bb8d1ce..f164203ddaca12e83f5ae5770dcb4848d366d3a4 100644
--- a/src/language/utils/BinaryOperatorRegisterForRnxn.cpp
+++ b/src/language/utils/BinaryOperatorRegisterForRnxn.cpp
@@ -1,8 +1,23 @@
#include <language/utils/BinaryOperatorRegisterForRnxn.hpp>
#include <language/utils/BinaryOperatorProcessorBuilder.hpp>
+#include <language/utils/DataHandler.hpp>
+#include <language/utils/OStream.hpp>
#include <language/utils/OperatorRepository.hpp>
+template <size_t Dimension>
+void
+BinaryOperatorRegisterForRnxn<Dimension>::_register_ostream()
+{
+ OperatorRepository& repository = OperatorRepository::instance();
+
+ using Rnxn = TinyMatrix<Dimension>;
+
+ repository.addBinaryOperator<language::shift_left_op>(
+ std::make_shared<BinaryOperatorProcessorBuilder<language::shift_left_op, std::shared_ptr<const OStream>,
+ std::shared_ptr<const OStream>, Rnxn>>());
+}
+
template <size_t Dimension>
void
BinaryOperatorRegisterForRnxn<Dimension>::_register_comparisons()
@@ -68,6 +83,8 @@ BinaryOperatorRegisterForRnxn<Dimension>::_register_arithmetic()
template <size_t Dimension>
BinaryOperatorRegisterForRnxn<Dimension>::BinaryOperatorRegisterForRnxn()
{
+ this->_register_ostream();
+
this->_register_comparisons();
this->_register_product_by_a_scalar();
diff --git a/src/language/utils/BinaryOperatorRegisterForRnxn.hpp b/src/language/utils/BinaryOperatorRegisterForRnxn.hpp
index 594740b629b0a7201d64139fd8e2ce4a12b38cd2..9edad436b51a5606fafc6aa980fd1578084f7e85 100644
--- a/src/language/utils/BinaryOperatorRegisterForRnxn.hpp
+++ b/src/language/utils/BinaryOperatorRegisterForRnxn.hpp
@@ -7,11 +7,10 @@ template <size_t Dimension>
class BinaryOperatorRegisterForRnxn
{
private:
+ void _register_ostream();
void _register_comparisons();
-
void _register_product_by_a_scalar();
void _register_product_by_a_vector();
-
template <typename OperatorT>
void _register_arithmetic();
diff --git a/src/language/utils/BinaryOperatorRegisterForString.cpp b/src/language/utils/BinaryOperatorRegisterForString.cpp
index 0e89cbe00502938d817405c3f3f3d6d991f0433e..999950cf1b1310bbe52614cc6823f37be53f82a4 100644
--- a/src/language/utils/BinaryOperatorRegisterForString.cpp
+++ b/src/language/utils/BinaryOperatorRegisterForString.cpp
@@ -2,8 +2,20 @@
#include <language/utils/BinaryOperatorProcessorBuilder.hpp>
#include <language/utils/ConcatExpressionProcessorBuilder.hpp>
+#include <language/utils/DataHandler.hpp>
+#include <language/utils/OStream.hpp>
#include <language/utils/OperatorRepository.hpp>
+void
+BinaryOperatorRegisterForString::_register_ostream()
+{
+ OperatorRepository& repository = OperatorRepository::instance();
+
+ repository.addBinaryOperator<language::shift_left_op>(
+ std::make_shared<BinaryOperatorProcessorBuilder<language::shift_left_op, std::shared_ptr<const OStream>,
+ std::shared_ptr<const OStream>, std::string>>());
+}
+
void
BinaryOperatorRegisterForString::_register_comparisons()
{
@@ -27,6 +39,8 @@ BinaryOperatorRegisterForString::_register_concat()
BinaryOperatorRegisterForString::BinaryOperatorRegisterForString()
{
+ this->_register_ostream();
+
this->_register_comparisons();
this->_register_concat<bool>();
diff --git a/src/language/utils/BinaryOperatorRegisterForString.hpp b/src/language/utils/BinaryOperatorRegisterForString.hpp
index 36af13ee6442086244ee5234ad496f3580f79c22..b7d4a74402092b60ae958c1beb20e61b4151592d 100644
--- a/src/language/utils/BinaryOperatorRegisterForString.hpp
+++ b/src/language/utils/BinaryOperatorRegisterForString.hpp
@@ -4,6 +4,8 @@
class BinaryOperatorRegisterForString
{
private:
+ void _register_ostream();
+
void _register_comparisons();
template <typename RHS_T>
diff --git a/src/language/utils/BinaryOperatorRegisterForZ.cpp b/src/language/utils/BinaryOperatorRegisterForZ.cpp
index cc8b22d117a8bebc557efff84083b8e45dcf8602..1d318850226beac44d27fe5f4fdba6c8e36a4e61 100644
--- a/src/language/utils/BinaryOperatorRegisterForZ.cpp
+++ b/src/language/utils/BinaryOperatorRegisterForZ.cpp
@@ -1,6 +1,18 @@
#include <language/utils/BinaryOperatorRegisterForZ.hpp>
#include <language/utils/BasicBinaryOperatorRegisterComparisonOf.hpp>
+#include <language/utils/DataHandler.hpp>
+#include <language/utils/OStream.hpp>
+
+void
+BinaryOperatorRegisterForZ::_register_ostream()
+{
+ OperatorRepository& repository = OperatorRepository::instance();
+
+ repository.addBinaryOperator<language::shift_left_op>(
+ std::make_shared<BinaryOperatorProcessorBuilder<language::shift_left_op, std::shared_ptr<const OStream>,
+ std::shared_ptr<const OStream>, int64_t>>());
+}
void
BinaryOperatorRegisterForZ::_register_comparisons()
@@ -36,6 +48,7 @@ BinaryOperatorRegisterForZ::_register_arithmetic()
BinaryOperatorRegisterForZ::BinaryOperatorRegisterForZ()
{
+ this->_register_ostream();
this->_register_comparisons();
this->_register_arithmetic<language::plus_op>();
this->_register_arithmetic<language::minus_op>();
diff --git a/src/language/utils/BinaryOperatorRegisterForZ.hpp b/src/language/utils/BinaryOperatorRegisterForZ.hpp
index 3f9a7c30c2f8e3674e6ebda7e9069153a6162f13..cb734843dc5013637ec597a8357a50d05a3b5373 100644
--- a/src/language/utils/BinaryOperatorRegisterForZ.hpp
+++ b/src/language/utils/BinaryOperatorRegisterForZ.hpp
@@ -4,6 +4,7 @@
class BinaryOperatorRegisterForZ
{
private:
+ void _register_ostream();
template <typename OperatorT>
void _register_arithmetic();
void _register_comparisons();
diff --git a/src/language/utils/CMakeLists.txt b/src/language/utils/CMakeLists.txt
index a86bd52d4a89b6102f262e6de7caae4355b99911..cd3f8af2314e96e87bacec25ac1dfc08d015192b 100644
--- a/src/language/utils/CMakeLists.txt
+++ b/src/language/utils/CMakeLists.txt
@@ -30,6 +30,7 @@ add_library(PugsLanguageUtils
IncDecOperatorRegisterForN.cpp
IncDecOperatorRegisterForR.cpp
IncDecOperatorRegisterForZ.cpp
+ OFStream.cpp
OperatorRepository.cpp
UnaryOperatorRegisterForB.cpp
UnaryOperatorRegisterForN.cpp
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index 52e69c29af17a994d1a51d8d87de48c02249645c..fa746269722229bff7fbd15f47edf926f14784cc 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -3,6 +3,7 @@
#include <language/utils/EmbeddedIDiscreteFunctionUtils.hpp>
#include <mesh/IMesh.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/IDiscreteFunction.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
@@ -318,45 +319,65 @@ template <size_t Dimension>
std::shared_ptr<const IDiscreteFunction>
dot(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
{
- Assert((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
- (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
- (f->dataType().dimension() == g->dataType().dimension()));
+ Assert(((f->descriptor().type() == DiscreteFunctionType::P0Vector) and
+ (g->descriptor().type() == DiscreteFunctionType::P0Vector)) or
+ ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+ (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
+ (f->dataType().dimension() == g->dataType().dimension())));
- using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+ if ((f->descriptor().type() == DiscreteFunctionType::P0Vector) and
+ (g->descriptor().type() == DiscreteFunctionType::P0Vector)) {
+ using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+ using DiscreteFunctionType = DiscreteFunctionP0Vector<Dimension, double>;
- switch (f->dataType().dimension()) {
- case 1: {
- using Rd = TinyVector<1>;
- using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
- return std::make_shared<const DiscreteFunctionResultType>(
- dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
- }
- case 2: {
- using Rd = TinyVector<2>;
- using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
- return std::make_shared<const DiscreteFunctionResultType>(
- dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
- }
- case 3: {
- using Rd = TinyVector<3>;
- using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
- return std::make_shared<const DiscreteFunctionResultType>(
- dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
- }
- // LCOV_EXCL_START
- default: {
- throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
- }
- // LCOV_EXCL_STOP
+ const DiscreteFunctionType& f_vector = dynamic_cast<const DiscreteFunctionType&>(*f);
+ const DiscreteFunctionType& g_vector = dynamic_cast<const DiscreteFunctionType&>(*g);
+
+ if (f_vector.size() != g_vector.size()) {
+ throw NormalError("operands have different dimension");
+ } else {
+ return std::make_shared<const DiscreteFunctionResultType>(dot(f_vector, g_vector));
+ }
+
+ } else {
+ using DiscreteFunctionResultType = DiscreteFunctionP0<Dimension, double>;
+
+ switch (f->dataType().dimension()) {
+ case 1: {
+ using Rd = TinyVector<1>;
+ using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
+ return std::make_shared<const DiscreteFunctionResultType>(
+ dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+ }
+ case 2: {
+ using Rd = TinyVector<2>;
+ using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
+ return std::make_shared<const DiscreteFunctionResultType>(
+ dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+ }
+ case 3: {
+ using Rd = TinyVector<3>;
+ using DiscreteFunctionType = DiscreteFunctionP0<Dimension, Rd>;
+ return std::make_shared<const DiscreteFunctionResultType>(
+ dot(dynamic_cast<const DiscreteFunctionType&>(*f), dynamic_cast<const DiscreteFunctionType&>(*g)));
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ // LCOV_EXCL_STOP
+ }
}
}
std::shared_ptr<const IDiscreteFunction>
dot(const std::shared_ptr<const IDiscreteFunction>& f, const std::shared_ptr<const IDiscreteFunction>& g)
{
- if ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
- (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
- (f->dataType().dimension() == g->dataType().dimension())) {
+ if (((f->descriptor().type() == DiscreteFunctionType::P0Vector) and
+ (g->descriptor().type() == DiscreteFunctionType::P0Vector)) or
+ ((f->dataType() == ASTNodeDataType::vector_t and f->descriptor().type() == DiscreteFunctionType::P0) and
+ (g->dataType() == ASTNodeDataType::vector_t and g->descriptor().type() == DiscreteFunctionType::P0) and
+ (f->dataType().dimension() == g->dataType().dimension()))) {
std::shared_ptr mesh = getCommonMesh({f, g});
if (mesh.use_count() == 0) {
diff --git a/src/language/utils/OFStream.cpp b/src/language/utils/OFStream.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..39b7cc8ecedee3ff76e3d6d54790d25f90d01569
--- /dev/null
+++ b/src/language/utils/OFStream.cpp
@@ -0,0 +1,17 @@
+#include <language/utils/OFStream.hpp>
+
+#include <utils/Messenger.hpp>
+
+OFStream::OFStream(const std::string& filename)
+{
+ if (parallel::rank() == 0) {
+ m_fstream.open(filename);
+ if (m_fstream.is_open()) {
+ m_ostream = &m_fstream;
+ } else {
+ std::stringstream error_msg;
+ error_msg << "cannot create file " << rang::fgB::yellow << filename << rang::style::reset;
+ throw NormalError(error_msg.str());
+ }
+ }
+}
diff --git a/src/language/utils/OFStream.hpp b/src/language/utils/OFStream.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8421906fca8746192b0d24730f7a731fe2640894
--- /dev/null
+++ b/src/language/utils/OFStream.hpp
@@ -0,0 +1,20 @@
+#ifndef OFSTREAM_HPP
+#define OFSTREAM_HPP
+
+#include <language/utils/OStream.hpp>
+
+#include <fstream>
+
+class OFStream final : public OStream
+{
+ private:
+ std::ofstream m_fstream;
+
+ public:
+ OFStream(const std::string& filename);
+
+ OFStream() = delete;
+ ~OFStream() = default;
+};
+
+#endif // OFSTREAM_HPP
diff --git a/src/language/utils/OStream.hpp b/src/language/utils/OStream.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d9410e25e35d2a3e9d5ac147911f6464d3b71eb7
--- /dev/null
+++ b/src/language/utils/OStream.hpp
@@ -0,0 +1,36 @@
+#ifndef OSTREAM_HPP
+#define OSTREAM_HPP
+
+#include <language/utils/ASTNodeDataTypeTraits.hpp>
+#include <utils/PugsAssert.hpp>
+
+#include <memory>
+
+class OStream
+{
+ protected:
+ mutable std::ostream* m_ostream = nullptr;
+
+ public:
+ template <typename DataT>
+ friend std::shared_ptr<const OStream>
+ operator<<(const std::shared_ptr<const OStream>& os, const DataT& t)
+ {
+ Assert(os.use_count() > 0, "non allocated stream");
+ if (os->m_ostream != nullptr) {
+ *os->m_ostream << t;
+ }
+ return os;
+ }
+
+ OStream(std::ostream& os) : m_ostream(&os) {}
+ OStream() = default;
+
+ virtual ~OStream() = default;
+};
+
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const OStream>> =
+ ASTNodeDataType::build<ASTNodeDataType::type_id_t>("ostream");
+
+#endif // OSTREAM_HPP
diff --git a/src/language/utils/ValueDescriptor.hpp b/src/language/utils/ValueDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..36bca35356b0b4a9759aee028b1b1ac02d9f7120
--- /dev/null
+++ b/src/language/utils/ValueDescriptor.hpp
@@ -0,0 +1,36 @@
+#ifndef VALUE_DESCRIPTOR_HPP
+#define VALUE_DESCRIPTOR_HPP
+
+#include <language/utils/ASTNodeDataType.hpp>
+#include <language/utils/DataVariant.hpp>
+
+#include <string>
+
+class ValueDescriptor
+{
+ private:
+ const ASTNodeDataType m_type;
+ const DataVariant m_value;
+
+ public:
+ const ASTNodeDataType
+ type() const
+ {
+ return m_type;
+ }
+
+ const DataVariant
+ value() const
+ {
+ return m_value;
+ }
+
+ ValueDescriptor(const ASTNodeDataType& type, const DataVariant& value) : m_type{type}, m_value{value} {}
+
+ ValueDescriptor(const ValueDescriptor&) = delete;
+ ValueDescriptor(ValueDescriptor&&) = delete;
+
+ ~ValueDescriptor() = default;
+};
+
+#endif // VALUE_DESCRIPTOR_HPP
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index b312e4a8e3296fceb7ca0b747bdf5841200b8274..4a3ab6939b02dffef76dbdda5e99273c6653bce9 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -85,7 +85,7 @@ class MeshData : public IMeshData
FaceValue<double> ll{m_mesh.connectivity()};
parallel_for(
- m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { ll[face_id] = L2Norm(Nl[face_id]); });
+ m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { ll[face_id] = l2Norm(Nl[face_id]); });
m_ll = ll;
}
@@ -103,11 +103,10 @@ class MeshData : public IMeshData
FaceValue<Rd> nl{m_mesh.connectivity()};
- const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
parallel_for(
m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { nl[face_id] = (1. / ll[face_id]) * Nl[face_id]; });
- m_Nl = Nl;
+ m_nl = nl;
}
}
diff --git a/src/output/OutputNamedItemValueSet.hpp b/src/output/OutputNamedItemValueSet.hpp
index bf90ddc3bad9636f2339695e87a446b11d81ae1a..69d6424d5bf3d4e9be9bbfb4ed3cee3281b3a786 100644
--- a/src/output/OutputNamedItemValueSet.hpp
+++ b/src/output/OutputNamedItemValueSet.hpp
@@ -67,9 +67,9 @@ class OutputNamedItemDataSet
NodeValue<const TinyVector<1, double>>,
NodeValue<const TinyVector<2, double>>,
NodeValue<const TinyVector<3, double>>,
- NodeValue<const TinyMatrix<1, double>>,
- NodeValue<const TinyMatrix<2, double>>,
- NodeValue<const TinyMatrix<3, double>>,
+ NodeValue<const TinyMatrix<1, 1, double>>,
+ NodeValue<const TinyMatrix<2, 2, double>>,
+ NodeValue<const TinyMatrix<3, 3, double>>,
CellValue<const bool>,
CellValue<const int>,
@@ -79,9 +79,9 @@ class OutputNamedItemDataSet
CellValue<const TinyVector<1, double>>,
CellValue<const TinyVector<2, double>>,
CellValue<const TinyVector<3, double>>,
- CellValue<const TinyMatrix<1, double>>,
- CellValue<const TinyMatrix<2, double>>,
- CellValue<const TinyMatrix<3, double>>,
+ CellValue<const TinyMatrix<1, 1, double>>,
+ CellValue<const TinyMatrix<2, 2, double>>,
+ CellValue<const TinyMatrix<3, 3, double>>,
NodeArray<const bool>,
NodeArray<const int>,
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 0bbdde090daab2cb49a086d359dbbe3b4e53739a..d01476cac3da67c1400efef7e3603e9bb91e096d 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -273,7 +273,7 @@ template <size_t N, typename DataType>
void
VTKWriter::_write_node_value(std::ofstream& os,
const std::string& name,
- const NodeValue<const TinyMatrix<N, DataType>>& item_value,
+ const NodeValue<const TinyMatrix<N, N, DataType>>& item_value,
SerializedDataList& serialized_data_list) const
{
os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N * N
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 7baede4556fb37e2b3b29f449dfbdaf5dfa8b143..913a05ef1de358cfe32f9b6faf1a9e6720585238 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -61,7 +61,7 @@ class VTKWriter : public WriterBase
template <size_t N, typename DataType>
void _write_node_value(std::ofstream& os,
const std::string& name,
- const NodeValue<const TinyMatrix<N, DataType>>& item_value,
+ const NodeValue<const TinyMatrix<N, N, DataType>>& item_value,
SerializedDataList& serialized_data_list) const;
template <typename DataType>
diff --git a/src/output/WriterBase.cpp b/src/output/WriterBase.cpp
index 21589e732ecde0a85491438ae7dba259b8ed3060..adf3f4a7f38216558602b5e5634c89364ef3cf6e 100644
--- a/src/output/WriterBase.cpp
+++ b/src/output/WriterBase.cpp
@@ -89,19 +89,19 @@ WriterBase::_register(const std::string& name,
switch (data_type.numberOfRows()) {
case 1: {
_register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<1, double>>&>(i_discrete_function),
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<1, 1, double>>&>(i_discrete_function),
named_item_data_set);
break;
}
case 2: {
_register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<2, double>>&>(i_discrete_function),
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<2, 2, double>>&>(i_discrete_function),
named_item_data_set);
break;
}
case 3: {
_register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<3, double>>&>(i_discrete_function),
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<3, 3, double>>&>(i_discrete_function),
named_item_data_set);
break;
}
diff --git a/src/scheme/DiscreteFunctionP0Vector.hpp b/src/scheme/DiscreteFunctionP0Vector.hpp
index 103718b125e123692241715e697ab89869d078eb..5f00d5964d3b0e060549a6337b62e40fe2b97f21 100644
--- a/src/scheme/DiscreteFunctionP0Vector.hpp
+++ b/src/scheme/DiscreteFunctionP0Vector.hpp
@@ -3,6 +3,7 @@
#include <scheme/IDiscreteFunction.hpp>
+#include <algebra/Vector.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/ItemArray.hpp>
#include <mesh/Mesh.hpp>
@@ -194,6 +195,19 @@ class DiscreteFunctionP0Vector : public IDiscreteFunction
return product;
}
+ PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+ dot(const DiscreteFunctionP0Vector& f, const DiscreteFunctionP0Vector& g)
+ {
+ Assert(f.mesh() == g.mesh(), "functions are nor defined on the same mesh");
+ Assert(f.size() == g.size());
+ DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
+ parallel_for(
+ f.m_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(Vector{f[cell_id]}, Vector{g[cell_id]}); });
+
+ return result;
+ }
+
DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& mesh, size_t size)
: m_mesh{mesh}, m_cell_arrays{mesh->connectivity(), size}
{}
diff --git a/src/utils/BacktraceManager.cpp b/src/utils/BacktraceManager.cpp
index 26219cbd7b0eb3eb6a451e3ba0587af14445676a..0ba95d6a6edd0e52b059e0b28638384eec13b602 100644
--- a/src/utils/BacktraceManager.cpp
+++ b/src/utils/BacktraceManager.cpp
@@ -8,6 +8,14 @@
#include <rang.hpp>
#include <regex>
+bool BacktraceManager::s_show = false;
+
+void
+BacktraceManager::setShow(bool show_backtrace)
+{
+ s_show = show_backtrace;
+}
+
BacktraceManager::BacktraceManager()
{
const int size = 100;
@@ -26,28 +34,33 @@ BacktraceManager::BacktraceManager()
std::ostream&
operator<<(std::ostream& os, const BacktraceManager& btm)
{
- const std::vector<std::string>& lines = btm.m_lines;
-
- const std::regex mangled_function(R"%(\(.*\+)%");
- const int width = std::log10(lines.size()) + 1;
-
- for (size_t i_line = 0; i_line < lines.size(); ++i_line) {
- const auto& line = lines[i_line];
- os << rang::fg::green << "[" << std::setw(width) << i_line + 1 << '/' << lines.size() << "] " << rang::fg::reset;
- std::smatch matchex;
- if (std::regex_search(line, matchex, mangled_function)) {
- std::string prefix = matchex.prefix().str();
- std::string function = line.substr(matchex.position() + 1, matchex.length() - 2);
- std::string suffix = matchex.suffix().str();
-
- os << prefix << '(';
- if (function.size() > 0) {
- os << rang::style::bold << demangle(function) << rang::style::reset;
+ if (BacktraceManager::s_show) {
+ const std::vector<std::string>& lines = btm.m_lines;
+
+ const std::regex mangled_function(R"%(\(.*\+)%");
+ const int width = std::log10(lines.size()) + 1;
+
+ for (size_t i_line = 0; i_line < lines.size(); ++i_line) {
+ const auto& line = lines[i_line];
+ os << rang::fg::green << "[" << std::setw(width) << i_line + 1 << '/' << lines.size() << "] " << rang::fg::reset;
+ std::smatch matchex;
+ if (std::regex_search(line, matchex, mangled_function)) {
+ std::string prefix = matchex.prefix().str();
+ std::string function = line.substr(matchex.position() + 1, matchex.length() - 2);
+ std::string suffix = matchex.suffix().str();
+
+ os << prefix << '(';
+ if (function.size() > 0) {
+ os << rang::style::bold << demangle(function) << rang::style::reset;
+ }
+ os << '+' << suffix << '\n';
+ } else {
+ os << line << '\n';
}
- os << '+' << suffix << '\n';
- } else {
- os << line << '\n';
}
+ } else {
+ os << rang::fg::yellow << "\n[To display backtrace launch pugs with the --backtrace option]" << rang::style::reset
+ << '\n';
}
return os;
diff --git a/src/utils/BacktraceManager.hpp b/src/utils/BacktraceManager.hpp
index 0ff2d2d42575ee42498a33fa2fd1c235a3e20b22..e990ffd3b4060f71370a0f8b7386cdfcd6c03ea8 100644
--- a/src/utils/BacktraceManager.hpp
+++ b/src/utils/BacktraceManager.hpp
@@ -7,9 +7,12 @@
class BacktraceManager
{
private:
+ static bool s_show;
std::vector<std::string> m_lines;
public:
+ static void setShow(bool show_backtrace);
+
BacktraceManager();
friend std::ostream& operator<<(std::ostream& os, const BacktraceManager& btm);
diff --git a/src/utils/PETScWrapper.cpp b/src/utils/PETScWrapper.cpp
index f88bbfd8a6bfd9de92d153e54c68e111f0b692b6..18b16aa1435c2f27dc4d850c0727ddb4e259acb0 100644
--- a/src/utils/PETScWrapper.cpp
+++ b/src/utils/PETScWrapper.cpp
@@ -13,6 +13,7 @@ initialize([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[])
{
#ifdef PUGS_HAS_PETSC
PetscOptionsSetValue(NULL, "-no_signal_handler", "true");
+ PetscOptionsSetValue(NULL, "-fp_trap", "false");
PetscInitialize(&argc, &argv, 0, 0);
#endif // PUGS_HAS_PETSC
}
diff --git a/src/utils/PugsTraits.hpp b/src/utils/PugsTraits.hpp
index 504d694ad86a08811a4a87e2c8098dab77cd607d..46a2b1d6dc5923d8ee12668d88484875f404944c 100644
--- a/src/utils/PugsTraits.hpp
+++ b/src/utils/PugsTraits.hpp
@@ -12,7 +12,7 @@
template <size_t N, typename T>
class TinyVector;
-template <size_t N, typename T>
+template <size_t M, size_t N, typename T>
class TinyMatrix;
template <typename DataType, ItemType item_type, typename ConnectivityPtr>
@@ -30,10 +30,10 @@ inline constexpr bool is_trivially_castable<TinyVector<N, T>> = is_trivially_cas
template <size_t N, typename T>
inline constexpr bool is_trivially_castable<const TinyVector<N, T>> = is_trivially_castable<T>;
-template <size_t N, typename T>
-inline constexpr bool is_trivially_castable<TinyMatrix<N, T>> = is_trivially_castable<T>;
-template <size_t N, typename T>
-inline constexpr bool is_trivially_castable<const TinyMatrix<N, T>> = is_trivially_castable<T>;
+template <size_t M, size_t N, typename T>
+inline constexpr bool is_trivially_castable<TinyMatrix<M, N, T>> = is_trivially_castable<T>;
+template <size_t M, size_t N, typename T>
+inline constexpr bool is_trivially_castable<const TinyMatrix<M, N, T>> = is_trivially_castable<T>;
// Traits is_false
@@ -105,8 +105,8 @@ inline constexpr bool is_tiny_vector_v<TinyVector<N, T>> = true;
template <typename T>
inline constexpr bool is_tiny_matrix_v = false;
-template <size_t N, typename T>
-inline constexpr bool is_tiny_matrix_v<TinyMatrix<N, T>> = true;
+template <size_t M, size_t N, typename T>
+inline constexpr bool is_tiny_matrix_v<TinyMatrix<M, N, T>> = true;
// Trais is ItemValue
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index d35fa95f43cb03810814ec6b62ed039842de1289..c4372c47dd3b826fdff2be76704070a04b2f9fcf 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -1,5 +1,6 @@
#include <utils/PugsUtils.hpp>
+#include <utils/BacktraceManager.hpp>
#include <utils/BuildInfo.hpp>
#include <utils/ConsoleManager.hpp>
#include <utils/FPEManager.hpp>
@@ -79,6 +80,8 @@ std::string
initialize(int& argc, char* argv[])
{
parallel::Messenger::create(argc, argv);
+ bool enable_fpe = true;
+ bool enable_signals = true;
std::string filename;
{
@@ -100,10 +103,11 @@ initialize(int& argc, char* argv[])
bool enable_color = true;
app.add_flag("--color,!--no-color", enable_color, "Colorize console output [default: true]");
- bool enable_fpe = true;
app.add_flag("--fpe,!--no-fpe", enable_fpe, "Trap floating point exceptions [default: true]");
- bool enable_signals = true;
+ bool show_backtrace = false;
+ app.add_flag("-b,--backtrace,!--no-backtrace", show_backtrace, "Show backtrace on failure [default: false]");
+
app.add_flag("--signal,!--no-signal", enable_signals, "Catches signals [default: true]");
bool pause_on_error = false;
@@ -118,15 +122,17 @@ initialize(int& argc, char* argv[])
std::exit(app.exit(e, std::cout, std::cerr));
}
+ BacktraceManager::setShow(show_backtrace);
ConsoleManager::init(enable_color);
- FPEManager::init(enable_fpe);
SignalManager::setPauseForDebug(pause_on_error);
- SignalManager::init(enable_signals);
}
PETScWrapper::initialize(argc, argv);
SLEPcWrapper::initialize(argc, argv);
+ FPEManager::init(enable_fpe);
+ SignalManager::init(enable_signals);
+
setDefaultOMPEnvironment();
Kokkos::initialize(argc, argv);
std::cout << "----------------- " << rang::fg::green << "pugs exec info" << rang::fg::reset
diff --git a/src/utils/SmallArray.hpp b/src/utils/SmallArray.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b5186a7f6341606a92084fd7ee5521b5d252b4ee
--- /dev/null
+++ b/src/utils/SmallArray.hpp
@@ -0,0 +1,144 @@
+#ifndef SMALL_ARRAY_HPP
+#define SMALL_ARRAY_HPP
+
+#include <utils/InvalidData.hpp>
+#include <utils/NaNHelper.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+#include <utils/PugsUtils.hpp>
+
+#include <algorithm>
+
+template <typename DataType>
+class [[nodiscard]] SmallArray
+{
+ public:
+ using data_type = DataType;
+ using index_type = size_t;
+
+ private:
+ index_type m_size;
+ std::shared_ptr<DataType[]> m_values;
+
+ // Allows const version to access our data
+ friend SmallArray<std::add_const_t<DataType>>;
+
+ public:
+ PUGS_INLINE size_t size() const noexcept
+ {
+ return m_size;
+ }
+
+ friend PUGS_INLINE SmallArray<std::remove_const_t<DataType>> copy(const SmallArray<DataType>& source)
+ {
+ SmallArray<std::remove_const_t<DataType>> image(source.m_size);
+ std::copy(source.m_values.get(), source.m_values.get() + source.m_size, image.m_values.get());
+
+ return image;
+ }
+
+ friend PUGS_INLINE void copy_to(const SmallArray<DataType>& source,
+ const SmallArray<std::remove_const_t<DataType>>& destination)
+ {
+ Assert(source.size() == destination.size());
+ std::copy(source.m_values.get(), source.m_values.get() + source.m_size, destination.m_values.get());
+ }
+
+ PUGS_INLINE DataType& operator[](index_type i) const noexcept(NO_ASSERT)
+ {
+ Assert(i < m_size);
+ return m_values[i];
+ }
+
+ PUGS_INLINE
+ void fill(const DataType& data) const
+ {
+ static_assert(not std::is_const_v<DataType>, "Cannot modify SmallArray of const");
+ std::fill(m_values.get(), m_values.get() + m_size, data);
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE SmallArray& operator=(const SmallArray<DataType2>& array) noexcept
+ {
+ // ensures that DataType is the same as source DataType2
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot assign SmallArray of different type");
+ // ensures that const is not lost through copy
+ static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
+ "Cannot assign SmallArray of const to SmallArray of non-const");
+
+ m_size = array.size();
+ m_values = array.m_values;
+ return *this;
+ }
+
+ PUGS_INLINE
+ SmallArray& operator=(const SmallArray&) = default;
+
+ PUGS_INLINE
+ SmallArray& operator=(SmallArray&&) = default;
+
+ PUGS_INLINE
+ explicit SmallArray(size_t size) : m_size{size}, m_values{std::shared_ptr<DataType[]>(new DataType[size])}
+ {
+ static_assert(not std::is_const<DataType>(), "Cannot allocate SmallArray of const data: only view is "
+ "supported");
+
+#ifndef NDEBUG
+ if constexpr (not std::is_const_v<DataType>) {
+ using T = std::decay_t<DataType>;
+ if constexpr (std::is_arithmetic_v<T>) {
+ this->fill(invalid_data_v<T>);
+ } else if constexpr (is_tiny_vector_v<T>) {
+ this->fill(T{});
+ } else if constexpr (is_tiny_matrix_v<T>) {
+ this->fill(T{});
+ }
+ }
+#endif // NDEBUG
+ }
+
+ friend std::ostream& operator<<(std::ostream& os, const SmallArray& x)
+ {
+ if (x.size() > 0) {
+ os << 0 << ':' << NaNHelper(x[0]);
+ }
+ for (size_t i = 1; i < x.size(); ++i) {
+ os << ' ' << i << ':' << NaNHelper(x[i]);
+ }
+ return os;
+ }
+
+ PUGS_INLINE
+ SmallArray() = default;
+
+ PUGS_INLINE
+ SmallArray(const SmallArray&) = default;
+
+ template <typename DataType2>
+ PUGS_INLINE SmallArray(const SmallArray<DataType2>& array) noexcept
+ {
+ this->operator=(array);
+ }
+
+ PUGS_INLINE
+ SmallArray(SmallArray &&) = default;
+
+ PUGS_INLINE
+ ~SmallArray() = default;
+};
+
+// map, multimap, unordered_map and stack cannot be copied this way
+template <typename Container>
+PUGS_INLINE SmallArray<std::remove_const_t<typename Container::value_type>>
+convert_to_small_array(const Container& given_container)
+{
+ using DataType = typename Container::value_type;
+ SmallArray<std::remove_const_t<DataType>> array(given_container.size());
+ if (given_container.size() > 0) {
+ std::copy(begin(given_container), end(given_container), &(array[0]));
+ }
+ return array;
+}
+
+#endif // SMALL_ARRAY_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 63a95ae67d58d3028c9af1996955588f96992346..5902fa80c022ef4efea13be21c6865ad32204df1 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -59,7 +59,6 @@ add_executable (unit_tests
test_CRSMatrix.cpp
test_CRSMatrixDescriptor.cpp
test_DataVariant.cpp
- test_DenseMatrix.cpp
test_Demangle.cpp
test_DiscreteFunctionDescriptorP0.cpp
test_DiscreteFunctionDescriptorP0Vector.cpp
@@ -88,13 +87,14 @@ add_executable (unit_tests
test_MathModule.cpp
test_NameProcessor.cpp
test_NaNHelper.cpp
- test_OStreamProcessor.cpp
+ test_OStream.cpp
test_ParseError.cpp
test_PETScUtils.cpp
test_PugsAssert.cpp
test_PugsFunctionAdapter.cpp
test_PugsUtils.cpp
test_RevisionInfo.cpp
+ test_SmallArray.cpp
test_SymbolTable.cpp
test_Table.cpp
test_Timer.cpp
@@ -120,6 +120,7 @@ add_executable (mpi_unit_tests
test_ItemValue.cpp
test_ItemValueUtils.cpp
test_Messenger.cpp
+ test_OFStream.cpp
test_Partitioner.cpp
test_RandomEngine.cpp
test_SubItemValuePerItem.cpp
diff --git a/tests/test_ASTBuilder.cpp b/tests/test_ASTBuilder.cpp
index 387e84b75b359b1cd4f95dd594f798e56b02208c..3ce2855eff1b059c0ac0e4215b8b1410e45af165 100644
--- a/tests/test_ASTBuilder.cpp
+++ b/tests/test_ASTBuilder.cpp
@@ -600,7 +600,7 @@ for(let i:N, i=0; i<10; ++i) {
CHECK_AST(data, result);
}
- SECTION("ostream simplifications")
+ SECTION("ostream")
{
std::string_view data = R"(
cout << 1+2 << "\n";
@@ -610,16 +610,22 @@ clog << "log " << l << "\n";
std::string_view result = R"(
(root)
- +-(language::cout_kw)
- | +-(language::plus_op)
- | | +-(language::integer:1)
- | | `-(language::integer:2)
+ +-(language::shift_left_op)
+ | +-(language::shift_left_op)
+ | | +-(language::name:cout)
+ | | `-(language::plus_op)
+ | | +-(language::integer:1)
+ | | `-(language::integer:2)
| `-(language::literal:"\n")
- +-(language::cerr_kw)
+ +-(language::shift_left_op)
+ | +-(language::name:cerr)
| `-(language::literal:"error?\n")
- `-(language::clog_kw)
- +-(language::literal:"log ")
- +-(language::name:l)
+ `-(language::shift_left_op)
+ +-(language::shift_left_op)
+ | +-(language::shift_left_op)
+ | | +-(language::name:clog)
+ | | `-(language::literal:"log ")
+ | `-(language::name:l)
`-(language::literal:"\n")
)";
CHECK_AST(data, result);
diff --git a/tests/test_ASTNodeDataTypeBuilder.cpp b/tests/test_ASTNodeDataTypeBuilder.cpp
index 96af38a0b905014ea518653795ba95850e486055..2e3835024a270781d0dccb37edd914abfdd63a51 100644
--- a/tests/test_ASTNodeDataTypeBuilder.cpp
+++ b/tests/test_ASTNodeDataTypeBuilder.cpp
@@ -1413,11 +1413,14 @@ clog << "clog\n";
std::string_view result = R"(
(root:void)
- +-(language::cout_kw:void)
+ +-(language::shift_left_op:ostream)
+ | +-(language::name:cout:ostream)
| `-(language::literal:"cout\n":string)
- +-(language::cerr_kw:void)
+ +-(language::shift_left_op:ostream)
+ | +-(language::name:cerr:ostream)
| `-(language::literal:"cerr\n":string)
- `-(language::clog_kw:void)
+ `-(language::shift_left_op:ostream)
+ +-(language::name:clog:ostream)
`-(language::literal:"clog\n":string)
)";
@@ -1445,8 +1448,10 @@ for (let i : N, i=0; i<3; ++i){
| `-(language::integer:3:Z)
+-(language::unary_plusplus:N)
| `-(language::name:i:N)
- `-(language::cout_kw:void)
- +-(language::name:i:N)
+ `-(language::shift_left_op:ostream)
+ +-(language::shift_left_op:ostream)
+ | +-(language::name:cout:ostream)
+ | `-(language::name:i:N)
`-(language::literal:"\n":string)
)";
diff --git a/tests/test_ASTNodeExpressionBuilder.cpp b/tests/test_ASTNodeExpressionBuilder.cpp
index cb6fdbfef85884a013e112a978207729ba965835..25a3d3a93f21294bbd04ae238ccbfcca5f4738fb 100644
--- a/tests/test_ASTNodeExpressionBuilder.cpp
+++ b/tests/test_ASTNodeExpressionBuilder.cpp
@@ -684,7 +684,7 @@ cout;
std::string result = R"(
(root:ASTNodeListProcessor)
- `-(language::cout_kw:OStreamProcessor)
+ `-(language::name:cout:NameProcessor)
)";
CHECK_AST(data, result);
@@ -698,7 +698,7 @@ cerr;
std::string result = R"(
(root:ASTNodeListProcessor)
- `-(language::cerr_kw:OStreamProcessor)
+ `-(language::name:cerr:NameProcessor)
)";
CHECK_AST(data, result);
@@ -712,7 +712,7 @@ clog;
std::string result = R"(
(root:ASTNodeListProcessor)
- `-(language::clog_kw:OStreamProcessor)
+ `-(language::name:clog:NameProcessor)
)";
CHECK_AST(data, result);
diff --git a/tests/test_ASTNodeFunctionExpressionBuilder.cpp b/tests/test_ASTNodeFunctionExpressionBuilder.cpp
index b37dd10869bb5665c9292b22e363af09be4b9ead..4bf775a679490306efecc2617547cb424c6a256b 100644
--- a/tests/test_ASTNodeFunctionExpressionBuilder.cpp
+++ b/tests/test_ASTNodeFunctionExpressionBuilder.cpp
@@ -655,7 +655,7 @@ f(1);
std::string_view result = R"(
(root:ASTNodeListProcessor)
- `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<1ul, double>, ZeroType>)
+ `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<1ul, 1ul, double>, ZeroType>)
+-(language::name:f:NameProcessor)
`-(language::integer:1:ValueProcessor)
)";
@@ -672,7 +672,7 @@ f(1);
std::string_view result = R"(
(root:ASTNodeListProcessor)
- `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<2ul, double>, ZeroType>)
+ `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<2ul, 2ul, double>, ZeroType>)
+-(language::name:f:NameProcessor)
`-(language::integer:1:ValueProcessor)
)";
@@ -689,7 +689,7 @@ f(1);
std::string_view result = R"(
(root:ASTNodeListProcessor)
- `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<3ul, double>, ZeroType>)
+ `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<3ul, 3ul, double>, ZeroType>)
+-(language::name:f:NameProcessor)
`-(language::integer:1:ValueProcessor)
)";
diff --git a/tests/test_ASTNodeListAffectationExpressionBuilder.cpp b/tests/test_ASTNodeListAffectationExpressionBuilder.cpp
index 6eab14fe5162b64a1d883cd54d8b9a70e20a5d66..8d495053e0d0e09b67c103a4f387ecf6afc018bf 100644
--- a/tests/test_ASTNodeListAffectationExpressionBuilder.cpp
+++ b/tests/test_ASTNodeListAffectationExpressionBuilder.cpp
@@ -239,17 +239,17 @@ let (x1,x2,x3,x) : R^1x1*R^2x2*R^3x3*R,
std::string_view result = R"(
(root:ASTNodeListProcessor)
- +-(language::eq_op:AffectationProcessor<language::eq_op, TinyMatrix<1ul, double>, long>)
+ +-(language::eq_op:AffectationProcessor<language::eq_op, TinyMatrix<1ul, 1ul, double>, long>)
| +-(language::name:a:NameProcessor)
| `-(language::integer:0:ValueProcessor)
- +-(language::eq_op:AffectationToTinyMatrixFromListProcessor<language::eq_op, TinyMatrix<2ul, double> >)
+ +-(language::eq_op:AffectationToTinyMatrixFromListProcessor<language::eq_op, TinyMatrix<2ul, 2ul, double> >)
| +-(language::name:b:NameProcessor)
| `-(language::expression_list:ASTNodeExpressionListProcessor)
| +-(language::integer:1:ValueProcessor)
| +-(language::integer:2:ValueProcessor)
| +-(language::integer:3:ValueProcessor)
| `-(language::integer:4:ValueProcessor)
- +-(language::eq_op:AffectationToTinyMatrixFromListProcessor<language::eq_op, TinyMatrix<3ul, double> >)
+ +-(language::eq_op:AffectationToTinyMatrixFromListProcessor<language::eq_op, TinyMatrix<3ul, 3ul, double> >)
| +-(language::name:c:NameProcessor)
| `-(language::expression_list:ASTNodeExpressionListProcessor)
| +-(language::integer:9:ValueProcessor)
diff --git a/tests/test_ASTNodeTypeCleaner.cpp b/tests/test_ASTNodeTypeCleaner.cpp
index e1abf88d600d45ba9ae77d3c25b4904adf6992d4..d78d4ce17fba3046d75e3ea1108d3874a3636dec 100644
--- a/tests/test_ASTNodeTypeCleaner.cpp
+++ b/tests/test_ASTNodeTypeCleaner.cpp
@@ -42,9 +42,12 @@ cout << "two=" << 2 << "\n";
std::string_view result = R"(
(root)
- `-(language::cout_kw)
- +-(language::literal:"two=")
- +-(language::integer:2)
+ `-(language::shift_left_op)
+ +-(language::shift_left_op)
+ | +-(language::shift_left_op)
+ | | +-(language::name:cout)
+ | | `-(language::literal:"two=")
+ | `-(language::integer:2)
`-(language::literal:"\n")
)";
diff --git a/tests/test_AffectationToTupleProcessor.cpp b/tests/test_AffectationToTupleProcessor.cpp
index 650e60ca61bb6bc1d8b2e239c9e4c67713c3212a..d832ac2739cb6b8417e6dfd3818e11252639b7b4 100644
--- a/tests/test_AffectationToTupleProcessor.cpp
+++ b/tests/test_AffectationToTupleProcessor.cpp
@@ -87,7 +87,7 @@ let s :(R^1); s = 1.3;
const std::string A_string = []() -> std::string {
std::ostringstream os;
- os << TinyMatrix<2, double>{1, 2, 3, 4};
+ os << TinyMatrix<2>{1, 2, 3, 4};
return os.str();
}();
@@ -162,7 +162,7 @@ let t :(R^1); t = (x,2);
const std::string A_string = []() -> std::string {
std::ostringstream os;
- os << TinyMatrix<2, double>{1, 2, 3, 4};
+ os << TinyMatrix<2>{1, 2, 3, 4};
return os.str();
}();
@@ -216,7 +216,7 @@ let s :(string); s = x;
const std::string A_string = []() -> std::string {
std::ostringstream os;
- os << TinyMatrix<3, double>{1, 2, 3, 4, 5, 6, 7, 8, 9};
+ os << TinyMatrix<3>{1, 2, 3, 4, 5, 6, 7, 8, 9};
return os.str();
}();
diff --git a/tests/test_ArrayUtils.cpp b/tests/test_ArrayUtils.cpp
index 79df610d909873a84328b59403ec861269b4f8b3..eb976b40df67c640bbe1a376671706835ca4c2e0 100644
--- a/tests/test_ArrayUtils.cpp
+++ b/tests/test_ArrayUtils.cpp
@@ -64,7 +64,7 @@ TEST_CASE("ArrayUtils", "[utils]")
SECTION("TinyMatrix Sum")
{
- using N22 = TinyMatrix<2, int>;
+ using N22 = TinyMatrix<2, 2, int>;
Array<N22> b(10);
b[0] = {13, 2, 0, 1};
b[1] = {1, 3, 6, 3};
diff --git a/tests/test_CRSMatrix.cpp b/tests/test_CRSMatrix.cpp
index d0d9476ab0e47ff90874391c63f39a6321adbc3b..cf6d65b87b10c362e89935ce89dc7ac79f490cd4 100644
--- a/tests/test_CRSMatrix.cpp
+++ b/tests/test_CRSMatrix.cpp
@@ -165,6 +165,177 @@ TEST_CASE("CRSMatrix", "[algebra]")
REQUIRE(y[3] == 150);
}
+ SECTION("square matrices sum")
+ {
+ Array<int> non_zeros_DA{4};
+ non_zeros_DA.fill(2);
+ CRSMatrixDescriptor<int> DA{4, 4, non_zeros_DA};
+ DA(0, 0) = 1;
+ DA(0, 1) = 2;
+ DA(1, 0) = 5;
+ DA(1, 1) = 6;
+ DA(2, 0) = 9;
+ DA(2, 3) = 12;
+ DA(3, 0) = 13;
+ DA(3, 3) = 16;
+
+ CRSMatrix<int> A{DA.getCRSMatrix()};
+
+ Array<int> non_zeros_DB{4};
+ non_zeros_DB.fill(2);
+ CRSMatrixDescriptor<int> DB{4, 4, non_zeros_DB};
+ DB(0, 0) = 1;
+ DB(0, 2) = 2;
+ DB(1, 0) = 5;
+ DB(1, 2) = 6;
+ DB(2, 0) = 9;
+ DB(2, 2) = 12;
+ DB(3, 0) = 13;
+ DB(3, 3) = 16;
+
+ CRSMatrix<int> B{DB.getCRSMatrix()};
+
+ std::ostringstream ost;
+ ost << A + B;
+
+ std::string ref = R"(0| 0:2 1:2 2:2
+1| 0:10 1:6 2:6
+2| 0:18 2:12 3:12
+3| 0:26 3:32
+)";
+
+ REQUIRE(ost.str() == ref);
+ }
+
+ SECTION("square matrices difference")
+ {
+ Array<int> non_zeros_DA{4};
+ non_zeros_DA.fill(2);
+ CRSMatrixDescriptor<int> DA{4, 4, non_zeros_DA};
+ DA(0, 0) = 1;
+ DA(0, 1) = 2;
+ DA(1, 0) = 5;
+ DA(1, 1) = 6;
+ DA(2, 0) = 9;
+ DA(2, 3) = 12;
+ DA(3, 0) = 13;
+ DA(3, 3) = 16;
+
+ CRSMatrix<int> A{DA.getCRSMatrix()};
+
+ Array<int> non_zeros_DB{4};
+ non_zeros_DB.fill(2);
+ CRSMatrixDescriptor<int> DB{4, 4, non_zeros_DB};
+ DB(0, 0) = -1;
+ DB(0, 2) = 3;
+ DB(1, 0) = 7;
+ DB(1, 2) = 3;
+ DB(2, 0) = 7;
+ DB(2, 2) = 11;
+ DB(3, 0) = 3;
+ DB(3, 3) = 8;
+
+ CRSMatrix<int> B{DB.getCRSMatrix()};
+
+ std::ostringstream ost;
+ ost << A - B;
+
+ std::string ref = R"(0| 0:2 1:2 2:-3
+1| 0:-2 1:6 2:-3
+2| 0:2 2:-11 3:12
+3| 0:10 3:8
+)";
+
+ REQUIRE(ost.str() == ref);
+ }
+
+ SECTION("rectangle matrices sum")
+ {
+ Array<int> non_zeros_DA{4};
+ non_zeros_DA.fill(2);
+ CRSMatrixDescriptor<int> DA{4, 5, non_zeros_DA};
+ DA(0, 0) = 1;
+ DA(0, 1) = 2;
+ DA(0, 4) = 2;
+ DA(1, 0) = 5;
+ DA(1, 1) = 6;
+ DA(2, 0) = 9;
+ DA(2, 3) = 12;
+ DA(3, 0) = 13;
+ DA(3, 3) = 16;
+ DA(3, 4) = 16;
+
+ CRSMatrix<int> A{DA.getCRSMatrix()};
+
+ Array<int> non_zeros_DB{4};
+ non_zeros_DB.fill(2);
+ CRSMatrixDescriptor<int> DB{4, 5, non_zeros_DB};
+ DB(0, 0) = 1;
+ DB(0, 2) = 2;
+ DB(1, 0) = 5;
+ DB(1, 2) = 6;
+ DB(1, 4) = 3;
+ DB(2, 0) = 9;
+ DB(2, 2) = 12;
+ DB(3, 0) = 13;
+ DB(3, 3) = 16;
+
+ CRSMatrix<int> B{DB.getCRSMatrix()};
+
+ std::ostringstream ost;
+ ost << A + B;
+
+ std::string ref = R"(0| 0:2 1:2 2:2 4:2
+1| 0:10 1:6 2:6 4:3
+2| 0:18 2:12 3:12
+3| 0:26 3:32 4:16
+)";
+
+ REQUIRE(ost.str() == ref);
+ }
+
+ SECTION("rectangle matrices difference")
+ {
+ Array<int> non_zeros_DA{4};
+ non_zeros_DA.fill(2);
+ CRSMatrixDescriptor<int> DA{4, 3, non_zeros_DA};
+ DA(0, 0) = 1;
+ DA(0, 1) = 2;
+ DA(1, 0) = 5;
+ DA(1, 1) = 6;
+ DA(2, 0) = 9;
+ DA(2, 2) = 12;
+ DA(3, 0) = 13;
+ DA(3, 2) = 16;
+
+ CRSMatrix<int> A{DA.getCRSMatrix()};
+
+ Array<int> non_zeros_DB{4};
+ non_zeros_DB.fill(2);
+ CRSMatrixDescriptor<int> DB{4, 3, non_zeros_DB};
+ DB(0, 0) = -1;
+ DB(0, 2) = 3;
+ DB(1, 0) = 7;
+ DB(1, 2) = 3;
+ DB(2, 0) = 7;
+ DB(2, 2) = 11;
+ DB(3, 0) = 3;
+ DB(3, 1) = 8;
+
+ CRSMatrix<int> B{DB.getCRSMatrix()};
+
+ std::ostringstream ost;
+ ost << A - B;
+
+ std::string ref = R"(0| 0:2 1:2 2:-3
+1| 0:-2 1:6 2:-3
+2| 0:2 2:1
+3| 0:10 1:-8 2:16
+)";
+
+ REQUIRE(ost.str() == ref);
+ }
+
SECTION("check values")
{
Array<int> non_zeros{4};
@@ -237,14 +408,45 @@ TEST_CASE("CRSMatrix", "[algebra]")
}
#ifndef NDEBUG
- SECTION("incompatible runtime matrix/vector product")
+ SECTION("runtime incompatible matrix/vector product")
{
Array<int> non_zeros(2);
non_zeros.fill(0);
CRSMatrixDescriptor<int> S{2, 4, non_zeros};
CRSMatrix<int> A{S.getCRSMatrix()};
Vector<int> x{2};
- REQUIRE_THROWS_AS(A * x, AssertError);
+ REQUIRE_THROWS_WITH(A * x, "cannot compute product: incompatible matrix and vector sizes");
}
+
+ SECTION("runtime incompatible matrices sum")
+ {
+ Array<int> A_non_zeros(2);
+ A_non_zeros.fill(0);
+ CRSMatrixDescriptor<int> DA{2, 4, A_non_zeros};
+ CRSMatrix<int> A{DA.getCRSMatrix()};
+
+ Array<int> B_non_zeros(3);
+ B_non_zeros.fill(0);
+ CRSMatrixDescriptor<int> DB{3, 4, B_non_zeros};
+ CRSMatrix<int> B{DB.getCRSMatrix()};
+
+ REQUIRE_THROWS_WITH(A + B, "cannot apply inner binary operator on matrices of different sizes");
+ }
+
+ SECTION("runtime incompatible matrices difference")
+ {
+ Array<int> A_non_zeros(2);
+ A_non_zeros.fill(0);
+ CRSMatrixDescriptor<int> DA{2, 4, A_non_zeros};
+ CRSMatrix<int> A{DA.getCRSMatrix()};
+
+ Array<int> B_non_zeros(2);
+ B_non_zeros.fill(0);
+ CRSMatrixDescriptor<int> DB{2, 3, B_non_zeros};
+ CRSMatrix<int> B{DB.getCRSMatrix()};
+
+ REQUIRE_THROWS_WITH(A - B, "cannot apply inner binary operator on matrices of different sizes");
+ }
+
#endif // NDEBUG
}
diff --git a/tests/test_DenseMatrix.cpp b/tests/test_DenseMatrix.cpp
deleted file mode 100644
index b794689a09793f0bd844144c9b7ccdb39338f648..0000000000000000000000000000000000000000
--- a/tests/test_DenseMatrix.cpp
+++ /dev/null
@@ -1,470 +0,0 @@
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_all.hpp>
-
-#include <utils/PugsAssert.hpp>
-
-#include <algebra/DenseMatrix.hpp>
-#include <algebra/Vector.hpp>
-
-// Instantiate to ensure full coverage is performed
-template class DenseMatrix<int>;
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("DenseMatrix", "[algebra]")
-{
- SECTION("size")
- {
- DenseMatrix<int> A{2, 3};
- REQUIRE(A.numberOfRows() == 2);
- REQUIRE(A.numberOfColumns() == 3);
- }
-
- SECTION("is square")
- {
- REQUIRE(not DenseMatrix<int>{2, 3}.isSquare());
- REQUIRE(not DenseMatrix<int>{4, 3}.isSquare());
- REQUIRE(DenseMatrix<int>{1, 1}.isSquare());
- REQUIRE(DenseMatrix<int>{2, 2}.isSquare());
- REQUIRE(DenseMatrix<int>{3, 3}.isSquare());
- }
-
- SECTION("write access")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 0;
- A(0, 1) = 1;
- A(0, 2) = 2;
- A(1, 0) = 3;
- A(1, 1) = 4;
- A(1, 2) = 5;
- REQUIRE(A(0, 0) == 0);
- REQUIRE(A(0, 1) == 1);
- REQUIRE(A(0, 2) == 2);
- REQUIRE(A(1, 0) == 3);
- REQUIRE(A(1, 1) == 4);
- REQUIRE(A(1, 2) == 5);
- DenseMatrix<const int> const_A = A;
- REQUIRE(const_A(0, 0) == 0);
- REQUIRE(const_A(0, 1) == 1);
- REQUIRE(const_A(0, 2) == 2);
- REQUIRE(const_A(1, 0) == 3);
- REQUIRE(const_A(1, 1) == 4);
- REQUIRE(const_A(1, 2) == 5);
- }
-
- SECTION("fill")
- {
- DenseMatrix<int> A{2, 3};
- A.fill(2);
- REQUIRE(A(0, 0) == 2);
- REQUIRE(A(0, 1) == 2);
- REQUIRE(A(0, 2) == 2);
- REQUIRE(A(1, 0) == 2);
- REQUIRE(A(1, 1) == 2);
- REQUIRE(A(1, 2) == 2);
- }
-
- SECTION("zero matrix")
- {
- DenseMatrix<int> A{2, 3};
- A = zero;
- for (size_t i = 0; i < A.numberOfRows(); ++i) {
- for (size_t j = 0; j < A.numberOfColumns(); ++j) {
- REQUIRE(A(i, j) == 0);
- }
- }
- }
-
- SECTION("identity matrix")
- {
- DenseMatrix<int> A{3, 3};
- A = identity;
- for (size_t i = 0; i < A.numberOfRows(); ++i) {
- for (size_t j = 0; j < A.numberOfColumns(); ++j) {
- REQUIRE(A(i, j) == (i == j));
- }
- }
- }
-
- SECTION("copy constructor (shallow)")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 0;
- A(0, 1) = 1;
- A(0, 2) = 2;
- A(1, 0) = 3;
- A(1, 1) = 4;
- A(1, 2) = 5;
-
- const DenseMatrix<int> B = A;
- REQUIRE(B(0, 0) == 0);
- REQUIRE(B(0, 1) == 1);
- REQUIRE(B(0, 2) == 2);
- REQUIRE(B(1, 0) == 3);
- REQUIRE(B(1, 1) == 4);
- REQUIRE(B(1, 2) == 5);
- }
-
- SECTION("copy (deep)")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 0;
- A(0, 1) = 1;
- A(0, 2) = 2;
- A(1, 0) = 3;
- A(1, 1) = 4;
- A(1, 2) = 5;
-
- const DenseMatrix<int> B = copy(A);
-
- A(0, 0) = 10;
- A(0, 1) = 11;
- A(0, 2) = 12;
- A(1, 0) = 13;
- A(1, 1) = 14;
- A(1, 2) = 15;
-
- REQUIRE(B(0, 0) == 0);
- REQUIRE(B(0, 1) == 1);
- REQUIRE(B(0, 2) == 2);
- REQUIRE(B(1, 0) == 3);
- REQUIRE(B(1, 1) == 4);
- REQUIRE(B(1, 2) == 5);
-
- REQUIRE(A(0, 0) == 10);
- REQUIRE(A(0, 1) == 11);
- REQUIRE(A(0, 2) == 12);
- REQUIRE(A(1, 0) == 13);
- REQUIRE(A(1, 1) == 14);
- REQUIRE(A(1, 2) == 15);
- }
-
- SECTION("self scalar multiplication")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 0;
- A(0, 1) = 1;
- A(0, 2) = 2;
- A(1, 0) = 3;
- A(1, 1) = 4;
- A(1, 2) = 5;
-
- A *= 2;
- REQUIRE(A(0, 0) == 0);
- REQUIRE(A(0, 1) == 2);
- REQUIRE(A(0, 2) == 4);
- REQUIRE(A(1, 0) == 6);
- REQUIRE(A(1, 1) == 8);
- REQUIRE(A(1, 2) == 10);
- }
-
- SECTION("left scalar multiplication")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 0;
- A(0, 1) = 1;
- A(0, 2) = 2;
- A(1, 0) = 3;
- A(1, 1) = 4;
- A(1, 2) = 5;
-
- const DenseMatrix<int> B = 2 * A;
-
- REQUIRE(B(0, 0) == 0);
- REQUIRE(B(0, 1) == 2);
- REQUIRE(B(0, 2) == 4);
- REQUIRE(B(1, 0) == 6);
- REQUIRE(B(1, 1) == 8);
- REQUIRE(B(1, 2) == 10);
- }
-
- SECTION("product matrix vector")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 6;
- A(0, 1) = 1;
- A(0, 2) = 2;
- A(1, 0) = 3;
- A(1, 1) = 4;
- A(1, 2) = 5;
-
- Vector<int> x{3};
- x[0] = 7;
- x[1] = 3;
- x[2] = 4;
-
- Vector y = A * x;
- REQUIRE(y[0] == 53);
- REQUIRE(y[1] == 53);
- }
-
- SECTION("self scalar division")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 0;
- A(0, 1) = 2;
- A(0, 2) = 4;
- A(1, 0) = 6;
- A(1, 1) = 8;
- A(1, 2) = 10;
-
- A /= 2;
-
- REQUIRE(A(0, 0) == 0);
- REQUIRE(A(0, 1) == 1);
- REQUIRE(A(0, 2) == 2);
- REQUIRE(A(1, 0) == 3);
- REQUIRE(A(1, 1) == 4);
- REQUIRE(A(1, 2) == 5);
- }
-
- SECTION("self minus")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 0;
- A(0, 1) = 1;
- A(0, 2) = 2;
- A(1, 0) = 3;
- A(1, 1) = 4;
- A(1, 2) = 5;
-
- DenseMatrix<int> B{2, 3};
- B(0, 0) = 5;
- B(0, 1) = 6;
- B(0, 2) = 4;
- B(1, 0) = 2;
- B(1, 1) = 1;
- B(1, 2) = 3;
-
- A -= B;
-
- REQUIRE(A(0, 0) == -5);
- REQUIRE(A(0, 1) == -5);
- REQUIRE(A(0, 2) == -2);
- REQUIRE(A(1, 0) == 1);
- REQUIRE(A(1, 1) == 3);
- REQUIRE(A(1, 2) == 2);
- }
-
- SECTION("self sum")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 0;
- A(0, 1) = 1;
- A(0, 2) = 2;
- A(1, 0) = 3;
- A(1, 1) = 4;
- A(1, 2) = 5;
-
- DenseMatrix<int> B{2, 3};
- B(0, 0) = 5;
- B(0, 1) = 6;
- B(0, 2) = 4;
- B(1, 0) = 2;
- B(1, 1) = 1;
- B(1, 2) = 3;
-
- A += B;
-
- REQUIRE(A(0, 0) == 5);
- REQUIRE(A(0, 1) == 7);
- REQUIRE(A(0, 2) == 6);
- REQUIRE(A(1, 0) == 5);
- REQUIRE(A(1, 1) == 5);
- REQUIRE(A(1, 2) == 8);
- }
-
- SECTION("sum")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 6;
- A(0, 1) = 5;
- A(0, 2) = 4;
- A(1, 0) = 3;
- A(1, 1) = 2;
- A(1, 2) = 1;
-
- DenseMatrix<int> B{2, 3};
- B(0, 0) = 0;
- B(0, 1) = 1;
- B(0, 2) = 2;
- B(1, 0) = 3;
- B(1, 1) = 4;
- B(1, 2) = 5;
-
- DenseMatrix C = A + B;
- REQUIRE(C(0, 0) == 6);
- REQUIRE(C(0, 1) == 6);
- REQUIRE(C(0, 2) == 6);
- REQUIRE(C(1, 0) == 6);
- REQUIRE(C(1, 1) == 6);
- REQUIRE(C(1, 2) == 6);
- }
-
- SECTION("difference")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 6;
- A(0, 1) = 5;
- A(0, 2) = 4;
- A(1, 0) = 3;
- A(1, 1) = 2;
- A(1, 2) = 1;
-
- DenseMatrix<int> B{2, 3};
- B(0, 0) = 0;
- B(0, 1) = 1;
- B(0, 2) = 2;
- B(1, 0) = 3;
- B(1, 1) = 4;
- B(1, 2) = 5;
-
- DenseMatrix C = A - B;
- REQUIRE(C(0, 0) == 6);
- REQUIRE(C(0, 1) == 4);
- REQUIRE(C(0, 2) == 2);
- REQUIRE(C(1, 0) == 0);
- REQUIRE(C(1, 1) == -2);
- REQUIRE(C(1, 2) == -4);
- }
-
- SECTION("transpose")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 0;
- A(0, 1) = 1;
- A(0, 2) = 2;
- A(1, 0) = 3;
- A(1, 1) = 4;
- A(1, 2) = 5;
-
- DenseMatrix B = transpose(A);
- REQUIRE(B(0, 0) == 0);
- REQUIRE(B(0, 1) == 3);
- REQUIRE(B(1, 0) == 1);
- REQUIRE(B(1, 1) == 4);
- REQUIRE(B(2, 0) == 2);
- REQUIRE(B(2, 1) == 5);
- }
-
- SECTION("product matrix vector")
- {
- DenseMatrix<int> A{2, 3};
- A(0, 0) = 1;
- A(0, 1) = 2;
- A(0, 2) = 3;
- A(1, 0) = 4;
- A(1, 1) = 5;
- A(1, 2) = 6;
-
- DenseMatrix<int> B{3, 2};
- B(0, 0) = 2;
- B(0, 1) = 8;
- B(1, 0) = 4;
- B(1, 1) = 9;
- B(2, 0) = 6;
- B(2, 1) = 10;
-
- DenseMatrix C = A * B;
- REQUIRE(C(0, 0) == 28);
- REQUIRE(C(0, 1) == 56);
- REQUIRE(C(1, 0) == 64);
- REQUIRE(C(1, 1) == 137);
- }
-
- SECTION("output")
- {
- DenseMatrix<double> A(3, 2);
- A(0, 0) = 1;
- A(0, 1) = 3;
- A(1, 0) = -8;
- A(1, 1) = -2;
- A(2, 0) = 4;
- A(2, 1) = -5;
- std::ostringstream A_ost;
- A_ost << A;
-
- std::ostringstream ref_ost;
- ref_ost << "0| " << 0 << ':' << 1. << ' ' << 1 << ':' << 3. << '\n';
- ref_ost << "1| " << 0 << ':' << -8. << ' ' << 1 << ':' << -2. << '\n';
- ref_ost << "2| " << 0 << ':' << 4. << ' ' << 1 << ':' << -5. << '\n';
-
- REQUIRE(A_ost.str() == ref_ost.str());
- }
-
-#ifndef NDEBUG
- SECTION("output with signaling NaN")
- {
- DenseMatrix<double> A(3, 2);
- A(0, 0) = 1;
- A(1, 1) = -2;
- A(2, 1) = -5;
- std::ostringstream A_ost;
- A_ost << A;
-
- std::ostringstream ref_ost;
- ref_ost << "0| " << 0 << ':' << 1. << ' ' << 1 << ":nan\n";
- ref_ost << "1| " << 0 << ":nan " << 1 << ':' << -2. << '\n';
- ref_ost << "2| " << 0 << ":nan " << 1 << ':' << -5. << '\n';
-
- REQUIRE(A_ost.str() == ref_ost.str());
- }
-
- SECTION("non square identity matrix")
- {
- DenseMatrix<int> A(2, 3);
-
- REQUIRE_THROWS_WITH(A = identity, "identity must be a square matrix");
- }
-
- SECTION("invalid matrix vector product")
- {
- DenseMatrix<int> A(2, 3);
- Vector<int> u(2);
-
- REQUIRE_THROWS_WITH(A * u, "cannot compute matrix-vector product: incompatible sizes");
- }
-
- SECTION("invalid matrix product")
- {
- DenseMatrix<int> A(2, 3);
- DenseMatrix<int> B(2, 3);
-
- REQUIRE_THROWS_WITH(A * B, "cannot compute matrix product: incompatible sizes");
- }
-
- SECTION("invalid matrix substract")
- {
- DenseMatrix<int> A(2, 3);
- DenseMatrix<int> B(3, 2);
-
- REQUIRE_THROWS_WITH(A -= B, "cannot substract matrix: incompatible sizes");
- }
-
- SECTION("invalid matrix sum")
- {
- DenseMatrix<int> A(2, 3);
- DenseMatrix<int> B(3, 2);
-
- REQUIRE_THROWS_WITH(A + B, "cannot compute matrix sum: incompatible sizes");
- }
-
- SECTION("invalid matrix difference")
- {
- DenseMatrix<int> A(2, 3);
- DenseMatrix<int> B(3, 2);
-
- REQUIRE_THROWS_WITH(A - B, "cannot compute matrix difference: incompatible sizes");
- }
-
- SECTION("invalid indices")
- {
- DenseMatrix<int> A(2, 3);
-
- REQUIRE_THROWS_WITH(A(2, 0), "cannot access element: invalid indices");
- REQUIRE_THROWS_WITH(A(0, 3), "cannot access element: invalid indices");
- }
-
-#endif // NDEBUG
-}
diff --git a/tests/test_DiscreteFunctionP0.cpp b/tests/test_DiscreteFunctionP0.cpp
index acfa7a0064d5368c251af31ff3cb8b9a6a9a749a..f8d7bb6ef78d027287cf7d721299f051ba7c7a98 100644
--- a/tests/test_DiscreteFunctionP0.cpp
+++ b/tests/test_DiscreteFunctionP0.cpp
@@ -302,9 +302,9 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
SECTION("matrix")
{
- const TinyMatrix<3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
- const TinyMatrix<3, size_t> zero{ZeroType{}};
- DiscreteFunctionP0<Dimension, TinyMatrix<3, size_t>> f{mesh};
+ const TinyMatrix<3, 3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
+ const TinyMatrix<3, 3, size_t> zero{ZeroType{}};
+ DiscreteFunctionP0<Dimension, TinyMatrix<3, 3, size_t>> f{mesh};
f.fill(value);
REQUIRE(all_values_equal(f, value));
@@ -318,7 +318,7 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
copy_to(g, f);
g.fill(zero);
- DiscreteFunctionP0<Dimension, const TinyMatrix<3, size_t>> h = copy(f);
+ DiscreteFunctionP0<Dimension, const TinyMatrix<3, 3, size_t>> h = copy(f);
REQUIRE(all_values_equal(f, value));
REQUIRE(all_values_equal(g, zero));
@@ -397,9 +397,9 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
SECTION("matrix")
{
- const TinyMatrix<3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
- const TinyMatrix<3, size_t> zero{ZeroType{}};
- DiscreteFunctionP0<Dimension, TinyMatrix<3, size_t>> f{mesh};
+ const TinyMatrix<3, 3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
+ const TinyMatrix<3, 3, size_t> zero{ZeroType{}};
+ DiscreteFunctionP0<Dimension, TinyMatrix<3, 3, size_t>> f{mesh};
f.fill(value);
REQUIRE(all_values_equal(f, value));
@@ -413,7 +413,7 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
copy_to(g, f);
g.fill(zero);
- DiscreteFunctionP0<Dimension, const TinyMatrix<3, size_t>> h = copy(f);
+ DiscreteFunctionP0<Dimension, const TinyMatrix<3, 3, size_t>> h = copy(f);
REQUIRE(all_values_equal(f, value));
REQUIRE(all_values_equal(g, zero));
@@ -492,9 +492,9 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
SECTION("matrix")
{
- const TinyMatrix<3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
- const TinyMatrix<3, size_t> zero{ZeroType{}};
- DiscreteFunctionP0<Dimension, TinyMatrix<3, size_t>> f{mesh};
+ const TinyMatrix<3, 3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
+ const TinyMatrix<3, 3, size_t> zero{ZeroType{}};
+ DiscreteFunctionP0<Dimension, TinyMatrix<3, 3, size_t>> f{mesh};
f.fill(value);
REQUIRE(all_values_equal(f, value));
@@ -508,7 +508,7 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
copy_to(g, f);
g.fill(zero);
- DiscreteFunctionP0<Dimension, const TinyMatrix<3, size_t>> h = copy(f);
+ DiscreteFunctionP0<Dimension, const TinyMatrix<3, 3, size_t>> h = copy(f);
REQUIRE(all_values_equal(f, value));
REQUIRE(all_values_equal(g, zero));
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 3ca8f7f5136d72f7003357113b7b597d1b84ef54..fb7bf07907bf6c9e1db9333c9713365350918b06 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -3,8 +3,8 @@
#include <algebra/CRSMatrix.hpp>
#include <algebra/CRSMatrixDescriptor.hpp>
-#include <algebra/DenseMatrix.hpp>
#include <algebra/EigenvalueSolver.hpp>
+#include <algebra/SmallMatrix.hpp>
#include <utils/pugs_config.hpp>
@@ -36,7 +36,7 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
SECTION("eigenvalues")
{
- Array<double> eigenvalues;
+ SmallArray<double> eigenvalues;
#ifdef PUGS_HAS_SLEPC
EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalues);
@@ -51,8 +51,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
SECTION("eigenvalues and eigenvectors")
{
- Array<double> eigenvalue_list;
- std::vector<Vector<double>> eigenvector_list;
+ SmallArray<double> eigenvalue_list;
+ std::vector<SmallVector<double>> eigenvector_list;
#ifdef PUGS_HAS_SLEPC
EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
@@ -61,9 +61,9 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
- DenseMatrix<double> P{3};
- DenseMatrix<double> PT{3};
- DenseMatrix<double> D{3};
+ SmallMatrix<double> P{3};
+ SmallMatrix<double> PT{3};
+ SmallMatrix<double> D{3};
D = zero;
for (size_t i = 0; i < 3; ++i) {
@@ -74,7 +74,7 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
D(i, i) = eigenvalue_list[i];
}
- DenseMatrix PDPT = P * D * PT;
+ SmallMatrix PDPT = P * D * PT;
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
@@ -88,8 +88,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
SECTION("eigenvalues and passage matrix")
{
- Array<double> eigenvalue_list;
- DenseMatrix<double> P{};
+ SmallArray<double> eigenvalue_list;
+ SmallMatrix<double> P{};
#ifdef PUGS_HAS_SLEPC
EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P);
@@ -98,14 +98,14 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
- DenseMatrix<double> D{3};
+ SmallMatrix<double> D{3};
D = zero;
for (size_t i = 0; i < 3; ++i) {
D(i, i) = eigenvalue_list[i];
}
- DenseMatrix PT = transpose(P);
+ SmallMatrix PT = transpose(P);
- DenseMatrix PDPT = P * D * PT;
+ SmallMatrix PDPT = P * D * PT;
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
diff --git a/tests/test_EmbeddedIDiscreteFunctionUtils.cpp b/tests/test_EmbeddedIDiscreteFunctionUtils.cpp
index 72ce3374369ace3e93551c94094000d117cb336e..0ded74fd6b5764bfb67917a48a1a24ab8b1db873 100644
--- a/tests/test_EmbeddedIDiscreteFunctionUtils.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionUtils.cpp
@@ -15,9 +15,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionUtils", "[language]")
using R2 = TinyVector<2, double>;
using R3 = TinyVector<3, double>;
- using R1x1 = TinyMatrix<1, double>;
- using R2x2 = TinyMatrix<2, double>;
- using R3x3 = TinyMatrix<3, double>;
+ using R1x1 = TinyMatrix<1, 1, double>;
+ using R2x2 = TinyMatrix<2, 2, double>;
+ using R3x3 = TinyMatrix<3, 3, double>;
SECTION("operand type name")
{
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index f8bcd9972ca5b52a6ea5fe48d70eb1b8c2064107..11195faab3dad2f94e82dba78b9afc368d0edcd0 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -46,8 +46,8 @@ TEST_CASE("LinearSolver", "[algebra]")
options.method() = LSMethod::lu;
REQUIRE_THROWS_WITH(LinearSolver{options}, "error: LU is not a builtin linear solver!");
- options.method() = LSMethod::choleski;
- REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Choleski is not a builtin linear solver!");
+ options.method() = LSMethod::cholesky;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Cholesky is not a builtin linear solver!");
options.method() = LSMethod::gmres;
REQUIRE_THROWS_WITH(LinearSolver{options}, "error: GMRES is not a builtin linear solver!");
@@ -67,8 +67,8 @@ TEST_CASE("LinearSolver", "[algebra]")
options.precond() = LSPrecond::incomplete_LU;
REQUIRE_THROWS_WITH(LinearSolver{options}, "error: ILU is not a builtin preconditioner!");
- options.precond() = LSPrecond::incomplete_choleski;
- REQUIRE_THROWS_WITH(LinearSolver{options}, "error: ICholeski is not a builtin preconditioner!");
+ options.precond() = LSPrecond::incomplete_cholesky;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: ICholesky is not a builtin preconditioner!");
options.precond() = LSPrecond::amg;
REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not a builtin preconditioner!");
@@ -101,7 +101,7 @@ TEST_CASE("LinearSolver", "[algebra]")
options.method() = LSMethod::lu;
REQUIRE_NOTHROW(linear_solver.checkOptions(options));
- options.method() = LSMethod::choleski;
+ options.method() = LSMethod::cholesky;
REQUIRE_NOTHROW(linear_solver.checkOptions(options));
options.method() = LSMethod::gmres;
@@ -122,7 +122,7 @@ TEST_CASE("LinearSolver", "[algebra]")
options.precond() = LSPrecond::incomplete_LU;
REQUIRE_NOTHROW(linear_solver.checkOptions(options));
- options.precond() = LSPrecond::incomplete_choleski;
+ options.precond() = LSPrecond::incomplete_cholesky;
REQUIRE_NOTHROW(linear_solver.checkOptions(options));
options.precond() = LSPrecond::amg;
@@ -245,9 +245,9 @@ TEST_CASE("LinearSolver", "[algebra]")
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("CG ICholeski")
+ SECTION("CG ICholesky")
{
- options.precond() = LSPrecond::incomplete_choleski;
+ options.precond() = LSPrecond::incomplete_cholesky;
Vector<double> x{5};
x = zero;
@@ -274,11 +274,11 @@ TEST_CASE("LinearSolver", "[algebra]")
}
}
- SECTION("Choleski")
+ SECTION("Cholesky")
{
LinearSolverOptions options;
options.library() = LSLibrary::petsc;
- options.method() = LSMethod::choleski;
+ options.method() = LSMethod::cholesky;
options.precond() = LSPrecond::none;
Vector<double> x{5};
@@ -545,7 +545,7 @@ TEST_CASE("LinearSolver", "[algebra]")
{
SECTION("symmetric system")
{
- DenseMatrix<double> A{5};
+ SmallMatrix<double> A{5};
A = zero;
A(0, 0) = 2;
@@ -566,8 +566,8 @@ TEST_CASE("LinearSolver", "[algebra]")
A(4, 3) = -1;
A(4, 4) = 2;
- Vector<const double> x_exact = [] {
- Vector<double> y{5};
+ SmallVector<const double> x_exact = [] {
+ SmallVector<double> y{5};
y[0] = 1;
y[1] = 3;
y[2] = 2;
@@ -576,7 +576,7 @@ TEST_CASE("LinearSolver", "[algebra]")
return y;
}();
- Vector<double> b = A * x_exact;
+ SmallVector<double> b = A * x_exact;
SECTION("builtin")
{
@@ -587,13 +587,13 @@ TEST_CASE("LinearSolver", "[algebra]")
options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
}
@@ -614,13 +614,13 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::none;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
@@ -628,27 +628,27 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::diagonal;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("CG ICholeski")
+ SECTION("CG ICholesky")
{
- options.precond() = LSPrecond::incomplete_choleski;
+ options.precond() = LSPrecond::incomplete_cholesky;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
@@ -656,31 +656,31 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::amg;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
}
- SECTION("Choleski")
+ SECTION("Cholesky")
{
LinearSolverOptions options;
options.library() = LSLibrary::petsc;
- options.method() = LSMethod::choleski;
+ options.method() = LSMethod::cholesky;
options.precond() = LSPrecond::none;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
@@ -700,7 +700,7 @@ TEST_CASE("LinearSolver", "[algebra]")
SECTION("none symmetric system")
{
- DenseMatrix<double> A{5};
+ SmallMatrix<double> A{5};
A = zero;
A(0, 0) = 2;
@@ -721,8 +721,8 @@ TEST_CASE("LinearSolver", "[algebra]")
A(4, 3) = 1;
A(4, 4) = 3;
- Vector<const double> x_exact = [] {
- Vector<double> y{5};
+ SmallVector<const double> x_exact = [] {
+ SmallVector<double> y{5};
y[0] = 1;
y[1] = 3;
y[2] = 2;
@@ -731,7 +731,7 @@ TEST_CASE("LinearSolver", "[algebra]")
return y;
}();
- Vector<double> b = A * x_exact;
+ SmallVector<double> b = A * x_exact;
SECTION("builtin")
{
@@ -742,13 +742,13 @@ TEST_CASE("LinearSolver", "[algebra]")
options.method() = LSMethod::bicgstab;
options.precond() = LSPrecond::none;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
}
@@ -769,13 +769,13 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::none;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
@@ -783,13 +783,13 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::diagonal;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
@@ -797,13 +797,13 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::incomplete_LU;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
}
@@ -819,13 +819,13 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::none;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
@@ -833,13 +833,13 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::diagonal;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
}
@@ -855,13 +855,13 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::none;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
@@ -869,13 +869,13 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::diagonal;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
@@ -883,13 +883,13 @@ TEST_CASE("LinearSolver", "[algebra]")
{
options.precond() = LSPrecond::incomplete_LU;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
}
@@ -901,13 +901,13 @@ TEST_CASE("LinearSolver", "[algebra]")
options.method() = LSMethod::lu;
options.precond() = LSPrecond::none;
- Vector<double> x{5};
+ SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
+ SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
diff --git a/tests/test_LinearSolverOptions.cpp b/tests/test_LinearSolverOptions.cpp
index 3da7688892a736fe5f69c2adbab9c16cb6cd4fd9..90519ddf9797f0d6fa745cbc34e2f80443eed33d 100644
--- a/tests/test_LinearSolverOptions.cpp
+++ b/tests/test_LinearSolverOptions.cpp
@@ -17,7 +17,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]")
options.verbose() = true;
options.library() = LSLibrary::builtin;
options.method() = LSMethod::cg;
- options.precond() = LSPrecond::incomplete_choleski;
+ options.precond() = LSPrecond::incomplete_cholesky;
options.epsilon() = 1E-3;
options.maximumIteration() = 100;
@@ -28,7 +28,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]")
expected_output << R"(
library: builtin
method : CG
- precond: ICholeski
+ precond: ICholesky
epsilon: )" << 1E-3 << R"(
maxiter: 100
verbose: true
@@ -79,7 +79,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]")
REQUIRE(name(LSMethod::bicgstab2) == "BICGStab2");
REQUIRE(name(LSMethod::gmres) == "GMRES");
REQUIRE(name(LSMethod::lu) == "LU");
- REQUIRE(name(LSMethod::choleski) == "Choleski");
+ REQUIRE(name(LSMethod::cholesky) == "Cholesky");
REQUIRE_THROWS_WITH(name(LSMethod::LS__end), "unexpected error: Linear system method name is not defined!");
}
@@ -87,7 +87,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]")
{
REQUIRE(name(LSPrecond::none) == "none");
REQUIRE(name(LSPrecond::diagonal) == "diagonal");
- REQUIRE(name(LSPrecond::incomplete_choleski) == "ICholeski");
+ REQUIRE(name(LSPrecond::incomplete_cholesky) == "ICholesky");
REQUIRE(name(LSPrecond::incomplete_LU) == "ILU");
REQUIRE(name(LSPrecond::amg) == "AMG");
REQUIRE_THROWS_WITH(name(LSPrecond::LS__end),
@@ -109,7 +109,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]")
REQUIRE(LSMethod::bicgstab == getLSEnumFromName<LSMethod>("BICGStab"));
REQUIRE(LSMethod::bicgstab2 == getLSEnumFromName<LSMethod>("BICGStab2"));
REQUIRE(LSMethod::lu == getLSEnumFromName<LSMethod>("LU"));
- REQUIRE(LSMethod::choleski == getLSEnumFromName<LSMethod>("Choleski"));
+ REQUIRE(LSMethod::cholesky == getLSEnumFromName<LSMethod>("Cholesky"));
REQUIRE(LSMethod::gmres == getLSEnumFromName<LSMethod>("GMRES"));
REQUIRE_THROWS_WITH(getLSEnumFromName<LSMethod>("__invalid_method"),
@@ -120,7 +120,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]")
{
REQUIRE(LSPrecond::none == getLSEnumFromName<LSPrecond>("none"));
REQUIRE(LSPrecond::diagonal == getLSEnumFromName<LSPrecond>("diagonal"));
- REQUIRE(LSPrecond::incomplete_choleski == getLSEnumFromName<LSPrecond>("ICholeski"));
+ REQUIRE(LSPrecond::incomplete_cholesky == getLSEnumFromName<LSPrecond>("ICholesky"));
REQUIRE(LSPrecond::incomplete_LU == getLSEnumFromName<LSPrecond>("ILU"));
REQUIRE_THROWS_WITH(getLSEnumFromName<LSPrecond>("__invalid_precond"),
@@ -153,7 +153,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]")
- BICGStab2
- GMRES
- LU
- - Choleski
+ - Cholesky
)";
REQUIRE(os.str() == library_list);
@@ -168,7 +168,7 @@ TEST_CASE("LinearSolverOptions", "[algebra]")
const std::string library_list = R"(
- none
- diagonal
- - ICholeski
+ - ICholesky
- ILU
- AMG
)";
diff --git a/tests/test_OFStream.cpp b/tests/test_OFStream.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..32b4ff21b8e725dcaff822b586f9810d617d3a4d
--- /dev/null
+++ b/tests/test_OFStream.cpp
@@ -0,0 +1,58 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/utils/OFStream.hpp>
+#include <utils/Messenger.hpp>
+
+#include <filesystem>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("OFStream", "[language]")
+{
+ SECTION("ofstream")
+ {
+ const std::string basename = "ofstream_";
+ const std::string filename = basename + std::to_string(parallel::rank());
+
+ // Ensures that the file is closed after this line
+ std::make_shared<OFStream>(filename) << "foo" << 3 << " bar\n";
+
+ if (parallel::rank() == 0) {
+ REQUIRE(std::filesystem::is_regular_file(filename));
+
+ std::ifstream is(filename);
+
+ char file_content[10];
+ for (size_t i = 0; i < 10; ++i) {
+ char c = is.get();
+
+ file_content[i] = c;
+ if (c == '\n') {
+ file_content[i + 1] = '\0';
+ REQUIRE(i == 8);
+
+ c = is.get();
+ REQUIRE(is.eof());
+ break;
+ }
+ }
+
+ std::string content = file_content;
+ REQUIRE(content == "foo3 bar\n");
+
+ std::filesystem::remove(filename);
+ }
+
+ REQUIRE(not std::filesystem::exists(filename));
+ }
+
+ SECTION("invalid filename")
+ {
+ if (parallel::rank() == 0) {
+ const std::string filename = "badpath/invalidpath/ofstream";
+
+ REQUIRE_THROWS_WITH(std::make_shared<OFStream>(filename), "error: cannot create file " + filename);
+ }
+ }
+}
diff --git a/tests/test_OStream.cpp b/tests/test_OStream.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..417db374ee366531623af397e2d62e78d1c83475
--- /dev/null
+++ b/tests/test_OStream.cpp
@@ -0,0 +1,36 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/utils/OStream.hpp>
+
+#include <sstream>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("OStream", "[language]")
+{
+ SECTION("null ostream")
+ {
+ std::shared_ptr os = std::make_shared<OStream>();
+
+ REQUIRE_NOTHROW(os << "foo" << 3 << " bar");
+ }
+
+ SECTION("valid ostream")
+ {
+ std::stringstream sstr;
+
+ std::shared_ptr os = std::make_shared<OStream>(sstr);
+ os << "foo" << 3 << " bar";
+
+ REQUIRE(sstr.str() == "foo3 bar");
+ }
+
+#ifndef NDEBUG
+ SECTION("non allocated stream")
+ {
+ std::shared_ptr<OStream> bad_os;
+ REQUIRE_THROWS_WITH((bad_os << 2), "non allocated stream");
+ }
+#endif // NDEBUG
+}
diff --git a/tests/test_OStreamProcessor.cpp b/tests/test_OStreamProcessor.cpp
deleted file mode 100644
index 22213d64bcbb70adf6ea6fab93f82e641870da29..0000000000000000000000000000000000000000
--- a/tests/test_OStreamProcessor.cpp
+++ /dev/null
@@ -1,100 +0,0 @@
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_all.hpp>
-
-#include <language/ast/ASTBuilder.hpp>
-#include <language/ast/ASTNodeDataTypeBuilder.hpp>
-#include <language/ast/ASTNodeDeclarationToAffectationConverter.hpp>
-#include <language/ast/ASTNodeExpressionBuilder.hpp>
-#include <language/ast/ASTNodeTypeCleaner.hpp>
-#include <language/ast/ASTSymbolTableBuilder.hpp>
-#include <language/node_processor/OStreamProcessor.hpp>
-#include <utils/Demangle.hpp>
-
-#include <pegtl/string_input.hpp>
-
-#include <sstream>
-
-void
-_replaceOStream(ASTNode& node, std::ostringstream& sout)
-{
- if (node.is_type<language::cout_kw>() or node.is_type<language::cerr_kw>()) {
- node.m_node_processor = std::make_unique<OStreamProcessor>(node, sout);
- } else {
- for (auto& child_node : node.children) {
- _replaceOStream(*child_node, sout);
- }
- }
-}
-
-#define CHECK_OSTREAM_EXPRESSION_RESULT(data, expected_value) \
- { \
- TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"}; \
- auto ast = ASTBuilder::build(input); \
- \
- ASTSymbolTableBuilder{*ast}; \
- ASTNodeDataTypeBuilder{*ast}; \
- \
- ASTNodeDeclarationToAffectationConverter{*ast}; \
- ASTNodeTypeCleaner<language::var_declaration>{*ast}; \
- \
- ASTNodeExpressionBuilder{*ast}; \
- ExecutionPolicy exec_policy; \
- \
- std::ostringstream sout; \
- _replaceOStream(*ast, sout); \
- \
- ast->execute(exec_policy); \
- \
- REQUIRE(sout.str() == expected_value); \
- }
-
-#define CHECK_OSTREAM_EXPRESSION_THROWS(data, expected_error) \
- { \
- static_assert(std::is_same_v<std::decay_t<decltype(data)>, std::string_view>); \
- static_assert((std::is_same_v<std::decay_t<decltype(expected_error)>, std::string_view>) or \
- (std::is_same_v<std::decay_t<decltype(expected_error)>, std::string>)); \
- \
- TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"}; \
- auto ast = ASTBuilder::build(input); \
- \
- ASTSymbolTableBuilder{*ast}; \
- ASTNodeDataTypeBuilder{*ast}; \
- \
- ASTNodeDeclarationToAffectationConverter{*ast}; \
- ASTNodeTypeCleaner<language::var_declaration>{*ast}; \
- ASTNodeTypeCleaner<language::fct_declaration>{*ast}; \
- \
- REQUIRE_THROWS_WITH(ASTNodeExpressionBuilder{*ast}, expected_error); \
- }
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("OStreamProcessor", "[language]")
-{
- SECTION("cout")
- {
- CHECK_OSTREAM_EXPRESSION_RESULT(R"(cout << 2;)", "2");
- CHECK_OSTREAM_EXPRESSION_RESULT(R"(cout << true;)", "true");
- CHECK_OSTREAM_EXPRESSION_RESULT(R"(cout << false;)", "false");
- CHECK_OSTREAM_EXPRESSION_RESULT(R"(cout << "x=" << 2 << "\n";)", "x=2\n");
- }
-
- SECTION("cerr")
- {
- CHECK_OSTREAM_EXPRESSION_RESULT(R"(cerr << 2;)", "2");
- CHECK_OSTREAM_EXPRESSION_RESULT(R"(cerr << true;)", "true");
- CHECK_OSTREAM_EXPRESSION_RESULT(R"(cerr << false;)", "false");
- CHECK_OSTREAM_EXPRESSION_RESULT(R"(cerr << "x=" << 2 << "\n";)", "x=2\n");
- }
-
- SECTION("runtime error")
- {
- std::string_view data_type = R"(
-let f : R -> R, x -> 2;
-cerr << "f=" << f << "\n";
-)";
-
- std::string error_msg = "invalid argument, cannot print a 'function'";
- CHECK_OSTREAM_EXPRESSION_THROWS(data_type, error_msg);
- }
-}
diff --git a/tests/test_PETScUtils.cpp b/tests/test_PETScUtils.cpp
index 2e0daa8519241fde2f7fd1749b40398977b839de..493d2dd445fd62f9712748137c806c4cd206f9b1 100644
--- a/tests/test_PETScUtils.cpp
+++ b/tests/test_PETScUtils.cpp
@@ -29,9 +29,9 @@ TEST_CASE("PETScUtils", "[algebra]")
}
}
- SECTION("from DenseMatrix")
+ SECTION("from SmallMatrix")
{
- DenseMatrix<double> A(3, 3);
+ SmallMatrix<double> A(3, 3);
for (size_t i = 0; i < 3; ++i) {
for (size_t j = 0; j < 3; ++j) {
A(i, j) = 1 + i * 3 + j;
@@ -52,9 +52,9 @@ TEST_CASE("PETScUtils", "[algebra]")
}
}
- SECTION("from DenseMatrix [non-square]")
+ SECTION("from SmallMatrix [non-square]")
{
- DenseMatrix<double> A(4, 3);
+ SmallMatrix<double> A(4, 3);
for (size_t i = 0; i < 4; ++i) {
for (size_t j = 0; j < 3; ++j) {
A(i, j) = 1 + i * 3 + j;
diff --git a/tests/test_SmallArray.cpp b/tests/test_SmallArray.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9b5d482791f9df0a19c41c64cb55070c5bf8be24
--- /dev/null
+++ b/tests/test_SmallArray.cpp
@@ -0,0 +1,285 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <utils/PugsAssert.hpp>
+#include <utils/SmallArray.hpp>
+#include <utils/Types.hpp>
+
+#include <deque>
+#include <list>
+#include <set>
+#include <unordered_set>
+#include <valarray>
+#include <vector>
+
+// Instantiate to ensure full coverage is performed
+template class SmallArray<int>;
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SmallArray", "[utils]")
+{
+ SmallArray<int> a(10);
+ REQUIRE(a.size() == 10);
+
+ for (size_t i = 0; i < a.size(); ++i) {
+ a[i] = 2 * i;
+ }
+
+ REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
+ (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));
+
+ SECTION("checking for copies")
+ {
+ SmallArray<const int> b{a};
+
+ REQUIRE(((b[0] == 0) and (b[1] == 2) and (b[2] == 4) and (b[3] == 6) and (b[4] == 8) and (b[5] == 10) and
+ (b[6] == 12) and (b[7] == 14) and (b[8] == 16) and (b[9] == 18)));
+
+ SmallArray<int> c{a};
+
+ REQUIRE(((c[0] == 0) and (c[1] == 2) and (c[2] == 4) and (c[3] == 6) and (c[4] == 8) and (c[5] == 10) and
+ (c[6] == 12) and (c[7] == 14) and (c[8] == 16) and (c[9] == 18)));
+
+ SmallArray<int> d = std::move(c);
+
+ REQUIRE(((d[0] == 0) and (d[1] == 2) and (d[2] == 4) and (d[3] == 6) and (d[4] == 8) and (d[5] == 10) and
+ (d[6] == 12) and (d[7] == 14) and (d[8] == 16) and (d[9] == 18)));
+ }
+
+ SECTION("checking for fill")
+ {
+ SmallArray<int> b(10);
+ b.fill(3);
+
+ REQUIRE(((b[0] == 3) and (b[1] == 3) and (b[2] == 3) and (b[3] == 3) and (b[4] == 3) and (b[5] == 3) and
+ (b[6] == 3) and (b[7] == 3) and (b[8] == 3) and (b[9] == 3)));
+ }
+
+ SECTION("checking for affectations (shallow copy)")
+ {
+ SmallArray<const int> b;
+ b = a;
+
+ REQUIRE(((b[0] == 0) and (b[1] == 2) and (b[2] == 4) and (b[3] == 6) and (b[4] == 8) and (b[5] == 10) and
+ (b[6] == 12) and (b[7] == 14) and (b[8] == 16) and (b[9] == 18)));
+
+ SmallArray<int> c;
+ c = a;
+
+ REQUIRE(((c[0] == 0) and (c[1] == 2) and (c[2] == 4) and (c[3] == 6) and (c[4] == 8) and (c[5] == 10) and
+ (c[6] == 12) and (c[7] == 14) and (c[8] == 16) and (c[9] == 18)));
+
+ SmallArray<int> d;
+ d = std::move(c);
+
+ REQUIRE(((d[0] == 0) and (d[1] == 2) and (d[2] == 4) and (d[3] == 6) and (d[4] == 8) and (d[5] == 10) and
+ (d[6] == 12) and (d[7] == 14) and (d[8] == 16) and (d[9] == 18)));
+ }
+
+ SECTION("checking for affectations (deep copy)")
+ {
+ SmallArray<int> b(copy(a));
+
+ REQUIRE(((b[0] == 0) and (b[1] == 2) and (b[2] == 4) and (b[3] == 6) and (b[4] == 8) and (b[5] == 10) and
+ (b[6] == 12) and (b[7] == 14) and (b[8] == 16) and (b[9] == 18)));
+
+ b.fill(2);
+
+ REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
+ (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));
+
+ REQUIRE(((b[0] == 2) and (b[1] == 2) and (b[2] == 2) and (b[3] == 2) and (b[4] == 2) and (b[5] == 2) and
+ (b[6] == 2) and (b[7] == 2) and (b[8] == 2) and (b[9] == 2)));
+
+ SmallArray<int> c;
+ c = a;
+
+ REQUIRE(((c[0] == 0) and (c[1] == 2) and (c[2] == 4) and (c[3] == 6) and (c[4] == 8) and (c[5] == 10) and
+ (c[6] == 12) and (c[7] == 14) and (c[8] == 16) and (c[9] == 18)));
+
+ c = copy(b);
+
+ REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
+ (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));
+
+ REQUIRE(((c[0] == 2) and (c[1] == 2) and (c[2] == 2) and (c[3] == 2) and (c[4] == 2) and (c[5] == 2) and
+ (c[6] == 2) and (c[7] == 2) and (c[8] == 2) and (c[9] == 2)));
+
+ SmallArray<int> d{a.size()};
+ copy_to(a, d);
+
+ REQUIRE(((d[0] == 0) and (d[1] == 2) and (d[2] == 4) and (d[3] == 6) and (d[4] == 8) and (d[5] == 10) and
+ (d[6] == 12) and (d[7] == 14) and (d[8] == 16) and (d[9] == 18)));
+
+ REQUIRE(((c[0] == 2) and (c[1] == 2) and (c[2] == 2) and (c[3] == 2) and (c[4] == 2) and (c[5] == 2) and
+ (c[6] == 2) and (c[7] == 2) and (c[8] == 2) and (c[9] == 2)));
+
+ copy_to(c, d);
+
+ REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
+ (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));
+
+ REQUIRE(((d[0] == 2) and (d[1] == 2) and (d[2] == 2) and (d[3] == 2) and (d[4] == 2) and (d[5] == 2) and
+ (d[6] == 2) and (d[7] == 2) and (d[8] == 2) and (d[9] == 2)));
+ }
+
+ SECTION("checking for std container conversion")
+ {
+ {
+ std::vector<int> v{1, 2, 5, 3};
+ {
+ SmallArray<int> v_array = convert_to_small_array(v);
+
+ REQUIRE(v_array.size() == v.size());
+ REQUIRE(((v_array[0] == 1) and (v_array[1] == 2) and (v_array[2] == 5) and (v_array[3] == 3)));
+ }
+
+ {
+ SmallArray<const int> v_array = convert_to_small_array(v);
+
+ REQUIRE(v_array.size() == v.size());
+ REQUIRE(((v_array[0] == 1) and (v_array[1] == 2) and (v_array[2] == 5) and (v_array[3] == 3)));
+ }
+ }
+
+ {
+ std::vector<int> w;
+ {
+ SmallArray<int> w_array = convert_to_small_array(w);
+ REQUIRE(w_array.size() == 0);
+ }
+ {
+ SmallArray<const int> w_array = convert_to_small_array(w);
+ REQUIRE(w_array.size() == 0);
+ }
+ }
+
+ {
+ std::valarray<int> v{1, 2, 5, 3};
+ SmallArray<int> v_array = convert_to_small_array(v);
+
+ REQUIRE(v_array.size() == v.size());
+ REQUIRE(((v_array[0] == 1) and (v_array[1] == 2) and (v_array[2] == 5) and (v_array[3] == 3)));
+ }
+
+ {
+ std::set<int> s{4, 2, 5, 3, 1, 3, 2};
+ SmallArray<int> s_array = convert_to_small_array(s);
+
+ REQUIRE(s_array.size() == s.size());
+ REQUIRE(
+ ((s_array[0] == 1) and (s_array[1] == 2) and (s_array[2] == 3) and (s_array[3] == 4) and (s_array[4] == 5)));
+ }
+
+ {
+ std::unordered_set<int> us{4, 2, 5, 3, 1, 3, 2};
+ SmallArray<int> us_array = convert_to_small_array(us);
+
+ REQUIRE(us_array.size() == us.size());
+
+ std::set<int> s;
+ for (size_t i = 0; i < us_array.size(); ++i) {
+ REQUIRE((us.find(us_array[i]) != us.end()));
+ s.insert(us_array[i]);
+ }
+ REQUIRE(s.size() == us_array.size());
+ }
+
+ {
+ std::multiset<int> ms{4, 2, 5, 3, 1, 3, 2};
+ SmallArray<int> ms_array = convert_to_small_array(ms);
+
+ REQUIRE(ms_array.size() == ms.size());
+ REQUIRE(((ms_array[0] == 1) and (ms_array[1] == 2) and (ms_array[2] == 2) and (ms_array[3] == 3) and
+ (ms_array[4] == 3) and (ms_array[5] == 4) and (ms_array[6] == 5)));
+ }
+
+ {
+ std::list<int> l{1, 3, 5, 6, 2};
+ SmallArray<int> l_array = convert_to_small_array(l);
+
+ REQUIRE(l_array.size() == l.size());
+ REQUIRE(
+ ((l_array[0] == 1) and (l_array[1] == 3) and (l_array[2] == 5) and (l_array[3] == 6) and (l_array[4] == 2)));
+ }
+
+ {
+ std::deque<int> q{1, 3, 5, 6, 2};
+ q.push_front(2);
+ SmallArray<int> q_array = convert_to_small_array(q);
+
+ REQUIRE(q_array.size() == q.size());
+ REQUIRE(((q_array[0] == 2) and (q_array[1] == 1) and (q_array[2] == 3) and (q_array[3] == 5) and
+ (q_array[4] == 6) and (q_array[5] == 2)));
+ }
+ }
+
+ SECTION("output")
+ {
+ SmallArray<int> x{5};
+ x[0] = 2;
+ x[1] = 6;
+ x[2] = 2;
+ x[3] = 3;
+ x[4] = 7;
+
+ std::ostringstream array_ost;
+ array_ost << x;
+ std::ostringstream ref_ost;
+ ref_ost << 0 << ':' << x[0];
+ for (size_t i = 1; i < x.size(); ++i) {
+ ref_ost << ' ' << i << ':' << x[i];
+ }
+ REQUIRE(array_ost.str() == ref_ost.str());
+ }
+
+#ifndef NDEBUG
+
+ SECTION("output with signaling NaN")
+ {
+ SmallArray<double> x{5};
+ x[0] = 2;
+ x[2] = 3;
+
+ std::ostringstream array_ost;
+ array_ost << x;
+ std::ostringstream ref_ost;
+ ref_ost << 0 << ':' << 2 << ' ';
+ ref_ost << 1 << ":nan ";
+ ref_ost << 2 << ':' << 3 << ' ';
+ ref_ost << 3 << ":nan ";
+ ref_ost << 4 << ":nan";
+ REQUIRE(array_ost.str() == ref_ost.str());
+ }
+
+ SECTION("checking for bounds violation")
+ {
+ REQUIRE_THROWS_AS(a[10], AssertError);
+ }
+
+ SECTION("invalid copy_to")
+ {
+ SmallArray<int> b{2 * a.size()};
+ REQUIRE_THROWS_AS(copy_to(a, b), AssertError);
+ }
+
+ SECTION("checking for nan initialization")
+ {
+ SmallArray<double> array(10);
+
+ for (size_t i = 0; i < array.size(); ++i) {
+ REQUIRE(std::isnan(array[i]));
+ }
+ }
+
+ SECTION("checking for bad initialization")
+ {
+ SmallArray<int> array(10);
+
+ for (size_t i = 0; i < array.size(); ++i) {
+ REQUIRE(array[i] == std::numeric_limits<int>::max() / 2);
+ }
+ }
+#endif // NDEBUG
+}
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index f50d2abc61f26bdbebbac4ba14ce884b352ec0aa..c9b72d5de9ecb6b047ae8c16ad832f8b29b38157 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -12,39 +12,41 @@
#include <sstream>
// Instantiate to ensure full coverage is performed
-template class TinyMatrix<1, int>;
-template class TinyMatrix<2, int>;
-template class TinyMatrix<3, int>;
-template class TinyMatrix<4, double>;
+template class TinyMatrix<1, 1, int>;
+template class TinyMatrix<2, 2, int>;
+template class TinyMatrix<3, 4, int>;
+template class TinyMatrix<4, 4, double>;
// clazy:excludeall=non-pod-global-static
TEST_CASE("TinyMatrix", "[algebra]")
{
- TinyMatrix<3, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9);
- REQUIRE(((A(0, 0) == 1) and (A(0, 1) == 2) and (A(0, 2) == 3) and (A(1, 0) == 4) and (A(1, 1) == 5) and
- (A(1, 2) == 6) and (A(2, 0) == 7) and (A(2, 1) == 8) and (A(2, 2) == 9)));
+ TinyMatrix<3, 4, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
+ REQUIRE(((A(0, 0) == 1) and (A(0, 1) == 2) and (A(0, 2) == 3) and (A(0, 3) == 4) and //
+ (A(1, 0) == 5) and (A(1, 1) == 6) and (A(1, 2) == 7) and (A(1, 3) == 8) and //
+ (A(2, 0) == 9) and (A(2, 1) == 10) and (A(2, 2) == 11) and (A(2, 3) == 12)));
- TinyMatrix<3, int> B(6, 5, 3, 8, 34, 6, 35, 6, 7);
+ TinyMatrix<3, 4, int> B(6, 5, 3, 8, 34, 6, 35, 6, 7, 1, 3, 6);
SECTION("checking for opposed matrix")
{
- const TinyMatrix<3, int> minus_A = -A;
- REQUIRE(((minus_A(0, 0) == -1) and (minus_A(0, 1) == -2) and (minus_A(0, 2) == -3) and (minus_A(1, 0) == -4) and
- (minus_A(1, 1) == -5) and (minus_A(1, 2) == -6) and (minus_A(2, 0) == -7) and (minus_A(2, 1) == -8) and
- (minus_A(2, 2) == -9)));
+ const TinyMatrix minus_A = -A;
+ REQUIRE(
+ ((minus_A(0, 0) == -1) and (minus_A(0, 1) == -2) and (minus_A(0, 2) == -3) and (minus_A(0, 3) == -4) and //
+ (minus_A(1, 0) == -5) and (minus_A(1, 1) == -6) and (minus_A(1, 2) == -7) and (minus_A(1, 3) == -8) and //
+ (minus_A(2, 0) == -9) and (minus_A(2, 1) == -10) and (minus_A(2, 2) == -11) and (minus_A(2, 3) == -12)));
}
SECTION("checking for equality and difference tests")
{
- const TinyMatrix<3, int> copy_A = A;
- REQUIRE(((copy_A(0, 0) == 1) and (copy_A(0, 1) == 2) and (copy_A(0, 2) == 3) and (copy_A(1, 0) == 4) and
- (copy_A(1, 1) == 5) and (copy_A(1, 2) == 6) and (copy_A(2, 0) == 7) and (copy_A(2, 1) == 8) and
- (copy_A(2, 2) == 9)));
+ const TinyMatrix copy_A = A;
+ REQUIRE(((copy_A(0, 0) == 1) and (copy_A(0, 1) == 2) and (copy_A(0, 2) == 3) and (copy_A(0, 3) == 4) and //
+ (copy_A(1, 0) == 5) and (copy_A(1, 1) == 6) and (copy_A(1, 2) == 7) and (copy_A(1, 3) == 8) and //
+ (copy_A(2, 0) == 9) and (copy_A(2, 1) == 10) and (copy_A(2, 2) == 11) and (copy_A(2, 3) == 12)));
REQUIRE(copy_A == A);
REQUIRE_FALSE(copy_A != A);
- TinyMatrix<3, int> affected_A;
+ TinyMatrix<3, 4, int> affected_A;
affected_A = A;
REQUIRE(affected_A == A);
REQUIRE_FALSE(affected_A != A);
@@ -55,101 +57,109 @@ TEST_CASE("TinyMatrix", "[algebra]")
SECTION("checking for scalar left product")
{
- const int a = 2;
- const TinyMatrix<3, int> aA = a * A;
+ const int a = 2;
+ const TinyMatrix<3, 4, int> aA = a * A;
- REQUIRE(aA == TinyMatrix<3, int>(2, 4, 6, 8, 10, 12, 14, 16, 18));
+ REQUIRE(aA == TinyMatrix<3, 4, int>(2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24));
}
SECTION("checking for scalar seft product")
{
- const int a = 2;
- TinyMatrix<3, int> copy_A = A;
+ const int a = 2;
- REQUIRE((copy_A *= a) == TinyMatrix<3, int>(2, 4, 6, 8, 10, 12, 14, 16, 18));
+ TinyMatrix copy_A = A;
+
+ REQUIRE((copy_A *= a) == TinyMatrix<3, 4, int>(2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24));
}
SECTION("checking for null matrix management")
{
- TinyMatrix<3, int> Z = zero;
- REQUIRE(Z == TinyMatrix<3, int>(0, 0, 0, 0, 0, 0, 0, 0, 0));
+ TinyMatrix<4, 3, int> Z = zero;
+ REQUIRE(Z == TinyMatrix<4, 3, int>(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
- TinyMatrix<3, int> affected_Z;
+ TinyMatrix<4, 3, int> affected_Z;
affected_Z = zero;
REQUIRE(affected_Z == Z);
}
SECTION("checking for identity management")
{
- TinyMatrix<3, int> I = identity;
- REQUIRE(I == TinyMatrix<3, int>(1, 0, 0, 0, 1, 0, 0, 0, 1));
+ TinyMatrix<3, 3, int> I = identity;
+ REQUIRE(I == TinyMatrix<3, 3, int>(1, 0, 0, 0, 1, 0, 0, 0, 1));
- TinyMatrix<3, int> affected_I;
+ TinyMatrix<3, 3, int> affected_I;
affected_I = identity;
REQUIRE(affected_I == I);
}
SECTION("checking for matrices sum")
{
- REQUIRE(A + B == TinyMatrix<3, int>(7, 7, 6, 12, 39, 12, 42, 14, 16));
+ REQUIRE(A + B == TinyMatrix<3, 4, int>(7, 7, 6, 12, 39, 12, 42, 14, 16, 11, 14, 18));
- TinyMatrix<3, int> ApB = A;
+ TinyMatrix<3, 4, int> ApB = A;
ApB += B;
REQUIRE(ApB == A + B);
- TinyMatrix<3, int> Ap2B = A + 2 * B;
+ TinyMatrix<3, 4, int> Ap2B = A + 2 * B;
REQUIRE(Ap2B == ApB + B);
}
SECTION("checking for matrices difference ")
{
- REQUIRE(A - B == TinyMatrix<3, int>(-5, -3, 0, -4, -29, 0, -28, 2, 2));
+ REQUIRE(A - B == TinyMatrix<3, 4, int>(-5, -3, 0, -4, -29, 0, -28, 2, 2, 9, 8, 6));
- TinyMatrix<3, int> AmB = A;
+ TinyMatrix<3, 4, int> AmB = A;
AmB -= B;
REQUIRE(AmB == A - B);
- TinyMatrix<3, int> Am2B = A - 2 * B;
+ TinyMatrix<3, 4, int> Am2B = A - 2 * B;
REQUIRE(Am2B == AmB - B);
}
SECTION("checking for matrices product")
{
- REQUIRE(A * B == TinyMatrix<3, int>(127, 91, 36, 274, 226, 84, 421, 361, 132));
+ TinyMatrix<4, 2, int> C{3, -2, 2, 6, -2, 5, 7, 2};
+ REQUIRE(A * C == TinyMatrix<3, 2, int>(29, 33, 69, 77, 109, 121));
+
+ TinyMatrix<2, 3, int> D{-3, 2, 3, 2, 5, -7};
+ REQUIRE(D * A == TinyMatrix<2, 4, int>(34, 36, 38, 40, -36, -36, -36, -36));
}
SECTION("checking for matrix-vector product")
{
- REQUIRE(A * TinyVector<3, int>(2, -3, 5) == TinyVector<3, int>(11, 23, 35));
+ REQUIRE(A * TinyVector<4, int>(2, -3, 5, 2) == TinyVector<3, int>(19, 43, 67));
}
SECTION("checking for tensor product")
{
- const TinyVector<3, int> u(1, 3, 7);
+ const TinyVector<4, int> u(1, 3, 7, 5);
const TinyVector<3, int> v(6, 2, -3);
- REQUIRE(tensorProduct(u, v) == TinyMatrix<3, int>(6, 2, -3, 18, 6, -9, 42, 14, -21));
+ REQUIRE(tensorProduct(u, v) == TinyMatrix<4, 3, int>(6, 2, -3, 18, 6, -9, 42, 14, -21, 30, 10, -15));
}
SECTION("checking for minor calculation")
{
- TinyMatrix<3, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9);
- REQUIRE(getMinor(A, 0, 0) == TinyMatrix<2, int>(5, 6, 8, 9));
- REQUIRE(getMinor(A, 1, 0) == TinyMatrix<2, int>(2, 3, 8, 9));
- REQUIRE(getMinor(A, 2, 0) == TinyMatrix<2, int>(2, 3, 5, 6));
-
- REQUIRE(getMinor(A, 0, 1) == TinyMatrix<2, int>(4, 6, 7, 9));
- REQUIRE(getMinor(A, 1, 1) == TinyMatrix<2, int>(1, 3, 7, 9));
- REQUIRE(getMinor(A, 2, 1) == TinyMatrix<2, int>(1, 3, 4, 6));
-
- REQUIRE(getMinor(A, 0, 2) == TinyMatrix<2, int>(4, 5, 7, 8));
- REQUIRE(getMinor(A, 1, 2) == TinyMatrix<2, int>(1, 2, 7, 8));
- REQUIRE(getMinor(A, 2, 2) == TinyMatrix<2, int>(1, 2, 4, 5));
+ REQUIRE(getMinor(A, 0, 0) == TinyMatrix<2, 3, int>(6, 7, 8, 10, 11, 12));
+ REQUIRE(getMinor(A, 1, 0) == TinyMatrix<2, 3, int>(2, 3, 4, 10, 11, 12));
+ REQUIRE(getMinor(A, 2, 0) == TinyMatrix<2, 3, int>(2, 3, 4, 6, 7, 8));
+
+ REQUIRE(getMinor(A, 0, 1) == TinyMatrix<2, 3, int>(5, 7, 8, 9, 11, 12));
+ REQUIRE(getMinor(A, 1, 1) == TinyMatrix<2, 3, int>(1, 3, 4, 9, 11, 12));
+ REQUIRE(getMinor(A, 2, 1) == TinyMatrix<2, 3, int>(1, 3, 4, 5, 7, 8));
+
+ REQUIRE(getMinor(A, 0, 2) == TinyMatrix<2, 3, int>(5, 6, 8, 9, 10, 12));
+ REQUIRE(getMinor(A, 1, 2) == TinyMatrix<2, 3, int>(1, 2, 4, 9, 10, 12));
+ REQUIRE(getMinor(A, 2, 2) == TinyMatrix<2, 3, int>(1, 2, 4, 5, 6, 8));
+
+ REQUIRE(getMinor(A, 0, 3) == TinyMatrix<2, 3, int>(5, 6, 7, 9, 10, 11));
+ REQUIRE(getMinor(A, 1, 3) == TinyMatrix<2, 3, int>(1, 2, 3, 9, 10, 11));
+ REQUIRE(getMinor(A, 2, 3) == TinyMatrix<2, 3, int>(1, 2, 3, 5, 6, 7));
}
SECTION("checking for cofactors")
{
- TinyMatrix<3, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9);
+ TinyMatrix<3, 3, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9);
REQUIRE(cofactor(A, 0, 0) == (5 * 9 - 8 * 6));
REQUIRE(cofactor(A, 1, 0) == -(2 * 9 - 8 * 3));
REQUIRE(cofactor(A, 2, 0) == (2 * 6 - 5 * 3));
@@ -165,26 +175,25 @@ TEST_CASE("TinyMatrix", "[algebra]")
SECTION("checking for determinant calculations")
{
- REQUIRE(det(TinyMatrix<1, int>(6)) == 6);
- REQUIRE(det(TinyMatrix<2, int>(3, 1, -3, 6)) == 21);
- REQUIRE(det(TinyMatrix<3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1);
- REQUIRE(det(B) == -1444);
- REQUIRE(det(TinyMatrix<4, double>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
+ REQUIRE(det(TinyMatrix<1, 1, int>(6)) == 6);
+ REQUIRE(det(TinyMatrix<2, 2, int>(3, 1, -3, 6)) == 21);
+ REQUIRE(det(TinyMatrix<3, 3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1);
+ REQUIRE(det(TinyMatrix<3, 3, int>(6, 5, 3, 8, 34, 6, 35, 6, 7)) == -1444);
+ REQUIRE(det(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
Catch::Approx(6661.455).epsilon(1E-14));
-
- REQUIRE(det(TinyMatrix<4, double>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 0);
+ REQUIRE(det(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 0);
}
SECTION("checking for inverse calculations")
{
{
- const TinyMatrix<1, double> A1(2);
+ const TinyMatrix<1> A1(2);
REQUIRE(inverse(A1)(0, 0) == Catch::Approx(0.5).epsilon(1E-14));
}
{
- const TinyMatrix<2, double> A2(2, 3, 4, 1);
- const TinyMatrix<2, double> I = inverse(A2) * A2;
+ const TinyMatrix<2> A2(2, 3, 4, 1);
+ const TinyMatrix<2> I = inverse(A2) * A2;
REQUIRE(I(0, 0) == Catch::Approx(1).epsilon(1E-14));
REQUIRE(I(0, 1) == Catch::Approx(0).margin(1E-14));
@@ -193,8 +202,8 @@ TEST_CASE("TinyMatrix", "[algebra]")
}
{
- const TinyMatrix<3, double> A3(2, 3, 1, 4, -1, 5, -2, 3, 4);
- const TinyMatrix<3, double> I = inverse(A3) * A3;
+ const TinyMatrix<3> A3(2, 3, 1, 4, -1, 5, -2, 3, 4);
+ const TinyMatrix<3> I = inverse(A3) * A3;
REQUIRE(I(0, 0) == Catch::Approx(1).epsilon(1E-14));
REQUIRE(I(0, 1) == Catch::Approx(0).margin(1E-14));
@@ -223,6 +232,11 @@ TEST_CASE("TinyMatrix", "[algebra]")
REQUIRE(TinyMatrix<3>{}.numberOfColumns() == 3);
REQUIRE(TinyMatrix<3>{}.dimension() == 9);
REQUIRE(TinyMatrix<3>{}.numberOfValues() == 9);
+
+ REQUIRE(TinyMatrix<3, 4>{}.numberOfRows() == 3);
+ REQUIRE(TinyMatrix<3, 4>{}.numberOfColumns() == 4);
+ REQUIRE(TinyMatrix<3, 4>{}.dimension() == 12);
+ REQUIRE(TinyMatrix<3, 4>{}.numberOfValues() == 12);
}
SECTION("is square matrix")
@@ -230,26 +244,40 @@ TEST_CASE("TinyMatrix", "[algebra]")
REQUIRE(TinyMatrix<1>{}.isSquare());
REQUIRE(TinyMatrix<2>{}.isSquare());
REQUIRE(TinyMatrix<3>{}.isSquare());
- REQUIRE(TinyMatrix<4>{}.isSquare());
+
+ REQUIRE_FALSE(TinyMatrix<3, 4>{}.isSquare());
+ }
+
+ SECTION("transpose")
+ {
+ TinyMatrix tA = transpose(A);
+
+ REQUIRE(((tA(0, 0) == 1) and (tA(1, 0) == 2) and (tA(2, 0) == 3) and (tA(3, 0) == 4) and //
+ (tA(0, 1) == 5) and (tA(1, 1) == 6) and (tA(2, 1) == 7) and (tA(3, 1) == 8) and //
+ (tA(0, 2) == 9) and (tA(1, 2) == 10) and (tA(2, 2) == 11) and (tA(3, 2) == 12)));
+
+ TinyMatrix ttA = transpose(tA);
+ REQUIRE(ttA == A);
}
SECTION("checking for matrices output")
{
- REQUIRE(Catch::Detail::stringify(A) == "[(1,2,3)(4,5,6)(7,8,9)]");
- REQUIRE(Catch::Detail::stringify(TinyMatrix<1, int>(7)) == "[(7)]");
+ REQUIRE(Catch::Detail::stringify(A) == "[(1,2,3,4)(5,6,7,8)(9,10,11,12)]");
+ REQUIRE(Catch::Detail::stringify(TinyMatrix<1, 1, int>(7)) == "[(7)]");
}
#ifndef NDEBUG
SECTION("output with signaling NaN")
{
- TinyMatrix<2> A;
+ TinyMatrix<2, 3> A;
A(0, 0) = 1;
+ A(0, 2) = 3;
A(1, 0) = 2;
std::ostringstream A_ost;
A_ost << A;
std::ostringstream ref_ost;
- ref_ost << "[(1,nan)(2,nan)]";
+ ref_ost << "[(1,nan,3)(2,nan,nan)]";
REQUIRE(A_ost.str() == ref_ost.str());
}
@@ -257,19 +285,19 @@ TEST_CASE("TinyMatrix", "[algebra]")
SECTION("checking for bounds violation")
{
REQUIRE_THROWS_AS(A(3, 0), AssertError);
- REQUIRE_THROWS_AS(A(0, 3), AssertError);
+ REQUIRE_THROWS_AS(A(0, 4), AssertError);
REQUIRE_THROWS_AS(getMinor(A, 3, 0), AssertError);
- REQUIRE_THROWS_AS(getMinor(A, 0, 3), AssertError);
+ REQUIRE_THROWS_AS(getMinor(A, 0, 4), AssertError);
- const TinyMatrix<3, int>& constA = A;
+ const TinyMatrix<3, 4, int>& constA = A;
REQUIRE_THROWS_AS(constA(3, 0), AssertError);
- REQUIRE_THROWS_AS(constA(0, 3), AssertError);
+ REQUIRE_THROWS_AS(constA(0, 4), AssertError);
}
SECTION("checking for nan initialization")
{
- TinyMatrix<3, double> B;
+ TinyMatrix<3, 4, double> B;
for (size_t i = 0; i < B.numberOfRows(); ++i) {
for (size_t j = 0; j < B.numberOfColumns(); ++j) {
@@ -280,7 +308,7 @@ TEST_CASE("TinyMatrix", "[algebra]")
SECTION("checking for bad initialization")
{
- TinyMatrix<3, int> B;
+ TinyMatrix<3, 4, int> B;
for (size_t i = 0; i < B.numberOfRows(); ++i) {
for (size_t j = 0; j < B.numberOfColumns(); ++j) {