Skip to content

Update MetaPhysicL #4300

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Nov 6, 2025
+102 −94
Merged

Update MetaPhysicL #4300

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
48c9f67
Update headers for compatibility with kokkos metaphysicl
lindsayad Oct 14, 2025
c9b7ad7
Include raw_type instead of numberarray bc it doesn't depend on config
lindsayad Oct 30, 2025
b843ea4
Run bootstrap
lindsayad Oct 30, 2025
91779e7
Update fparser
lindsayad Nov 3, 2025
d899143
Update MetaPhysicL
lindsayad Oct 29, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions configure
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63450,13 +63450,13 @@ ac_compiler_gnu=$ac_cv_cxx_compiler_gnu
cat confdefs.h - <<_ACEOF >conftest.$ac_ext
/* end confdefs.h. */

#include "metaphysicl/numberarray.h"
#include "metaphysicl/raw_type.h"

int
main (void)
{

MetaPhysicL::NumberArray<4, double> x;
MetaPhysicL::RawType<double> x;

;
return 0;
Expand Down
98 changes: 49 additions & 49 deletions contrib/fparser/fparser_ad.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -800,77 +800,77 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
case cDiv:
--sp; ccout << "s[" << sp << "] /= s[" << (sp+1) << "];\n"; break;
case cMod:
--sp; ccout << "s[" << sp << "] = std::fmod(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
--sp; ccout << "using std::fmod; s[" << sp << "] = fmod(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
case cRDiv:
--sp; ccout << "s[" << sp << "] = s[" << (sp+1) << "] / s[" << sp << "];\n"; break;

case cSin:
ccout << "s[" << sp << "] = std::sin(s[" << sp << "]);\n"; break;
ccout << "using std::sin; s[" << sp << "] = sin(s[" << sp << "]);\n"; break;
case cCos:
ccout << "s[" << sp << "] = std::cos(s[" << sp << "]);\n"; break;
ccout << "using std::cos; s[" << sp << "] = cos(s[" << sp << "]);\n"; break;
case cTan:
ccout << "s[" << sp << "] = std::tan(s[" << sp << "]);\n"; break;
ccout << "using std::tan; s[" << sp << "] = tan(s[" << sp << "]);\n"; break;
case cSinh:
ccout << "s[" << sp << "] = std::sinh(s[" << sp << "]);\n"; break;
ccout << "using std::sinh; s[" << sp << "] = sinh(s[" << sp << "]);\n"; break;
case cCosh:
ccout << "s[" << sp << "] = std::cosh(s[" << sp << "]);\n"; break;
ccout << "using std::cosh; s[" << sp << "] = cosh(s[" << sp << "]);\n"; break;
case cTanh:
ccout << "s[" << sp << "] = std::tanh(s[" << sp << "]);\n"; break;
ccout << "using std::tanh; s[" << sp << "] = tanh(s[" << sp << "]);\n"; break;
// TODO: div by zero -> this->mData->mEvalErrorType=1; return Value_t(0);
case cCsc:
ccout << "s[" << sp << "] = 1.0/std::sin(s[" << sp << "]);\n"; break;
ccout << "using std::sin; s[" << sp << "] = 1.0/sin(s[" << sp << "]);\n"; break;
case cSec:
ccout << "s[" << sp << "] = 1.0/std::cos(s[" << sp << "]);\n"; break;
ccout << "using std::cos; s[" << sp << "] = 1.0/cos(s[" << sp << "]);\n"; break;
case cCot:
ccout << "s[" << sp << "] = 1.0/std::tan(s[" << sp << "]);\n"; break;
ccout << "using std::tan; s[" << sp << "] = 1.0/tan(s[" << sp << "]);\n"; break;
case cSinCos:
ccout << "s[" << (sp+1) << "] = std::cos(s[" << sp << "]);\n";
ccout << "s[" << sp << "] = std::sin(s[" << sp << "]);\n";
ccout << "using std::cos; s[" << (sp+1) << "] = cos(s[" << sp << "]);\n";
ccout << "using std::sin; s[" << sp << "] = sin(s[" << sp << "]);\n";
++sp;
break;
case cSinhCosh:
ccout << "s[" << (sp+1) << "] = std::cosh(s[" << sp << "]);\n";
ccout << "s[" << sp << "] = std::sinh(s[" << sp << "]);\n";
ccout << "using std::cosh; s[" << (sp+1) << "] = cosh(s[" << sp << "]);\n";
ccout << "using std::sinh; s[" << sp << "] = sinh(s[" << sp << "]);\n";
++sp;
break;
case cAsin:
ccout << "s[" << sp << "] = std::asin(s[" << sp << "]);\n"; break;
ccout << " using std::asin; s[" << sp << "] = asin(s[" << sp << "]);\n"; break;
case cAcos:
ccout << "s[" << sp << "] = std::acos(s[" << sp << "]);\n"; break;
ccout << " using std::acos; s[" << sp << "] = acos(s[" << sp << "]);\n"; break;
case cAsinh:
ccout << "s[" << sp << "] = std::asinh(s[" << sp << "]);\n"; break;
ccout << " using std::asinh; s[" << sp << "] = asinh(s[" << sp << "]);\n"; break;
case cAcosh:
ccout << "s[" << sp << "] = std::acosh(s[" << sp << "]);\n"; break;
ccout << " using std::acosh; s[" << sp << "] = acosh(s[" << sp << "]);\n"; break;
case cAtan:
ccout << "s[" << sp << "] = std::atan(s[" << sp << "]);\n"; break;
ccout << " using std::atan; s[" << sp << "] = atan(s[" << sp << "]);\n"; break;
case cAtanh:
ccout << "s[" << sp << "] = std::atanh(s[" << sp << "]);\n"; break;
ccout << " using std::atanh; s[" << sp << "] = atanh(s[" << sp << "]);\n"; break;
case cAtan2:
--sp; ccout << "s[" << sp << "] = std::atan2(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
--sp; ccout << " using std::atan2; s[" << sp << "] = atan2(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
case cHypot:
--sp; ccout << "s[" << sp << "] = std::sqrt(s[" << sp << "]*s[" << sp << "] + s[" << (sp+1) << "]*s[" << (sp+1) << "]);\n"; break;
--sp; ccout << " using std::sqrt; s[" << sp << "] = sqrt(s[" << sp << "]*s[" << sp << "] + s[" << (sp+1) << "]*s[" << (sp+1) << "]);\n"; break;

case cAbs:
ccout << "s[" << sp << "] = std::abs(s[" << sp << "]);\n"; break;
ccout << "using std::abs; s[" << sp << "] = abs(s[" << sp << "]);\n"; break;
case cMax:
--sp; ccout << "s[" << sp << "] = s[" << sp << "] > s[" << (sp+1) << "] ? s[" << sp << "] : s[" << (sp+1) << "];\n"; break;
case cMin:
--sp; ccout << "s[" << sp << "] = s[" << sp << "] < s[" << (sp+1) << "] ? s[" << sp << "] : s[" << (sp+1) << "];\n"; break;
case cTrunc:
ccout << "s[" << sp << "] = s[" << sp << "] < 0 ? std::ceil(s[" << sp << "]) : std::floor(s[" << sp << "]);\n"; break;
ccout << "using std::ceil, std::floor; s[" << sp << "] = s[" << sp << "] < 0 ? ceil(s[" << sp << "]) : floor(s[" << sp << "]);\n"; break;
case cCeil:
ccout << "s[" << sp << "] = std::ceil(s[" << sp << "]);\n"; break;
ccout << "using std::ceil; s[" << sp << "] = ceil(s[" << sp << "]);\n"; break;
case cFloor:
ccout << "s[" << sp << "] = std::floor(s[" << sp << "]);\n"; break;
ccout << "using std::floor; s[" << sp << "] = floor(s[" << sp << "]);\n"; break;
case cInt:
ccout << "s[" << sp << "] = s[" << sp << "] < 0 ? std::ceil(s[" << sp << "] - 0.5) : std::floor(s[" << sp << "] + 0.5);\n"; break;
ccout << "using std::ceil, std::floor; s[" << sp << "] = s[" << sp << "] < 0 ? ceil(s[" << sp << "] - 0.5) : floor(s[" << sp << "] + 0.5);\n"; break;

case cEqual:
//--sp; ccout << "s[" << sp << "] = s[" << sp << "] == s[" << (sp+1) << "];\n"; break;
--sp; ccout << "s[" << sp << "] = std::abs(s[" << sp << "] - s[" << (sp+1) << "]) <= eps;\n"; break;
--sp; ccout << "using std::abs; s[" << sp << "] = abs(s[" << sp << "] - s[" << (sp+1) << "]) <= eps;\n"; break;
case cNEqual:
//--sp; ccout << "s[" << sp << "] = s[" << sp << "] != s[" << (sp+1) << "];\n"; break;
--sp; ccout << "s[" << sp << "] = std::abs(s[" << sp << "] - s[" << (sp+1) << "]) > eps;\n"; break;
--sp; ccout << "using std::abs; s[" << sp << "] = abs(s[" << sp << "] - s[" << (sp+1) << "]) > eps;\n"; break;
case cLess:
--sp; ccout << "s[" << sp << "] = s[" << sp << "] < (s[" << (sp+1) << "] - eps);\n"; break;
case cLessOrEq:
Expand All @@ -880,33 +880,33 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
case cGreaterOrEq:
--sp; ccout << "s[" << sp << "] = (s[" << sp << "] + eps) >= s[" << (sp+1) << "];\n"; break;
case cNot:
ccout << "s[" << sp << "] = std::abs(s[" << sp << "]) < 0.5;\n"; break;
ccout << "using std::abs; s[" << sp << "] = abs(s[" << sp << "]) < 0.5;\n"; break;
case cNotNot:
ccout << "s[" << sp << "] = std::abs(s[" << sp << "]) >= 0.5;\n"; break;
ccout << "using std::abs; s[" << sp << "] = abs(s[" << sp << "]) >= 0.5;\n"; break;
case cAbsNot:
ccout << "s[" << sp << "] = s[" << sp << "] < 0.5;\n"; break;
case cAbsNotNot:
ccout << "s[" << sp << "] = s[" << sp << "] >= 0.5;\n"; break;
case cOr:
--sp; ccout << "s[" << sp << "] = (std::abs(s[" << sp << "]) >= 0.5) || (std::abs(s[" << (sp+1) << "]) >= 0.5);\n"; break;
--sp; ccout << "using std::abs; s[" << sp << "] = (abs(s[" << sp << "]) >= 0.5) || (abs(s[" << (sp+1) << "]) >= 0.5);\n"; break;
case cAbsOr:
--sp; ccout << "s[" << sp << "] = (s[" << sp << "] >= 0.5) || (s[" << (sp+1) << "] >= 0.5);\n"; break;
case cAnd:
--sp; ccout << "s[" << sp << "] = (std::abs(s[" << sp << "]) >= 0.5) && (std::abs(s[" << (sp+1) << "]) >= 0.5);\n"; break;
--sp; ccout << "using std::abs; s[" << sp << "] = (abs(s[" << sp << "]) >= 0.5) && (abs(s[" << (sp+1) << "]) >= 0.5);\n"; break;
case cAbsAnd:
--sp; ccout << "s[" << sp << "] = (s[" << sp << "] >= 0.5) && (s[" << (sp+1) << "] >= 0.5);\n"; break;

case cLog:
ccout << "s[" << sp << "] = std::log(s[" << sp << "]);\n"; break;
ccout << "using std::log; s[" << sp << "] = log(s[" << sp << "]);\n"; break;
case cLog2:
#ifdef FP_SUPPORT_CPLUSPLUS11_MATH_FUNCS
ccout << "s[" << (sp-1) << "] = std::log2(s[" << (sp-1) << "]);\n";
ccout << "using std::log2; s[" << (sp-1) << "] = log2(s[" << (sp-1) << "]);\n";
#else
ccout << "s[" << sp << "] = std::log(s[" << sp << "])/log(2.0);\n";
ccout << "using std::log; s[" << sp << "] = log(s[" << sp << "])/log(2.0);\n";
#endif
break;
case cLog10:
ccout << "s[" << sp << "] = std::log10(s[" << sp << "]);\n"; break;
ccout << "using std::log10; s[" << sp << "] = log10(s[" << sp << "]);\n"; break;

case cNeg:
ccout << "s[" << sp << "] = -s[" << sp << "];\n"; break;
Expand All @@ -929,14 +929,14 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
{
// --sp; ccout << "s[" << sp << "] = s[" << sp << "] < s[" << (sp+1) << "] ? std::log(s[" << (sp+1) << "]) + (s[" << sp << "] - s[" << (sp+1) << "]) / s[" << (sp+1) << "] : std::log(s[" << sp << "]);\n";
// --sp; ccout << "s[" << sp << "] = s[" << sp << "] < s[" << (sp+1) << "] ? std::log(s[" << (sp+1) << "]) - 1.5 + 2.0/s[" << (sp+1) << "] * s[" << sp << "] - 0.5/(s[" << (sp+1) << "]*s[" << (sp+1) << "]) * s[" << sp << "]*s[" << sp << "] : std::log(s[" << sp << "]);\n";
--sp; ccout << "s[" << sp << "] = s[" << sp << "] < s[" << (sp+1) << "] ? std::log(s[" << (sp+1) << "]) + (s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "] - std::pow((s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "],2.0)/2.0 + std::pow((s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "],3.0)/3.0 : std::log(s[" << sp << "]);\n";
--sp; ccout << "using std::log, std::pow; s[" << sp << "] = s[" << sp << "] < s[" << (sp+1) << "] ? log(s[" << (sp+1) << "]) + (s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "] - pow((s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "],2.0)/2.0 + pow((s[" << sp << "]-s[" << (sp+1) << "])/s[" << (sp+1) << "],3.0)/3.0 : log(s[" << sp << "]);\n";
}
else if (function == mFErf)
{
#if LIBMESH_HAVE_CXX11_ERF
ccout << "s[" << sp << "] = std::erf(s[" << sp << "]);\n";
ccout << "using std::erf; s[" << sp << "] = erf(s[" << sp << "]);\n";
#else
std::cerr << "Libmesh is not compiled with c++11 so std::erf is not supported by JIT.\n";
std::cerr << "Libmesh is not compiled with c++11 so erf is not supported by JIT.\n";
return false;
#endif
}
Expand All @@ -958,28 +958,28 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
break;
}
case cLog2by:
--sp; ccout << "s[" << sp << "] = std::log(s[" << sp << "])/log(2.0) * s[" << (sp+1) << "];\n"; break;
--sp; ccout << "using std::log; s[" << sp << "] = log(s[" << sp << "])/log(2.0) * s[" << (sp+1) << "];\n"; break;
case cNop:
break;
#endif

case cSqr:
ccout << "s[" << sp << "] *= s[" << sp << "];\n"; break;
case cSqrt:
ccout << "s[" << sp << "] = std::sqrt(s[" << sp << "]);\n"; break;
ccout << "using std::sqrt; s[" << sp << "] = sqrt(s[" << sp << "]);\n"; break;
case cRSqrt:
ccout << "s[" << sp << "] = std::pow(s[" << sp << "], (-0.5));\n"; break;
ccout << "using std::pow; s[" << sp << "] = pow(s[" << sp << "], (-0.5));\n"; break;
case cPow:
--sp; ccout << "s[" << sp << "] = std::pow(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
--sp; ccout << "using std::pow; s[" << sp << "] = pow(s[" << sp << "], s[" << (sp+1) << "]);\n"; break;
case cExp:
ccout << "s[" << sp << "] = std::exp(s[" << sp << "]);\n"; break;
ccout << "using std::exp; s[" << sp << "] = exp(s[" << sp << "]);\n"; break;
case cExp2:
ccout << "s[" << sp << "] = std::pow(2.0, s[" << sp << "]);\n"; break;
ccout << "using std::pow; s[" << sp << "] = pow(2.0, s[" << sp << "]);\n"; break;
case cCbrt:
#ifdef FP_SUPPORT_CPLUSPLUS11_MATH_FUNCS
ccout << "s[" << sp << "] = std::cbrt(s[" << sp << "]);\n"; break;
ccout << "using std::cbrt; s[" << sp << "] = cbrt(s[" << sp << "]);\n"; break;
#else
ccout << "s[" << sp << "] = s[" << sp << "] == 0 ? 0 : (s[" << sp << "] > 0 ? std::exp(std::log(s[" << sp << "])/3.0) : -std::exp(std::log(-s[" << sp << "])/3.0));\n"; break;
ccout << "using std::exp, std::log; s[" << sp << "] = s[" << sp << "] == 0 ? 0 : (s[" << sp << "] > 0 ? exp(log(s[" << sp << "])/3.0) : -exp(log(-s[" << sp << "])/3.0));\n"; break;
#endif

case cJump:
Expand All @@ -991,7 +991,7 @@ bool FunctionParserADBase<Value_t>::JITCodeGen(std::ostream & ccout, const std::
if (op == cIf)
ccout << "if (s[" << sp-- << "] < 0.5) ";
if (op == cAbsIf)
ccout << "if (std::abs(s[" << sp-- << "]) < 0.5) ";
ccout << "using std::abs; if (abs(s[" << sp-- << "]) < 0.5) ";

if (ip >= ByteCode.size())
ccout << "*ret = s[" << sp << "]; return;\n";
Expand Down
2 changes: 1 addition & 1 deletion contrib/metaphysicl
Open in desktop
Submodule metaphysicl updated 81 files
+13 −4 Makefile.in
+20 −0 NEWS
+1 −0 aclocal.m4
+557 −19 configure
+29 −1 configure.ac
+13 −4 doxygen/Makefile.in
+14 −0 m4/config_summary.m4
+135 −0 m4/kokkos.m4
+8 −2 src/Makefile.am
+21 −6 src/Makefile.in
+41 −0 src/core/include/metaphysicl/metaphysicl_device.h
+2 −2 src/numerics/include/metaphysicl/dualderivatives.h
+152 −192 src/numerics/include/metaphysicl/dualnumber.h
+132 −152 src/numerics/include/metaphysicl/dualnumber_decl.h
+20 −20 src/numerics/include/metaphysicl/dualsemidynamicsparsenumberarray.h
+36 −33 src/numerics/include/metaphysicl/dualsemidynamicsparsenumberarray_decl.h
+4 −2 src/numerics/include/metaphysicl/dualshadow.h
+4 −1 src/numerics/include/metaphysicl/dualshadowdynamicsparsearray.h
+4 −1 src/numerics/include/metaphysicl/dualshadowdynamicsparsevector.h
+4 −2 src/numerics/include/metaphysicl/dualshadowsparsestruct.h
+4 −2 src/numerics/include/metaphysicl/dualshadowsparsevector.h
+2 −1 src/numerics/include/metaphysicl/dualshadowvector.h
+3 −1 src/numerics/include/metaphysicl/dualsparsenumberstruct.h
+232 −198 src/numerics/include/metaphysicl/dynamicsparsenumberbase.h
+84 −112 src/numerics/include/metaphysicl/dynamicsparsenumberbase_decl.h
+39 −71 src/numerics/include/metaphysicl/numberarray.h
+22 −19 src/numerics/include/metaphysicl/parallel_semidynamicsparsenumberarray.h
+40 −38 src/numerics/include/metaphysicl/raw_type.h
+34 −32 src/numerics/include/metaphysicl/semidynamicsparsenumberarray.h
+78 −48 src/numerics/include/metaphysicl/semidynamicsparsenumberarray_decl.h
+132 −186 src/numerics/include/metaphysicl/shadownumber.h
+314 −0 src/numerics/include/metaphysicl/shadownumber_decl.h
+40 −74 src/numerics/include/metaphysicl/sparsenumberarray.h
+29 −64 src/numerics/include/metaphysicl/sparsenumberstruct.h
+1 −1 src/numerics/include/metaphysicl/sparsenumberutils.h
+37 −27 src/numerics/include/metaphysicl/sparsenumbervector.h
+155 −0 src/utilities/include/metaphysicl/dynamic_array_wrapper.h
+40 −0 src/utilities/include/metaphysicl/dynamic_kokkos_array_wrapper.h
+4 −142 src/utilities/include/metaphysicl/dynamic_std_array_wrapper.h
+15 −14 src/utilities/include/metaphysicl/metaphysicl_asserts.h
+3 −0 src/utilities/include/metaphysicl/metaphysicl_config.h.tmp.in
+3 −1 src/utilities/include/metaphysicl/metaphysicl_exceptions.h
+153 −0 src/utilities/include/metaphysicl/metaphysicl_math.h
+78 −0 src/utilities/include/metaphysicl/metaphysicl_numeric_limits.h
+8 −8 src/utilities/include/metaphysicl/parallel_dynamic_array_wrapper.h
+9 −7 src/utilities/include/metaphysicl/testable.h
+25 −0 test/Makefile.am
+82 −34 test/Makefile.in
+18 −17 test/complex_derivs_unit.C
+5 −36 test/derivs_unit.C
+1 −1 test/divgrad_unit.C
+1 −1 test/dualnamedarray_unit.C
+1 −1 test/dynamic_sparse_vector_navier_unit.C
+1 −1 test/dynamic_sparse_vector_pde_unit.C
+2 −11 test/dynamic_std_array_wrapper_unit.C
+1 −1 test/identities_unit.C
+156 −0 test/kokkos_sparse_derivs_unit.K
+2 −1 test/math_structs.h
+2 −2 test/namedindexarray_unit.C
+21 −19 test/navier_unit.h
+1 −1 test/parallel_unit.C
+1 −1 test/pde_unit.C
+21 −19 test/pde_unit.h
+3 −0 test/run_unit_tests.sh.in
+1 −1 test/shadow_dynamic_sparse_vector_navier_unit.C
+1 −1 test/shadow_dynamic_sparse_vector_pde_unit.C
+1 −1 test/shadow_sparse_struct_navier_unit.C
+1 −1 test/shadow_sparse_struct_pde_unit.C
+1 −1 test/shadow_sparse_vector_navier_unit.C
+1 −1 test/shadow_sparse_vector_pde_unit.C
+1 −1 test/shadow_vector_navier_unit.C
+1 −1 test/shadow_vector_pde_unit.C
+25 −13 test/sparse_derivs_unit.C
+1 −1 test/sparse_identities_unit.C
+1 −1 test/sparse_struct_navier_unit.C
+1 −1 test/sparse_struct_pde_unit.C
+1 −1 test/sparse_vector_navier_unit.C
+1 −1 test/sparse_vector_pde_unit.C
+19 −17 test/testopt_unit.C
+1 −1 test/vector_navier_unit.C
+1 −1 test/vector_pde_unit.C
22 changes: 12 additions & 10 deletions include/numerics/dense_matrix.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -324,7 +324,7 @@ class DenseMatrix : public DenseMatrixBase<T>
* This is the natural matrix norm that is compatible to the l1-norm
* for vectors, i.e. \f$ |Mv|_1 \leq |M|_1 |v|_1 \f$.
*/
auto l1_norm () const -> decltype(std::abs(T(0)));
auto l1_norm () const;

/**
* \returns The linfty-norm of the matrix, that is, the max row sum:
Expand All @@ -334,7 +334,7 @@ class DenseMatrix : public DenseMatrixBase<T>
* This is the natural matrix norm that is compatible to the
* linfty-norm of vectors, i.e. \f$ |Mv|_\infty \leq |M|_\infty |v|_\infty \f$.
*/
auto linfty_norm () const -> decltype(std::abs(T(0)));
auto linfty_norm () const;

/**
* Left multiplies by the transpose of the matrix \p A.
Expand Down Expand Up @@ -1122,23 +1122,24 @@ auto DenseMatrix<T>::max () const -> decltype(libmesh_real(T(0)))

template<typename T>
inline
auto DenseMatrix<T>::l1_norm () const -> decltype(std::abs(T(0)))
auto DenseMatrix<T>::l1_norm () const
{
libmesh_assert (this->_m);
libmesh_assert (this->_n);

auto columnsum = std::abs(T(0));
using std::abs;
auto columnsum = abs(T(0));
for (unsigned int i=0; i!=this->_m; i++)
{
columnsum += std::abs((*this)(i,0));
columnsum += abs((*this)(i,0));
}
auto my_max = columnsum;
for (unsigned int j=1; j!=this->_n; j++)
{
columnsum = 0.;
for (unsigned int i=0; i!=this->_m; i++)
{
columnsum += std::abs((*this)(i,j));
columnsum += abs((*this)(i,j));
}
my_max = (my_max > columnsum? my_max : columnsum);
}
Expand All @@ -1149,23 +1150,24 @@ auto DenseMatrix<T>::l1_norm () const -> decltype(std::abs(T(0)))

template<typename T>
inline
auto DenseMatrix<T>::linfty_norm () const -> decltype(std::abs(T(0)))
auto DenseMatrix<T>::linfty_norm () const
{
libmesh_assert (this->_m);
libmesh_assert (this->_n);
using std::abs;

auto rowsum = std::abs(T(0));
auto rowsum = abs(T(0));
for (unsigned int j=0; j!=this->_n; j++)
{
rowsum += std::abs((*this)(0,j));
rowsum += abs((*this)(0,j));
}
auto my_max = rowsum;
for (unsigned int i=1; i!=this->_m; i++)
{
rowsum = 0.;
for (unsigned int j=0; j!=this->_n; j++)
{
rowsum += std::abs((*this)(i,j));
rowsum += abs((*this)(i,j));
}
my_max = (my_max > rowsum? my_max : rowsum);
}
Expand Down
12 changes: 8 additions & 4 deletions include/numerics/dense_matrix_impl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -630,6 +630,8 @@ void DenseMatrix<T>::_lu_back_substitute (const DenseVector<T> & b,
template<typename T>
void DenseMatrix<T>::_lu_decompose ()
{
using std::abs;

// If this function was called, there better not be any
// previous decomposition of the matrix.
libmesh_assert_equal_to (this->_decomposition_type, NONE);
Expand All @@ -648,11 +650,11 @@ void DenseMatrix<T>::_lu_decompose ()
// Find the pivot row by searching down the i'th column
_pivots[i] = i;

// std::abs(complex) must return a Real!
auto the_max = std::abs( A(i,i) );
// abs(complex) must return a Real!
auto the_max = abs( A(i,i) );
for (unsigned int j=i+1; j<n_rows; ++j)
{
auto candidate_max = std::abs( A(j,i) );
auto candidate_max = abs( A(j,i) );
if (the_max < candidate_max)
{
the_max = candidate_max;
Expand Down Expand Up @@ -872,6 +874,8 @@ void DenseMatrix<T>::cholesky_solve (const DenseVector<T2> & b,
template<typename T>
void DenseMatrix<T>::_cholesky_decompose ()
{
using std::sqrt;

// If we called this function, there better not be any
// previous decomposition of the matrix.
libmesh_assert_equal_to (this->_decomposition_type, NONE);
Expand Down Expand Up @@ -901,7 +905,7 @@ void DenseMatrix<T>::_cholesky_decompose ()
"Error! Can only use Cholesky decomposition with symmetric positive definite matrices.");
#endif

A(i,i) = std::sqrt(A(i,j));
A(i,i) = sqrt(A(i,j));
}
else
A(j,i) = A(i,j) / A(i,i);
Expand Down
Loading