diff --git a/dune/xt/common/fmatrix.hh b/dune/xt/common/fmatrix.hh
index 65c21efbb17f0eb96ace6752fc43d77ea6e13435..92efd3b3e28c3c46c104303223622249b2abc306 100644
--- a/dune/xt/common/fmatrix.hh
+++ b/dune/xt/common/fmatrix.hh
@@ -386,6 +386,16 @@ struct MatrixAbstraction<Dune::XT::Common::FieldMatrix<K, N, M>>
return MatrixType(rows, cols, val);
}
+ static inline std::unique_ptr<MatrixType> create_dynamic(const size_t rows, const size_t cols)
+ {
+ return std::make_unique<MatrixType>(rows, cols);
+ }
+
+ static inline std::unique_ptr<MatrixType> create_dynamic(const size_t rows, const size_t cols, const ScalarType& val)
+ {
+ return std::make_unique<MatrixType>(rows, cols, val);
+ }
+
static inline size_t rows(const MatrixType& /*mat*/)
{
return N;
@@ -575,6 +585,19 @@ Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS> operator*(const Dune::FieldMatr
return left.rightmultiplyany(right);
}
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+void rightmultiply(Dune::FieldMatrix<K, L_ROWS, R_COLS>& ret,
+ const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
+ const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
+{
+ for (size_t ii = 0; ii < L_ROWS; ++ii) {
+ for (size_t jj = 0; jj < R_COLS; ++jj) {
+ ret[ii][jj] = 0.;
+ for (size_t kk = 0; kk < L_COLS; ++kk)
+ ret[ii][jj] += left[ii][kk] * right[kk][jj];
+ }
+ }
+}
template <class L, int L_ROWS, int L_COLS, class R, int R_COLS>
typename std::enable_if<!std::is_same<L, R>::value,
@@ -587,6 +610,35 @@ typename std::enable_if<!std::is_same<L, R>::value,
return convert_to<Promoted>(left).rightmultiplyany(convert_to<Promoted>(right));
}
+// versions that do not allocate matrices on the stack (for large matrices)
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const std::unique_ptr<Dune::FieldMatrix<K, L_ROWS, L_COLS>>& left,
+ const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
+{
+ auto ret = std::make_unique<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>();
+ rightmultiply(*ret, *left, right);
+ return ret;
+}
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
+ const std::unique_ptr<Dune::FieldMatrix<K, L_COLS, R_COLS>>& right)
+{
+ auto ret = std::make_unique<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>();
+ rightmultiply(*ret, left, *right);
+ return ret;
+}
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const std::unique_ptr<Dune::FieldMatrix<K, L_ROWS, L_COLS>>& left,
+ const std::unique_ptr<Dune::FieldMatrix<K, L_COLS, R_COLS>>& right)
+{
+ return left * *right;
+}
+
} // namespace Dune
diff --git a/dune/xt/common/lapacke.cc b/dune/xt/common/lapacke.cc
index ffc053a7a2392414b4166316f9f2bfbc5cf12154..7cd9c8710fc26a55da32aad574ac155dd1213f8d 100644
--- a/dune/xt/common/lapacke.cc
+++ b/dune/xt/common/lapacke.cc
@@ -320,21 +320,25 @@ void dtrsm(const int layout,
const int ldb)
{
#if HAVE_MKL
- return cblas_dtrsm(static_cast<CBLAS_LAYOUT>(layout),
- static_cast<CBLAS_SIDE>(side),
- static_cast<CBLAS_UPLO>(uplo),
- static_cast<CBLAS_TRANSPOSE>(transa),
- static_cast<CBLAS_DIAG>(diag),
- m,
- n,
- alpha,
- a,
- lda,
- b,
- ldb);
+ cblas_dtrsm(static_cast<CBLAS_LAYOUT>(layout),
+ static_cast<CBLAS_SIDE>(side),
+ static_cast<CBLAS_UPLO>(uplo),
+ static_cast<CBLAS_TRANSPOSE>(transa),
+ static_cast<CBLAS_DIAG>(diag),
+ m,
+ n,
+ alpha,
+ a,
+ lda,
+ b,
+ ldb);
+#ifndef NDEBUG
+ for (size_t ii = 0; ii < m; ++ii)
+ if (std::isnan(b[ii]) || std::isinf(b[ii]))
+ DUNE_THROW(Dune::MathError, "Triangular solve using cblas_dtrsm failed!");
+#endif
#else
DUNE_THROW(Exceptions::dependency_missing, "You are missing CBLAS or the intel mkl, check available() first!");
- return 1;
#endif
}
diff --git a/dune/xt/common/matrix.hh b/dune/xt/common/matrix.hh
index 71252c19eadfad219ffe0595f962f8dde8dcc3ef..5aee8b0e1478d5d7a6393ef7ce03ec5490896f15 100644
--- a/dune/xt/common/matrix.hh
+++ b/dune/xt/common/matrix.hh
@@ -61,6 +61,12 @@ struct MatrixAbstraction
static_assert(AlwaysFalse<MatType>::value, "Do not call me if is_matrix is false!");
}
+ static inline /*std::unique_ptr<MatrixType>*/ void
+ create_dynamic(const size_t /*rows*/, const size_t /*cols*/, const ScalarType& /*val*/)
+ {
+ static_assert(AlwaysFalse<MatType>::value, "Do not call me if is_matrix is false!");
+ }
+
static inline /*size_t*/ void rows(const MatrixType& /*mat*/)
{
static_assert(AlwaysFalse<MatType>::value, "Do not call me if is_matrix is false!");
@@ -119,6 +125,16 @@ struct MatrixAbstraction<Dune::DynamicMatrix<K>>
return MatrixType(rows, cols, val);
}
+ static inline std::unique_ptr<MatrixType> create_dynamic(const size_t rows, const size_t cols)
+ {
+ return std::make_unique<MatrixType>(rows, cols);
+ }
+
+ static inline std::unique_ptr<MatrixType> create_dynamic(const size_t rows, const size_t cols, const ScalarType& val)
+ {
+ return std::make_unique<MatrixType>(rows, cols, val);
+ }
+
static inline size_t rows(const MatrixType& mat)
{
return mat.rows();
@@ -175,6 +191,24 @@ struct MatrixAbstraction<Dune::FieldMatrix<K, N, M>>
return MatrixType(val);
}
+ static inline std::unique_ptr<MatrixType> create_dynamic(const size_t rows, const size_t cols)
+ {
+ if (rows != N)
+ DUNE_THROW(Exceptions::shapes_do_not_match, "rows = " << rows << "\nN = " << int(N));
+ if (cols != M)
+ DUNE_THROW(Exceptions::shapes_do_not_match, "cols = " << cols << "\nM = " << int(M));
+ return std::make_unique<MatrixType>();
+ }
+
+ static inline std::unique_ptr<MatrixType> create_dynamic(const size_t rows, const size_t cols, const ScalarType& val)
+ {
+ if (rows != N)
+ DUNE_THROW(Exceptions::shapes_do_not_match, "rows = " << rows << "\nN = " << int(N));
+ if (cols != M)
+ DUNE_THROW(Exceptions::shapes_do_not_match, "cols = " << cols << "\nM = " << int(M));
+ return std::make_unique<MatrixType>(val);
+ }
+
static inline size_t rows(const MatrixType& /*mat*/)
{
return numeric_cast<size_t>(N);