From 954ef9f3d312c91f59f9eecf380a09fbb6abc95b Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 13 Aug 2026 18:21:50 -0500 Subject: [PATCH 01/15] Add unit tests for dense vector/matrix norms --- tests/numerics/dense_matrix_test.C | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/numerics/dense_matrix_test.C b/tests/numerics/dense_matrix_test.C index 0e7b916128..045ada719c 100644 --- a/tests/numerics/dense_matrix_test.C +++ b/tests/numerics/dense_matrix_test.C @@ -23,6 +23,7 @@ public: LIBMESH_CPPUNIT_TEST_SUITE(DenseMatrixTest); CPPUNIT_TEST(testClassifiers); + CPPUNIT_TEST(testNorms); CPPUNIT_TEST(testOuterProduct); CPPUNIT_TEST(testSVD); CPPUNIT_TEST(testEVDreal); @@ -97,6 +98,26 @@ private: CPPUNIT_ASSERT(isnan(diag)); } + void testNorms() + { + LOG_UNIT_TEST; + + DenseVector a = {3.0, -4.0}; + LIBMESH_ASSERT_FP_EQUAL(a.l1_norm(), 7, TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(a.l2_norm(), 5, TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(a.linfty_norm(), 4, TOLERANCE*TOLERANCE); + + DenseVector b = {1.0, -2.0}; + LIBMESH_ASSERT_FP_EQUAL(b.l1_norm(), 3, TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(b.l2_norm(), std::sqrt(Real(5)), TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(b.linfty_norm(), 2, TOLERANCE*TOLERANCE); + + DenseMatrix a_times_b; // [3 -6; -4 8] + a_times_b.outer_product(a, b); + LIBMESH_ASSERT_FP_EQUAL(a_times_b.l1_norm(), 14, TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(a_times_b.linfty_norm(), 12, TOLERANCE*TOLERANCE); + } + void testOuterProduct() { LOG_UNIT_TEST; From 2ab78014c2f06616002e9251ef3082eaac2933a1 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 13 Aug 2026 18:23:09 -0500 Subject: [PATCH 02/15] Use friend functions to avoid operator() This simplifies the code a bit in DenseVector, and might optimize it a bit in DenseMatrix --- include/numerics/dense_matrix.h | 120 +++++++++++++++----------------- include/numerics/dense_vector.h | 113 +++++++++++++++--------------- 2 files changed, 112 insertions(+), 121 deletions(-) diff --git a/include/numerics/dense_matrix.h b/include/numerics/dense_matrix.h index 044221b308..2d681e2bf4 100644 --- a/include/numerics/dense_matrix.h +++ b/include/numerics/dense_matrix.h @@ -564,6 +564,62 @@ class DenseMatrix : public DenseMatrixBase */ T det(); + /** + * Returns true iff every entry is finite. + */ + friend bool isfinite (const DenseMatrix & var) + { + using std::isfinite; + using libMesh::isfinite; // for T==complex + for (const T & v : var._val) + if (!isfinite(v)) + return false; + return true; + } + + /** + * Returns true iff no entry is NaN and any entry is infinite. + * + * This is arguably inconsistent with our std::complex overload (and + * the C99 Annex G recommendations for _Complex, and C++ + * std::complex arithmetic), which treats mixed (inf,NaN) pairs as + * infinite, but this is probably safer for users. + */ + friend bool isinf (const DenseMatrix & var) + { + using std::isinf; + using libMesh::isinf; // for T==complex + using std::isnan; + using libMesh::isnan; + bool has_inf = false; + for (const T & v : var._val) + { + // NaN anywhere makes us NaN, not inf + if (isnan(v)) + return false; + has_inf = has_inf || isinf(v); + } + return has_inf; + } + + /** + * Returns true iff any entry is NaN. + * + * This is arguably inconsistent with our std::complex overload (and + * the C99 Annex G recommendations for _Complex, and C++ + * std::complex arithmetic), which treats mixed (inf,NaN) pairs as + * infinite, but this is probably safer for users. + */ + friend bool isnan (const DenseMatrix & var) + { + using std::isnan; + using libMesh::isnan; // for T==complex + for (const T & v : var._val) + if (isnan(v)) + return true; + return false; + } + /** * Computes the inverse of the dense matrix (assuming it is invertible) * by first computing the LU decomposition and then performing multiple @@ -1218,70 +1274,6 @@ T DenseMatrix::transpose (const unsigned int i, // } -// A matrix is finite iff every component is -template -bool isfinite (const DenseMatrix & var) -{ - using std::isfinite; - using libMesh::isfinite; // for T==complex - const auto m = var.m(), n = var.n(); - for (auto i : make_range(m)) - for (auto j : make_range(n)) - if (!isfinite(var(i,j))) - return false; - return true; -} - - -// A matrix is infinite iff some component is infinite but no -// component is NaN. -// -// This is arguably inconsistent with our std::complex overload (and -// the C99 Annex G recommendations for _Complex, and C++ std::complex -// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but -// this is probably safer for users. -template -bool isinf (const DenseMatrix & var) -{ - using std::isinf; - using libMesh::isinf; // for T==complex - using std::isnan; - using libMesh::isnan; - const auto m = var.m(), n = var.n(); - bool has_inf = false; - for (auto i : make_range(m)) - for (auto j : make_range(n)) - { - // NaN anywhere makes us NaN, not inf - if (isnan(var(i,j))) - return false; - has_inf = has_inf || isinf(var(i,j)); - } - return has_inf; -} - - - -// A matrix is NaN iff some component is NaN -// -// This is arguably inconsistent with our std::complex overload (and -// the C99 Annex G recommendations for _Complex, and C++ std::complex -// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but -// this is probably safer for users. -template -bool isnan (const DenseMatrix & var) -{ - using std::isnan; - using libMesh::isnan; // for T==complex - const auto m = var.m(), n = var.n(); - for (auto i : make_range(m)) - for (auto j : make_range(n)) - if (isnan(var(i,j))) - return true; - return false; -} - - } // namespace libMesh #ifdef LIBMESH_HAVE_METAPHYSICL diff --git a/include/numerics/dense_vector.h b/include/numerics/dense_vector.h index f5d0f21df0..94fe1266fb 100644 --- a/include/numerics/dense_vector.h +++ b/include/numerics/dense_vector.h @@ -294,6 +294,62 @@ class DenseVector : public DenseVectorBase typename std::vector::const_iterator end() const { return _val.end(); } typename std::vector::iterator end() { return _val.end(); } + /** + * Returns true iff every entry is finite. + */ + friend bool isfinite (const DenseVector & var) + { + using std::isfinite; + using libMesh::isfinite; // for T==complex + for (const T & v : var._val) + if (!isfinite(v)) + return false; + return true; + } + + /** + * Returns true iff no entry is NaN and any entry is infinite. + * + * This is arguably inconsistent with our std::complex overload (and + * the C99 Annex G recommendations for _Complex, and C++ + * std::complex arithmetic), which treats mixed (inf,NaN) pairs as + * infinite, but this is probably safer for users. + */ + friend bool isinf (const DenseVector & var) + { + using std::isinf; + using libMesh::isinf; // for T==complex + using std::isnan; + using libMesh::isnan; + bool has_inf = false; + for (const T & v : var._val) + { + // NaN anywhere makes us NaN, not inf + if (isnan(v)) + return false; + has_inf = has_inf || isinf(v); + } + return has_inf; + } + + /** + * Returns true iff any entry is NaN. + * + * This is arguably inconsistent with our std::complex overload (and + * the C99 Annex G recommendations for _Complex, and C++ + * std::complex arithmetic), which treats mixed (inf,NaN) pairs as + * infinite, but this is probably safer for users. + */ + friend bool isnan (const DenseVector & var) + { + using std::isnan; + using libMesh::isnan; // for T==complex + for (const T & v : var._val) + if (isnan(v)) + return true; + return false; + } + private: /** @@ -722,63 +778,6 @@ void DenseVector::get_principal_subvector (unsigned int sub_n, } -// A vector is finite iff every component is -template -bool isfinite (const DenseVector & var) -{ - using std::isfinite; - using libMesh::isfinite; // for T==complex - for (auto i : index_range(var)) - if (!isfinite(var(i))) - return false; - return true; -} - - -// A vector is infinite iff some component is infinite but no -// component is NaN. -// -// This is arguably inconsistent with our std::complex overload (and -// the C99 Annex G recommendations for _Complex, and C++ std::complex -// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but -// this is probably safer for users. -template -bool isinf (const DenseVector & var) -{ - using std::isinf; - using libMesh::isinf; // for T==complex - using std::isnan; - using libMesh::isnan; - bool has_inf = false; - for (auto i : index_range(var)) - { - // NaN anywhere makes us NaN, not inf - if (isnan(var(i))) - return false; - has_inf = has_inf || isinf(var(i)); - } - return has_inf; -} - - -// A vector is NaN iff some component is NaN -// -// This is arguably inconsistent with our std::complex overload (and -// the C99 Annex G recommendations for _Complex, and C++ std::complex -// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but -// this is probably safer for users. -template -bool isnan (const DenseVector & var) -{ - using std::isnan; - using libMesh::isnan; // for T==complex - for (auto i : index_range(var)) - if (isnan(var(i))) - return true; - return false; -} - - } // namespace libMesh #ifdef LIBMESH_HAVE_METAPHYSICL From f96443168651bff465687a3a262bc16e936cc7e2 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 13 Aug 2026 18:24:27 -0500 Subject: [PATCH 03/15] Using declaration for std::swap We'd like to support user types as type T, but user types are supposed to use ADL for these things, not put things in namespace std. --- include/numerics/dense_matrix.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/numerics/dense_matrix.h b/include/numerics/dense_matrix.h index 2d681e2bf4..70efd2d3ac 100644 --- a/include/numerics/dense_matrix.h +++ b/include/numerics/dense_matrix.h @@ -920,8 +920,9 @@ template inline void DenseMatrix::swap(DenseMatrix & other_matrix) { - std::swap(this->_m, other_matrix._m); - std::swap(this->_n, other_matrix._n); + using std::swap; + swap(this->_m, other_matrix._m); + swap(this->_n, other_matrix._n); _val.swap(other_matrix._val); DecompositionType _temp = _decomposition_type; _decomposition_type = other_matrix._decomposition_type; From 70ae621d271f1e0868df435a61fb07cf92603134 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 13 Aug 2026 18:27:11 -0500 Subject: [PATCH 04/15] Use assign() for copy + assignment --- include/numerics/dense_vector.h | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/include/numerics/dense_vector.h b/include/numerics/dense_vector.h index 94fe1266fb..315fb0bde1 100644 --- a/include/numerics/dense_vector.h +++ b/include/numerics/dense_vector.h @@ -388,13 +388,7 @@ DenseVector::DenseVector (const DenseVector & other_vector) : { const std::vector & other_vals = other_vector.get_values(); - _val.clear(); - - const int N = cast_int(other_vals.size()); - _val.reserve(N); - - for (int i=0; i & DenseVector::operator = (const DenseVector & other_vector) { const std::vector & other_vals = other_vector.get_values(); - - _val.clear(); - - const int N = cast_int(other_vals.size()); - _val.reserve(N); - - for (int i=0; i Date: Thu, 13 Aug 2026 18:27:44 -0500 Subject: [PATCH 05/15] range-for saves a line here --- include/numerics/dense_vector.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/include/numerics/dense_vector.h b/include/numerics/dense_vector.h index 315fb0bde1..e69d40ebdb 100644 --- a/include/numerics/dense_vector.h +++ b/include/numerics/dense_vector.h @@ -494,9 +494,8 @@ template inline void DenseVector::scale (const T factor) { - const int N = cast_int(_val.size()); - for (int i=0; i Date: Thu, 13 Aug 2026 18:28:25 -0500 Subject: [PATCH 06/15] Use transform_reduce in DenseMatrix/DenseVector In theory it's possible for a seq (default) STL reduce to be slower than a hand-coded loop; in practice I've found one case (plain integer accumulation) where g++ slows down by 20% (though clang++ gets it right), and for pretty much everything else the STL is like 50% faster or more, even in sequential mode, if only because using "reduce()" tells the compiler that we don't care about element ordering and so it can do SIMD tricks that treat FP arithmetic as associative. --- include/numerics/dense_matrix.h | 33 +++++----------- include/numerics/dense_vector.h | 69 ++++++++++++--------------------- 2 files changed, 34 insertions(+), 68 deletions(-) diff --git a/include/numerics/dense_matrix.h b/include/numerics/dense_matrix.h index 70efd2d3ac..de380ce2ee 100644 --- a/include/numerics/dense_matrix.h +++ b/include/numerics/dense_matrix.h @@ -42,9 +42,10 @@ #endif // C++ includes -#include #include #include +#include +#include #ifdef LIBMESH_HAVE_METAPHYSICL #include "metaphysicl/dualnumber_decl.h" @@ -1141,17 +1142,10 @@ auto DenseMatrix::min () const -> decltype(libmesh_real(T(0))) { libmesh_assert (this->_m); libmesh_assert (this->_n); - auto my_min = libmesh_real((*this)(0,0)); - - for (unsigned int i=0; i!=this->_m; i++) - { - for (unsigned int j=0; j!=this->_n; j++) - { - auto current = libmesh_real((*this)(i,j)); - my_min = (my_min < current? my_min : current); - } - } - return my_min; + return std::transform_reduce + (_val.begin(), _val.end(), std::numeric_limits::max(), + [](auto a, auto b){using std::min; return min(a,b);}, + [](const T & v){return libmesh_real(v);}); } @@ -1162,17 +1156,10 @@ auto DenseMatrix::max () const -> decltype(libmesh_real(T(0))) { libmesh_assert (this->_m); libmesh_assert (this->_n); - auto my_max = libmesh_real((*this)(0,0)); - - for (unsigned int i=0; i!=this->_m; i++) - { - for (unsigned int j=0; j!=this->_n; j++) - { - auto current = libmesh_real((*this)(i,j)); - my_max = (my_max > current? my_max : current); - } - } - return my_max; + return std::transform_reduce + (_val.begin(), _val.end(), std::numeric_limits::lowest(), + [](auto a, auto b){using std::max; return max(a,b);}, + [](const T & v){return libmesh_real(v);}); } diff --git a/include/numerics/dense_vector.h b/include/numerics/dense_vector.h index e69d40ebdb..abea4510e9 100644 --- a/include/numerics/dense_vector.h +++ b/include/numerics/dense_vector.h @@ -38,8 +38,10 @@ #endif // C++ includes -#include +#include #include +#include +#include namespace libMesh { @@ -647,15 +649,11 @@ inline Real DenseVector::min () const { libmesh_assert (this->size()); - Real my_min = libmesh_real((*this)(0)); - - const int N = cast_int(_val.size()); - for (int i=1; i!=N; i++) - { - Real current = libmesh_real((*this)(i)); - my_min = (my_min < current? my_min : current); - } - return my_min; + typedef decltype(libmesh_real(_val[0])) realT; + return std::transform_reduce + (_val.begin(), _val.end(), std::numeric_limits::max(), + [](auto a, auto b){using std::min; return min(a,b);}, + [](const T & v){return libmesh_real(v);}); } @@ -665,15 +663,11 @@ inline Real DenseVector::max () const { libmesh_assert (this->size()); - Real my_max = libmesh_real((*this)(0)); - - const int N = cast_int(_val.size()); - for (int i=1; i!=N; i++) - { - Real current = libmesh_real((*this)(i)); - my_max = (my_max > current? my_max : current); - } - return my_max; + typedef decltype(libmesh_real(_val[0])) realT; + return std::transform_reduce + (_val.begin(), _val.end(), std::numeric_limits::lowest(), + [](auto a, auto b){using std::max; return max(a,b);}, + [](const T & v){return libmesh_real(v);}); } @@ -688,12 +682,9 @@ Real DenseVector::l1_norm () const #ifdef LIBMESH_HAVE_EIGEN return Eigen::Map>(_val.data(), _val.size()).template lpNorm<1>(); #else - Real my_norm = 0.; - const int N = cast_int(_val.size()); - for (int i=0; i!=N; i++) - my_norm += std::abs((*this)(i)); - - return my_norm; + return std::transform_reduce + (_val.begin(), _val.end(), Real(0), std::plus<>(), + [](const T & v){using std::abs; return abs(v);}); #endif } @@ -709,17 +700,10 @@ Real DenseVector::l2_norm () const #ifdef LIBMESH_HAVE_EIGEN return Eigen::Map>(_val.data(), _val.size()).norm(); #else - Real my_norm = 0.; - const int N = cast_int(_val.size()); - // The following pragma tells clang's vectorizer that it is safe to - // reorder floating point operations for this loop. -#ifdef __clang__ -#pragma clang loop vectorize(enable) -#endif - for (int i=0; i!=N; i++) - my_norm += TensorTools::norm_sq((*this)(i)); - - return sqrt(my_norm); + using std::sqrt; + return sqrt(std::transform_reduce + (_val.begin(), _val.end(), Real(0), std::plus<>(), + [](const T & v){return TensorTools::norm_sq(v);})); #endif } @@ -735,15 +719,10 @@ Real DenseVector::linfty_norm () const #ifdef LIBMESH_HAVE_EIGEN return Eigen::Map>(_val.data(), _val.size()).template lpNorm(); #else - Real my_norm = TensorTools::norm_sq((*this)(0)); - - const int N = cast_int(_val.size()); - for (int i=1; i!=N; i++) - { - Real current = TensorTools::norm_sq((*this)(i)); - my_norm = (my_norm > current? my_norm : current); - } - return sqrt(my_norm); + return std::transform_reduce + (_val.begin(), _val.end(), Real(0), + [](auto a, auto b){using std::max; return max(a,b);}, + [](const T & v){using std::abs; return abs(v);}); #endif } From f3fb6161450d40cf582f80d04d3c93032c1901fb Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 14 Aug 2026 12:11:59 -0500 Subject: [PATCH 07/15] Ignore Eigen warnings in parallel_eigen.h Otherwise our tests trip one for me when built with -march=native --- include/parallel/parallel_eigen.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/include/parallel/parallel_eigen.h b/include/parallel/parallel_eigen.h index 5d06a9fff2..ef5f28dd45 100644 --- a/include/parallel/parallel_eigen.h +++ b/include/parallel/parallel_eigen.h @@ -28,8 +28,10 @@ #ifdef LIBMESH_HAVE_EIGEN -// libEigen includes +// libEigen includes, avoiding warnings triggered by some versions +#include "libmesh/ignore_warnings.h" #include +#include "libmesh/restore_warnings.h" namespace libMesh { From 97e8514df2dd7a349dc6879d1a6d4d8b35d2e29a Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Sun, 16 Aug 2026 13:41:45 -0500 Subject: [PATCH 08/15] Pass by const ref in new lambdas This might be a big deal later if we use e.g. AD types for T. --- include/numerics/dense_matrix.h | 4 ++-- include/numerics/dense_vector.h | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/numerics/dense_matrix.h b/include/numerics/dense_matrix.h index de380ce2ee..3f12101c19 100644 --- a/include/numerics/dense_matrix.h +++ b/include/numerics/dense_matrix.h @@ -1144,7 +1144,7 @@ auto DenseMatrix::min () const -> decltype(libmesh_real(T(0))) libmesh_assert (this->_n); return std::transform_reduce (_val.begin(), _val.end(), std::numeric_limits::max(), - [](auto a, auto b){using std::min; return min(a,b);}, + [](const auto & a, const auto & b){using std::min; return min(a,b);}, [](const T & v){return libmesh_real(v);}); } @@ -1158,7 +1158,7 @@ auto DenseMatrix::max () const -> decltype(libmesh_real(T(0))) libmesh_assert (this->_n); return std::transform_reduce (_val.begin(), _val.end(), std::numeric_limits::lowest(), - [](auto a, auto b){using std::max; return max(a,b);}, + [](const auto & a, const auto & b){using std::max; return max(a,b);}, [](const T & v){return libmesh_real(v);}); } diff --git a/include/numerics/dense_vector.h b/include/numerics/dense_vector.h index abea4510e9..e926e94d4b 100644 --- a/include/numerics/dense_vector.h +++ b/include/numerics/dense_vector.h @@ -652,7 +652,7 @@ Real DenseVector::min () const typedef decltype(libmesh_real(_val[0])) realT; return std::transform_reduce (_val.begin(), _val.end(), std::numeric_limits::max(), - [](auto a, auto b){using std::min; return min(a,b);}, + [](const auto & a, const auto & b){using std::min; return min(a,b);}, [](const T & v){return libmesh_real(v);}); } @@ -666,7 +666,7 @@ Real DenseVector::max () const typedef decltype(libmesh_real(_val[0])) realT; return std::transform_reduce (_val.begin(), _val.end(), std::numeric_limits::lowest(), - [](auto a, auto b){using std::max; return max(a,b);}, + [](const auto & a, const auto & b){using std::max; return max(a,b);}, [](const T & v){return libmesh_real(v);}); } From 901d1023f86c6f73cfc0d846674e7683ec88d770 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Sun, 16 Aug 2026 14:40:21 -0500 Subject: [PATCH 09/15] Fix new algorithms for DenseMatrix --- include/numerics/dense_matrix.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/include/numerics/dense_matrix.h b/include/numerics/dense_matrix.h index 3f12101c19..1fc1b5c391 100644 --- a/include/numerics/dense_matrix.h +++ b/include/numerics/dense_matrix.h @@ -1142,8 +1142,9 @@ auto DenseMatrix::min () const -> decltype(libmesh_real(T(0))) { libmesh_assert (this->_m); libmesh_assert (this->_n); + typedef decltype(libmesh_real(T(0))) realfromT; return std::transform_reduce - (_val.begin(), _val.end(), std::numeric_limits::max(), + (_val.begin(), _val.end(), std::numeric_limits::max(), [](const auto & a, const auto & b){using std::min; return min(a,b);}, [](const T & v){return libmesh_real(v);}); } @@ -1156,8 +1157,9 @@ auto DenseMatrix::max () const -> decltype(libmesh_real(T(0))) { libmesh_assert (this->_m); libmesh_assert (this->_n); + typedef decltype(libmesh_real(T(0))) realfromT; return std::transform_reduce - (_val.begin(), _val.end(), std::numeric_limits::lowest(), + (_val.begin(), _val.end(), std::numeric_limits::lowest(), [](const auto & a, const auto & b){using std::max; return max(a,b);}, [](const T & v){return libmesh_real(v);}); } From d8b5db1517b899ad4bf38b213a1a49087c358fc8 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Sun, 16 Aug 2026 14:41:21 -0500 Subject: [PATCH 10/15] Use the same idiom for DenseVector I like this version better. --- include/numerics/dense_vector.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/numerics/dense_vector.h b/include/numerics/dense_vector.h index e926e94d4b..d9ceb74ea3 100644 --- a/include/numerics/dense_vector.h +++ b/include/numerics/dense_vector.h @@ -649,9 +649,9 @@ inline Real DenseVector::min () const { libmesh_assert (this->size()); - typedef decltype(libmesh_real(_val[0])) realT; + typedef decltype(libmesh_real(T(0))) realfromT; return std::transform_reduce - (_val.begin(), _val.end(), std::numeric_limits::max(), + (_val.begin(), _val.end(), std::numeric_limits::max(), [](const auto & a, const auto & b){using std::min; return min(a,b);}, [](const T & v){return libmesh_real(v);}); } @@ -663,9 +663,9 @@ inline Real DenseVector::max () const { libmesh_assert (this->size()); - typedef decltype(libmesh_real(_val[0])) realT; + typedef decltype(libmesh_real(T(0))) realfromT; return std::transform_reduce - (_val.begin(), _val.end(), std::numeric_limits::lowest(), + (_val.begin(), _val.end(), std::numeric_limits::lowest(), [](const auto & a, const auto & b){using std::max; return max(a,b);}, [](const T & v){return libmesh_real(v);}); } From 63350cddc3e7bc0a6103c2b605a4ea9890a375b1 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Sun, 16 Aug 2026 15:09:35 -0500 Subject: [PATCH 11/15] Add autoconf test for std::transform_reduce() --- m4/acsm_cxx_tests.m4 | 39 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/m4/acsm_cxx_tests.m4 b/m4/acsm_cxx_tests.m4 index 8cf688dfca..6ce1512027 100644 --- a/m4/acsm_cxx_tests.m4 +++ b/m4/acsm_cxx_tests.m4 @@ -168,11 +168,46 @@ AC_DEFUN([ACSM_TEST_CXX_ALL], [AC_MSG_WARN([libMesh requires C++11 support for std::isinf]) have_cxx_all=no]) -dnl Optional test here - requiring it would force us to bump up clang -dnl requirements too high for now. +dnl Optional tests here - requiring it would force us to bump up clang +dnl and/or gcc requirements too high for now. LIBMESH_TEST_CXX17_SPLICING AS_IF([test "x$have_cxx17_splicing" != "xyes"], [AC_MSG_WARN([libMesh prefers C++17 support for set/map merge])]) + + LIBMESH_TEST_CXX17_TRANSFORM_REDUCE + AS_IF([test "x$have_cxx17_transform_reduce" != "xyes"], + [AC_MSG_WARN([libMesh prefers C++17 support for std::transform_reduce])]) + ]) + + +dnl Test C++17 std::transform_reduce()" +AC_DEFUN([LIBMESH_TEST_CXX17_TRANSFORM_REDUCE], + [ + have_cxx17_transform_reduce=no + + AC_LANG_PUSH([C++]) + + old_CXXFLAGS="$CXXFLAGS" + CXXFLAGS="$CXXFLAGS $libmesh_CXXFLAGS" + + AC_MSG_CHECKING(for C++17 std::transform_reduce) + AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ + @%:@include + @%:@include + ]], [[ + std::vector v{7,8,9}; + auto sumsq = std::transform_reduce(v.cbegin(), v.cend(), 0, std::plus{}, [](auto x){return x*x;}); + ]])],[ + AC_MSG_RESULT(yes) + have_cxx17_transform_reduce=yes + AC_DEFINE(HAVE_CXX17_TRANSFORM_REDUCE, 1, [Flag indicating whether compiler supports std::transform_reduce]) + ],[ + AC_MSG_RESULT(no) + ]) + + dnl Reset the flags + CXXFLAGS="$old_CXXFLAGS" + AC_LANG_POP([C++]) ]) From 910f7e1c06cd7d7fa439dfe9ac8908fe632d4ece Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Sun, 16 Aug 2026 15:17:44 -0500 Subject: [PATCH 12/15] Add libmesh_transform_reduce shim --- include/base/libmesh_common.h | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/include/base/libmesh_common.h b/include/base/libmesh_common.h index a95b0c4d3d..085c38441d 100644 --- a/include/base/libmesh_common.h +++ b/include/base/libmesh_common.h @@ -569,8 +569,7 @@ struct casting_compare { template inline void libmesh_ignore( const Args&... ) { } -// A workaround for the lack of C++17 merge() support in some -// compilers +// Workarounds for incomplete C++17 support in some compilers/libs #ifdef LIBMESH_HAVE_CXX17_SPLICING template @@ -587,6 +586,26 @@ void libmesh_merge_move(T & target, T & source) } #endif // LIBMESH_HAVE_CXX17_SPLICING +#ifdef LIBMESH_HAVE_CXX17_TRANSFORM_REDUCE +template +T libmesh_transform_reduce(InputIt begin, InputIt end, T init, BinaryOp reduce, UnaryOp transform) +{ + return std::transform_reduce(begin, end, init, reduce, transform); +} +#else +template +T libmesh_transform_reduce(InputIt begin, InputIt end, T init, BinaryOp reduce, UnaryOp transform) +{ + // Reuse init as returnval. + // + // Don't try to do any fancy reduce() reordering/vectorization; + // people who want that can get a real compiler. + for (auto it = begin; it != end; ++it) + init = reduce(init, transform(*it)); + return init; +} +#endif // LIBMESH_HAVE_CXX17_SPLICING + /** * Mostly system independent demangler */ From 420bc6d6a4ed1501b0605d5e915c6c1754dd57cb Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Sun, 16 Aug 2026 15:18:55 -0500 Subject: [PATCH 13/15] Use libmesh_transform_reduce --- include/numerics/dense_matrix.h | 4 ++-- include/numerics/dense_vector.h | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/include/numerics/dense_matrix.h b/include/numerics/dense_matrix.h index 1fc1b5c391..18b18b71fa 100644 --- a/include/numerics/dense_matrix.h +++ b/include/numerics/dense_matrix.h @@ -1143,7 +1143,7 @@ auto DenseMatrix::min () const -> decltype(libmesh_real(T(0))) libmesh_assert (this->_m); libmesh_assert (this->_n); typedef decltype(libmesh_real(T(0))) realfromT; - return std::transform_reduce + return libmesh_transform_reduce (_val.begin(), _val.end(), std::numeric_limits::max(), [](const auto & a, const auto & b){using std::min; return min(a,b);}, [](const T & v){return libmesh_real(v);}); @@ -1158,7 +1158,7 @@ auto DenseMatrix::max () const -> decltype(libmesh_real(T(0))) libmesh_assert (this->_m); libmesh_assert (this->_n); typedef decltype(libmesh_real(T(0))) realfromT; - return std::transform_reduce + return libmesh_transform_reduce (_val.begin(), _val.end(), std::numeric_limits::lowest(), [](const auto & a, const auto & b){using std::max; return max(a,b);}, [](const T & v){return libmesh_real(v);}); diff --git a/include/numerics/dense_vector.h b/include/numerics/dense_vector.h index d9ceb74ea3..7d3154445f 100644 --- a/include/numerics/dense_vector.h +++ b/include/numerics/dense_vector.h @@ -650,7 +650,7 @@ Real DenseVector::min () const { libmesh_assert (this->size()); typedef decltype(libmesh_real(T(0))) realfromT; - return std::transform_reduce + return libmesh_transform_reduce (_val.begin(), _val.end(), std::numeric_limits::max(), [](const auto & a, const auto & b){using std::min; return min(a,b);}, [](const T & v){return libmesh_real(v);}); @@ -664,7 +664,7 @@ Real DenseVector::max () const { libmesh_assert (this->size()); typedef decltype(libmesh_real(T(0))) realfromT; - return std::transform_reduce + return libmesh_transform_reduce (_val.begin(), _val.end(), std::numeric_limits::lowest(), [](const auto & a, const auto & b){using std::max; return max(a,b);}, [](const T & v){return libmesh_real(v);}); @@ -682,7 +682,7 @@ Real DenseVector::l1_norm () const #ifdef LIBMESH_HAVE_EIGEN return Eigen::Map>(_val.data(), _val.size()).template lpNorm<1>(); #else - return std::transform_reduce + return libmesh_transform_reduce (_val.begin(), _val.end(), Real(0), std::plus<>(), [](const T & v){using std::abs; return abs(v);}); #endif @@ -701,7 +701,7 @@ Real DenseVector::l2_norm () const return Eigen::Map>(_val.data(), _val.size()).norm(); #else using std::sqrt; - return sqrt(std::transform_reduce + return sqrt(libmesh_transform_reduce (_val.begin(), _val.end(), Real(0), std::plus<>(), [](const T & v){return TensorTools::norm_sq(v);})); #endif @@ -719,7 +719,7 @@ Real DenseVector::linfty_norm () const #ifdef LIBMESH_HAVE_EIGEN return Eigen::Map>(_val.data(), _val.size()).template lpNorm(); #else - return std::transform_reduce + return libmesh_transform_reduce (_val.begin(), _val.end(), Real(0), [](auto a, auto b){using std::max; return max(a,b);}, [](const T & v){using std::abs; return abs(v);}); From b5e98fc04432821bda9ecf2b06ac747d98e57e08 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Sun, 16 Aug 2026 15:35:04 -0500 Subject: [PATCH 14/15] libmesh_transform_reduce -> parallel_algorithms.h We don't want to pull into every single libMesh file... yet! --- include/base/libmesh_common.h | 22 +---------- include/include_HEADERS | 1 + include/libmesh/Makefile.am | 4 ++ include/numerics/dense_matrix.h | 2 +- include/numerics/dense_vector.h | 4 +- include/parallel/parallel_algorithms.h | 52 ++++++++++++++++++++++++++ 6 files changed, 61 insertions(+), 24 deletions(-) create mode 100644 include/parallel/parallel_algorithms.h diff --git a/include/base/libmesh_common.h b/include/base/libmesh_common.h index 085c38441d..40c96202b6 100644 --- a/include/base/libmesh_common.h +++ b/include/base/libmesh_common.h @@ -569,7 +569,7 @@ struct casting_compare { template inline void libmesh_ignore( const Args&... ) { } -// Workarounds for incomplete C++17 support in some compilers/libs +// Workaround for incomplete C++17 support in some compilers/libs #ifdef LIBMESH_HAVE_CXX17_SPLICING template @@ -586,26 +586,6 @@ void libmesh_merge_move(T & target, T & source) } #endif // LIBMESH_HAVE_CXX17_SPLICING -#ifdef LIBMESH_HAVE_CXX17_TRANSFORM_REDUCE -template -T libmesh_transform_reduce(InputIt begin, InputIt end, T init, BinaryOp reduce, UnaryOp transform) -{ - return std::transform_reduce(begin, end, init, reduce, transform); -} -#else -template -T libmesh_transform_reduce(InputIt begin, InputIt end, T init, BinaryOp reduce, UnaryOp transform) -{ - // Reuse init as returnval. - // - // Don't try to do any fancy reduce() reordering/vectorization; - // people who want that can get a real compiler. - for (auto it = begin; it != end; ++it) - init = reduce(init, transform(*it)); - return init; -} -#endif // LIBMESH_HAVE_CXX17_SPLICING - /** * Mostly system independent demangler */ diff --git a/include/include_HEADERS b/include/include_HEADERS index 05232af00c..55a5e4dfb8 100644 --- a/include/include_HEADERS +++ b/include/include_HEADERS @@ -308,6 +308,7 @@ include_HEADERS = \ parallel/libmesh_call_mpi.h \ parallel/parallel.h \ parallel/parallel_algebra.h \ + parallel/parallel_algorithms.h \ parallel/parallel_bin_sorter.h \ parallel/parallel_eigen.h \ parallel/parallel_elem.h \ diff --git a/include/libmesh/Makefile.am b/include/libmesh/Makefile.am index 03d58b3881..cb337d5ce8 100644 --- a/include/libmesh/Makefile.am +++ b/include/libmesh/Makefile.am @@ -301,6 +301,7 @@ BUILT_SOURCES = \ libmesh_call_mpi.h \ parallel.h \ parallel_algebra.h \ + parallel_algorithms.h \ parallel_bin_sorter.h \ parallel_conversion_utils.h \ parallel_eigen.h \ @@ -1505,6 +1506,9 @@ parallel.h: $(top_srcdir)/include/parallel/parallel.h parallel_algebra.h: $(top_srcdir)/include/parallel/parallel_algebra.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ +parallel_algorithms.h: $(top_srcdir)/include/parallel/parallel_algorithms.h + $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ + parallel_bin_sorter.h: $(top_srcdir)/include/parallel/parallel_bin_sorter.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ diff --git a/include/numerics/dense_matrix.h b/include/numerics/dense_matrix.h index 18b18b71fa..0a7797454e 100644 --- a/include/numerics/dense_matrix.h +++ b/include/numerics/dense_matrix.h @@ -24,6 +24,7 @@ #include "libmesh/libmesh_common.h" #include "libmesh/dense_matrix_base.h" #include "libmesh/int_range.h" +#include "libmesh/parallel_algorithms.h" // For the definition of PetscBLASInt. #if (LIBMESH_HAVE_PETSC) @@ -44,7 +45,6 @@ // C++ includes #include #include -#include #include #ifdef LIBMESH_HAVE_METAPHYSICL diff --git a/include/numerics/dense_vector.h b/include/numerics/dense_vector.h index 7d3154445f..28f2881035 100644 --- a/include/numerics/dense_vector.h +++ b/include/numerics/dense_vector.h @@ -22,9 +22,10 @@ // Local Includes #include "libmesh/libmesh_common.h" -#include "libmesh/dense_vector_base.h" #include "libmesh/compare_types.h" +#include "libmesh/dense_vector_base.h" #include "libmesh/int_range.h" +#include "libmesh/parallel_algorithms.h" #include "libmesh/tensor_tools.h" #ifdef LIBMESH_HAVE_EIGEN @@ -40,7 +41,6 @@ // C++ includes #include #include -#include #include namespace libMesh diff --git a/include/parallel/parallel_algorithms.h b/include/parallel/parallel_algorithms.h new file mode 100644 index 0000000000..e66508306d --- /dev/null +++ b/include/parallel/parallel_algorithms.h @@ -0,0 +1,52 @@ + +// The libMesh Finite Element Library. +// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. + +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. + +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + + +#ifndef LIBMESH_PARALLEL_ALGORITHMS_H +#define LIBMESH_PARALLEL_ALGORITHMS_H + +// The library configuration options +#include "libmesh/libmesh_config.h" + +// C++ headers +#include + +// Workaround incomplete C++17 support in some compilers/libs + +#ifdef LIBMESH_HAVE_CXX17_TRANSFORM_REDUCE +template +T libmesh_transform_reduce(InputIt begin, InputIt end, T init, BinaryOp reduce, UnaryOp transform) +{ + return std::transform_reduce(begin, end, init, reduce, transform); +} +#else +template +T libmesh_transform_reduce(InputIt begin, InputIt end, T init, BinaryOp reduce, UnaryOp transform) +{ + // Reuse init as returnval. + // + // Don't try to do any fancy reduce() reordering/vectorization; + // people who want that can get a real compiler. + for (auto it = begin; it != end; ++it) + init = reduce(init, transform(*it)); + return init; +} +#endif // LIBMESH_HAVE_CXX17_TRANSFORM_REDUCE + +#endif // LIBMESH_PARALLEL_ALGORITHMS_H From 33cf8b9b977eae6bc753b0819b4506cd771e028f Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Sun, 16 Aug 2026 15:11:02 -0500 Subject: [PATCH 15/15] Re-bootstrap --- configure | 325 ++++++++++++++++++++++++++++++++++++ include/Makefile.in | 1 + include/libmesh/Makefile.in | 17 +- include/libmesh_config.h.in | 3 + 4 files changed, 339 insertions(+), 7 deletions(-) diff --git a/configure b/configure index 9aaa955b76..3eae7ecd07 100755 --- a/configure +++ b/configure @@ -15497,6 +15497,71 @@ then : printf "%s\n" "$as_me: WARNING: libMesh prefers C++17 support for set/map merge" >&2;} fi + + have_cxx17_transform_reduce=no + + ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + + old_CXXFLAGS="$CXXFLAGS" + CXXFLAGS="$CXXFLAGS $libmesh_CXXFLAGS" + + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking for C++17 std::transform_reduce" >&5 +printf %s "checking for C++17 std::transform_reduce... " >&6; } + cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + + #include + #include + +int +main (void) +{ + + std::vector v{7,8,9}; + auto sumsq = std::transform_reduce(v.cbegin(), v.cend(), 0, std::plus{}, [](auto x){return x*x;}); + + ; + return 0; +} +_ACEOF +if ac_fn_cxx_try_compile "$LINENO" +then : + + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: yes" >&5 +printf "%s\n" "yes" >&6; } + have_cxx17_transform_reduce=yes + +printf "%s\n" "#define HAVE_CXX17_TRANSFORM_REDUCE 1" >>confdefs.h + + +else case e in #( + e) + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: no" >&5 +printf "%s\n" "no" >&6; } + ;; +esac +fi +rm -f core conftest.err conftest.$ac_objext conftest.beam conftest.$ac_ext + + CXXFLAGS="$old_CXXFLAGS" + ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + + if test "x$have_cxx17_transform_reduce" != "xyes" +then : + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: WARNING: libMesh prefers C++17 support for std::transform_reduce" >&5 +printf "%s\n" "$as_me: WARNING: libMesh prefers C++17 support for std::transform_reduce" >&2;} +fi + fi if test "$HAVE_TESTED_CXX" = "1" -a "x$have_cxx_all" = xyes then : @@ -21580,6 +21645,71 @@ then : printf "%s\n" "$as_me: WARNING: libMesh prefers C++17 support for set/map merge" >&2;} fi + + have_cxx17_transform_reduce=no + + ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + + old_CXXFLAGS="$CXXFLAGS" + CXXFLAGS="$CXXFLAGS $libmesh_CXXFLAGS" + + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking for C++17 std::transform_reduce" >&5 +printf %s "checking for C++17 std::transform_reduce... " >&6; } + cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + + #include + #include + +int +main (void) +{ + + std::vector v{7,8,9}; + auto sumsq = std::transform_reduce(v.cbegin(), v.cend(), 0, std::plus{}, [](auto x){return x*x;}); + + ; + return 0; +} +_ACEOF +if ac_fn_cxx_try_compile "$LINENO" +then : + + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: yes" >&5 +printf "%s\n" "yes" >&6; } + have_cxx17_transform_reduce=yes + +printf "%s\n" "#define HAVE_CXX17_TRANSFORM_REDUCE 1" >>confdefs.h + + +else case e in #( + e) + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: no" >&5 +printf "%s\n" "no" >&6; } + ;; +esac +fi +rm -f core conftest.err conftest.$ac_objext conftest.beam conftest.$ac_ext + + CXXFLAGS="$old_CXXFLAGS" + ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + + if test "x$have_cxx17_transform_reduce" != "xyes" +then : + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: WARNING: libMesh prefers C++17 support for std::transform_reduce" >&5 +printf "%s\n" "$as_me: WARNING: libMesh prefers C++17 support for std::transform_reduce" >&2;} +fi + fi if test "$HAVE_TESTED_CXX" = "1" -a "x$have_cxx_all" = xyes then : @@ -27559,6 +27689,71 @@ then : printf "%s\n" "$as_me: WARNING: libMesh prefers C++17 support for set/map merge" >&2;} fi + + have_cxx17_transform_reduce=no + + ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + + old_CXXFLAGS="$CXXFLAGS" + CXXFLAGS="$CXXFLAGS $libmesh_CXXFLAGS" + + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking for C++17 std::transform_reduce" >&5 +printf %s "checking for C++17 std::transform_reduce... " >&6; } + cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + + #include + #include + +int +main (void) +{ + + std::vector v{7,8,9}; + auto sumsq = std::transform_reduce(v.cbegin(), v.cend(), 0, std::plus{}, [](auto x){return x*x;}); + + ; + return 0; +} +_ACEOF +if ac_fn_cxx_try_compile "$LINENO" +then : + + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: yes" >&5 +printf "%s\n" "yes" >&6; } + have_cxx17_transform_reduce=yes + +printf "%s\n" "#define HAVE_CXX17_TRANSFORM_REDUCE 1" >>confdefs.h + + +else case e in #( + e) + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: no" >&5 +printf "%s\n" "no" >&6; } + ;; +esac +fi +rm -f core conftest.err conftest.$ac_objext conftest.beam conftest.$ac_ext + + CXXFLAGS="$old_CXXFLAGS" + ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + + if test "x$have_cxx17_transform_reduce" != "xyes" +then : + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: WARNING: libMesh prefers C++17 support for std::transform_reduce" >&5 +printf "%s\n" "$as_me: WARNING: libMesh prefers C++17 support for std::transform_reduce" >&2;} +fi + fi if test "$HAVE_TESTED_CXX" = "1" -a "x$have_cxx_all" = xyes then : @@ -32026,6 +32221,71 @@ then : printf "%s\n" "$as_me: WARNING: libMesh prefers C++17 support for set/map merge" >&2;} fi + + have_cxx17_transform_reduce=no + + ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + + old_CXXFLAGS="$CXXFLAGS" + CXXFLAGS="$CXXFLAGS $libmesh_CXXFLAGS" + + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking for C++17 std::transform_reduce" >&5 +printf %s "checking for C++17 std::transform_reduce... " >&6; } + cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + + #include + #include + +int +main (void) +{ + + std::vector v{7,8,9}; + auto sumsq = std::transform_reduce(v.cbegin(), v.cend(), 0, std::plus{}, [](auto x){return x*x;}); + + ; + return 0; +} +_ACEOF +if ac_fn_cxx_try_compile "$LINENO" +then : + + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: yes" >&5 +printf "%s\n" "yes" >&6; } + have_cxx17_transform_reduce=yes + +printf "%s\n" "#define HAVE_CXX17_TRANSFORM_REDUCE 1" >>confdefs.h + + +else case e in #( + e) + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: no" >&5 +printf "%s\n" "no" >&6; } + ;; +esac +fi +rm -f core conftest.err conftest.$ac_objext conftest.beam conftest.$ac_ext + + CXXFLAGS="$old_CXXFLAGS" + ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + + if test "x$have_cxx17_transform_reduce" != "xyes" +then : + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: WARNING: libMesh prefers C++17 support for std::transform_reduce" >&5 +printf "%s\n" "$as_me: WARNING: libMesh prefers C++17 support for std::transform_reduce" >&2;} +fi + fi if test "$HAVE_TESTED_CXX" = "1" -a "x$have_cxx_all" = xyes then : @@ -36013,6 +36273,71 @@ then : printf "%s\n" "$as_me: WARNING: libMesh prefers C++17 support for set/map merge" >&2;} fi + + have_cxx17_transform_reduce=no + + ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + + old_CXXFLAGS="$CXXFLAGS" + CXXFLAGS="$CXXFLAGS $libmesh_CXXFLAGS" + + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: checking for C++17 std::transform_reduce" >&5 +printf %s "checking for C++17 std::transform_reduce... " >&6; } + cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + + #include + #include + +int +main (void) +{ + + std::vector v{7,8,9}; + auto sumsq = std::transform_reduce(v.cbegin(), v.cend(), 0, std::plus{}, [](auto x){return x*x;}); + + ; + return 0; +} +_ACEOF +if ac_fn_cxx_try_compile "$LINENO" +then : + + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: yes" >&5 +printf "%s\n" "yes" >&6; } + have_cxx17_transform_reduce=yes + +printf "%s\n" "#define HAVE_CXX17_TRANSFORM_REDUCE 1" >>confdefs.h + + +else case e in #( + e) + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: no" >&5 +printf "%s\n" "no" >&6; } + ;; +esac +fi +rm -f core conftest.err conftest.$ac_objext conftest.beam conftest.$ac_ext + + CXXFLAGS="$old_CXXFLAGS" + ac_ext=cpp +ac_cpp='$CXXCPP $CPPFLAGS' +ac_compile='$CXX -c $CXXFLAGS $CPPFLAGS conftest.$ac_ext >&5' +ac_link='$CXX -o conftest$ac_exeext $CXXFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5' +ac_compiler_gnu=$ac_cv_cxx_compiler_gnu + + + if test "x$have_cxx17_transform_reduce" != "xyes" +then : + { printf "%s\n" "$as_me:${as_lineno-$LINENO}: WARNING: libMesh prefers C++17 support for std::transform_reduce" >&5 +printf "%s\n" "$as_me: WARNING: libMesh prefers C++17 support for std::transform_reduce" >&2;} +fi + fi if test "$HAVE_TESTED_CXX" = "1" -a "x$have_cxx_all" = xyes then : diff --git a/include/Makefile.in b/include/Makefile.in index e4979ee3d7..07653a6afd 100644 --- a/include/Makefile.in +++ b/include/Makefile.in @@ -920,6 +920,7 @@ include_HEADERS = \ parallel/libmesh_call_mpi.h \ parallel/parallel.h \ parallel/parallel_algebra.h \ + parallel/parallel_algorithms.h \ parallel/parallel_bin_sorter.h \ parallel/parallel_eigen.h \ parallel/parallel_elem.h \ diff --git a/include/libmesh/Makefile.in b/include/libmesh/Makefile.in index 62bc6c82e1..6b8ae24205 100644 --- a/include/libmesh/Makefile.in +++ b/include/libmesh/Makefile.in @@ -629,13 +629,13 @@ BUILT_SOURCES = dirichlet_boundaries.h dof_map.h dof_map_base.h \ type_tensor.h type_vector.h vector_value.h wrapped_function.h \ wrapped_functor.h wrapped_petsc.h zero_function.h \ libmesh_call_mpi.h parallel.h parallel_algebra.h \ - parallel_bin_sorter.h parallel_conversion_utils.h \ - parallel_eigen.h parallel_elem.h parallel_fe_type.h \ - parallel_ghost_sync.h parallel_hilbert.h parallel_histogram.h \ - parallel_node.h parallel_object.h parallel_only.h \ - parallel_sort.h threads.h threads_allocators.h threads_none.h \ - threads_pthread.h threads_spin_mutex_forward.h threads_tbb.h \ - centroid_partitioner.h hilbert_sfc_partitioner.h \ + parallel_algorithms.h parallel_bin_sorter.h \ + parallel_conversion_utils.h parallel_eigen.h parallel_elem.h \ + parallel_fe_type.h parallel_ghost_sync.h parallel_hilbert.h \ + parallel_histogram.h parallel_node.h parallel_object.h \ + parallel_only.h parallel_sort.h threads.h threads_allocators.h \ + threads_none.h threads_pthread.h threads_spin_mutex_forward.h \ + threads_tbb.h centroid_partitioner.h hilbert_sfc_partitioner.h \ linear_partitioner.h mapped_subdomain_partitioner.h \ metis_csr_graph.h metis_partitioner.h morton_sfc_partitioner.h \ parmetis_helper.h parmetis_partitioner.h partitioner.h \ @@ -1844,6 +1844,9 @@ parallel.h: $(top_srcdir)/include/parallel/parallel.h parallel_algebra.h: $(top_srcdir)/include/parallel/parallel_algebra.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ +parallel_algorithms.h: $(top_srcdir)/include/parallel/parallel_algorithms.h + $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ + parallel_bin_sorter.h: $(top_srcdir)/include/parallel/parallel_bin_sorter.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ diff --git a/include/libmesh_config.h.in b/include/libmesh_config.h.in index a8318aa5e6..ec538f5460 100644 --- a/include/libmesh_config.h.in +++ b/include/libmesh_config.h.in @@ -343,6 +343,9 @@ /* Flag indicating whether compiler supports std::*::merge */ #undef HAVE_CXX17_SPLICING +/* Flag indicating whether compiler supports std::transform_reduce */ +#undef HAVE_CXX17_TRANSFORM_REDUCE + /* define if the compiler supports basic C++20 syntax */ #undef HAVE_CXX20