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 <numeric>
+    #include <vector>
+
+int
+main (void)
+{
+
+    std::vector<int> 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 <numeric>
+    #include <vector>
+
+int
+main (void)
+{
+
+    std::vector<int> 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 <numeric>
+    #include <vector>
+
+int
+main (void)
+{
+
+    std::vector<int> 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 <numeric>
+    #include <vector>
+
+int
+main (void)
+{
+
+    std::vector<int> 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 <numeric>
+    #include <vector>
+
+int
+main (void)
+{
+
+    std::vector<int> 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/base/libmesh_common.h b/include/base/libmesh_common.h
index a95b0c4d3d..40c96202b6 100644
--- a/include/base/libmesh_common.h
+++ b/include/base/libmesh_common.h
@@ -569,8 +569,7 @@ struct casting_compare {
 template<class ...Args> inline void libmesh_ignore( const Args&... ) { }
 
 
-// A workaround for the lack of C++17 merge() support in some
-// compilers
+// Workaround for incomplete C++17 support in some compilers/libs
 
 #ifdef LIBMESH_HAVE_CXX17_SPLICING
 template <typename T>
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/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
 
diff --git a/include/numerics/dense_matrix.h b/include/numerics/dense_matrix.h
index 044221b308..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)
@@ -42,9 +43,9 @@
 #endif
 
 // C++ includes
-#include <vector>
 #include <algorithm>
 #include <initializer_list>
+#include <vector>
 
 #ifdef LIBMESH_HAVE_METAPHYSICL
 #include "metaphysicl/dualnumber_decl.h"
@@ -564,6 +565,62 @@ class DenseMatrix : public DenseMatrixBase<T>
    */
   T det();
 
+  /**
+   * Returns true iff every entry is finite.
+   */
+  friend bool isfinite (const DenseMatrix<T> & 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<T> & 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<T> & 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
@@ -864,8 +921,9 @@ template<typename T>
 inline
 void DenseMatrix<T>::swap(DenseMatrix<T> & 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;
@@ -1084,17 +1142,11 @@ auto DenseMatrix<T>::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;
+  typedef decltype(libmesh_real(T(0))) realfromT;
+  return libmesh_transform_reduce
+    (_val.begin(), _val.end(), std::numeric_limits<realfromT>::max(),
+     [](const auto & a, const auto & b){using std::min; return min(a,b);},
+     [](const T & v){return libmesh_real(v);});
 }
 
 
@@ -1105,17 +1157,11 @@ auto DenseMatrix<T>::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;
+  typedef decltype(libmesh_real(T(0))) realfromT;
+  return libmesh_transform_reduce
+    (_val.begin(), _val.end(), std::numeric_limits<realfromT>::lowest(),
+     [](const auto & a, const auto & b){using std::max; return max(a,b);},
+     [](const T & v){return libmesh_real(v);});
 }
 
 
@@ -1218,70 +1264,6 @@ T DenseMatrix<T>::transpose (const unsigned int i,
 // }
 
 
-// A matrix is finite iff every component is
-template <typename T>
-bool isfinite (const DenseMatrix<T> & 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 <typename T>
-bool isinf (const DenseMatrix<T> & 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 <typename T>
-bool isnan (const DenseMatrix<T> & 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..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
@@ -38,8 +39,9 @@
 #endif
 
 // C++ includes
-#include <vector>
+#include <algorithm>
 #include <initializer_list>
+#include <vector>
 
 namespace libMesh
 {
@@ -294,6 +296,62 @@ class DenseVector : public DenseVectorBase<T>
   typename std::vector<T>::const_iterator end() const { return _val.end(); }
   typename std::vector<T>::iterator end() { return _val.end(); }
 
+  /**
+   * Returns true iff every entry is finite.
+   */
+  friend bool isfinite (const DenseVector<T> & 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<T> & 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<T> & var)
+  {
+    using std::isnan;
+    using libMesh::isnan; // for T==complex
+    for (const T & v : var._val)
+      if (isnan(v))
+        return true;
+    return false;
+  }
+
 private:
 
   /**
@@ -332,13 +390,7 @@ DenseVector<T>::DenseVector (const DenseVector<T2> & other_vector) :
 {
   const std::vector<T2> & other_vals = other_vector.get_values();
 
-  _val.clear();
-
-  const int N = cast_int<int>(other_vals.size());
-  _val.reserve(N);
-
-  for (int i=0; i<N; i++)
-    _val.push_back(other_vals[i]);
+  _val.assign(other_vals.begin(), other_vals.end());
 }
 
 
@@ -368,15 +420,7 @@ inline
 DenseVector<T> & DenseVector<T>::operator = (const DenseVector<T2> & other_vector)
 {
   const std::vector<T2> & other_vals = other_vector.get_values();
-
-  _val.clear();
-
-  const int N = cast_int<int>(other_vals.size());
-  _val.reserve(N);
-
-  for (int i=0; i<N; i++)
-    _val.push_back(other_vals[i]);
-
+  _val.assign(other_vals.begin(), other_vals.end());
   return *this;
 }
 
@@ -452,9 +496,8 @@ template<typename T>
 inline
 void DenseVector<T>::scale (const T factor)
 {
-  const int N = cast_int<int>(_val.size());
-  for (int i=0; i<N; i++)
-    _val[i] *= factor;
+  for (auto & v : _val)
+    v *= factor;
 }
 
 
@@ -606,15 +649,11 @@ inline
 Real DenseVector<T>::min () const
 {
   libmesh_assert (this->size());
-  Real my_min = libmesh_real((*this)(0));
-
-  const int N = cast_int<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(T(0))) realfromT;
+  return libmesh_transform_reduce
+    (_val.begin(), _val.end(), std::numeric_limits<realfromT>::max(),
+     [](const auto & a, const auto & b){using std::min; return min(a,b);},
+     [](const T & v){return libmesh_real(v);});
 }
 
 
@@ -624,15 +663,11 @@ inline
 Real DenseVector<T>::max () const
 {
   libmesh_assert (this->size());
-  Real my_max = libmesh_real((*this)(0));
-
-  const int N = cast_int<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(T(0))) realfromT;
+  return libmesh_transform_reduce
+    (_val.begin(), _val.end(), std::numeric_limits<realfromT>::lowest(),
+     [](const auto & a, const auto & b){using std::max; return max(a,b);},
+     [](const T & v){return libmesh_real(v);});
 }
 
 
@@ -647,12 +682,9 @@ Real DenseVector<T>::l1_norm () const
 #ifdef LIBMESH_HAVE_EIGEN
   return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).template lpNorm<1>();
 #else
-  Real my_norm = 0.;
-  const int N = cast_int<int>(_val.size());
-  for (int i=0; i!=N; i++)
-    my_norm += std::abs((*this)(i));
-
-  return my_norm;
+  return libmesh_transform_reduce
+    (_val.begin(), _val.end(), Real(0), std::plus<>(),
+     [](const T & v){using std::abs; return abs(v);});
 #endif
 }
 
@@ -668,17 +700,10 @@ Real DenseVector<T>::l2_norm () const
 #ifdef LIBMESH_HAVE_EIGEN
   return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).norm();
 #else
-  Real my_norm = 0.;
-  const int N = cast_int<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(libmesh_transform_reduce
+    (_val.begin(), _val.end(), Real(0), std::plus<>(),
+     [](const T & v){return TensorTools::norm_sq(v);}));
 #endif
 }
 
@@ -694,15 +719,10 @@ Real DenseVector<T>::linfty_norm () const
 #ifdef LIBMESH_HAVE_EIGEN
   return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).template lpNorm<Eigen::Infinity>();
 #else
-  Real my_norm = TensorTools::norm_sq((*this)(0));
-
-  const int N = cast_int<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 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);});
 #endif
 }
 
@@ -722,63 +742,6 @@ void DenseVector<T>::get_principal_subvector (unsigned int sub_n,
 }
 
 
-// A vector is finite iff every component is
-template <typename T>
-bool isfinite (const DenseVector<T> & 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 <typename T>
-bool isinf (const DenseVector<T> & 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 <typename T>
-bool isnan (const DenseVector<T> & 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
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 <numeric>
+
+// Workaround incomplete C++17 support in some compilers/libs
+
+#ifdef LIBMESH_HAVE_CXX17_TRANSFORM_REDUCE
+template <typename InputIt, class T, class BinaryOp, class UnaryOp>
+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 <typename InputIt, class T, class BinaryOp, class UnaryOp>
+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
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 <Eigen/Core>
+#include "libmesh/restore_warnings.h"
 
 namespace libMesh
 {
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 <numeric>
+    @%:@include <vector>
+    ]], [[
+    std::vector<int> 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++])
   ])
 
 
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<Real> 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<Real> 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<Real> 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;