diff --git a/src/TiledArray/math/linalg/basic.h b/src/TiledArray/math/linalg/basic.h index 2045c8a82c..e64b24af5a 100644 --- a/src/TiledArray/math/linalg/basic.h +++ b/src/TiledArray/math/linalg/basic.h @@ -155,6 +155,17 @@ inline void vec_multiply( a1.array() *= a2.array(); } +template +inline auto clone(const Eigen::MatrixBase& a) { + return a.eval(); +} + +template +inline auto clone( + const Eigen::Block& a) { + return a.eval(); +} + template inline void scale(Eigen::MatrixBase& a, S scaling_factor) { using numeric_type = typename Eigen::MatrixBase::value_type; @@ -239,6 +250,21 @@ inline auto norm2( return m.template lpNorm<2>(); } +template +inline auto volume(const Eigen::MatrixBase& m) { + return m.size(); +} + +template +inline auto abs_min(const Eigen::MatrixBase& m) { + return m.array().abs().minCoeff(); +} + +template +inline auto abs_max(const Eigen::MatrixBase& m) { + return m.array().abs().maxCoeff(); +} + } // namespace Eigen #ifndef TILEDARRAY_MATH_LINALG_DISPATCH_W_TTG @@ -253,12 +279,12 @@ inline auto norm2( return scalapack::FN; \ return non_distributed::FN; #elif (TILEDARRAY_HAS_TTG && !TILEDARRAY_HAS_SCALAPACK) -#define TILEDARRAY_MATH_LINALG_DISPATCH_W_TTG(FN, MATRIX) \ - TA_MAX_THREADS; \ - if (get_linalg_backend() == LinearAlgebraBackend::TTG || \ - TiledArray::math::linalg::detail::prefer_distributed(MATRIX)) \ - return TiledArray::math::linalg::ttg::FN; \ - if (get_linalg_backend() == LinearAlgebraBackend::ScaLAPACK) \ +#define TILEDARRAY_MATH_LINALG_DISPATCH_W_TTG(FN, MATRIX) \ + TA_MAX_THREADS; \ + if (get_linalg_backend() == LinearAlgebraBackend::TTG || \ + TiledArray::math::linalg::detail::prefer_distributed(MATRIX)) \ + return TiledArray::math::linalg::ttg::FN; \ + if (get_linalg_backend() == LinearAlgebraBackend::ScaLAPACK) \ TA_EXCEPTION("ScaLAPACK linear algebra backend is not available"); \ return non_distributed::FN; #elif !TILEDARRAY_HAS_TTG && TILEDARRAY_HAS_SCALAPACK @@ -271,11 +297,11 @@ inline auto norm2( return scalapack::FN; \ return non_distributed::FN; #else // !TILEDARRAY_HAS_TTG && !TILEDARRAY_HAS_SCALAPACK -#define TILEDARRAY_MATH_LINALG_DISPATCH_W_TTG(FN, MATRIX) \ - TA_MAX_THREADS; \ - if (get_linalg_backend() == LinearAlgebraBackend::TTG) \ - TA_EXCEPTION("TTG linear algebra backend is not available"); \ - if (get_linalg_backend() == LinearAlgebraBackend::ScaLAPACK) \ +#define TILEDARRAY_MATH_LINALG_DISPATCH_W_TTG(FN, MATRIX) \ + TA_MAX_THREADS; \ + if (get_linalg_backend() == LinearAlgebraBackend::TTG) \ + TA_EXCEPTION("TTG linear algebra backend is not available"); \ + if (get_linalg_backend() == LinearAlgebraBackend::ScaLAPACK) \ TA_EXCEPTION("ScaLAPACK linear algebra backend is not available"); \ return non_distributed::FN; #endif // !TILEDARRAY_HAS_TTG && !TILEDARRAY_HAS_SCALAPACK @@ -297,12 +323,12 @@ inline auto norm2( return scalapack::FN; \ return non_distributed::FN; #elif TILEDARRAY_HAS_TTG && !TILEDARRAY_HAS_SCALAPACK -#define TILEDARRAY_MATH_LINALG_DISPATCH_WO_TTG(FN, MATRIX) \ - TA_MAX_THREADS; \ - if (get_linalg_backend() == LinearAlgebraBackend::TTG) \ - TA_EXCEPTION(TILEDARRAY_MATH_LINALG_DISPATCH_WO_TTG_STRINGIFY( \ - FN) " is not provided by the TTG backend"); \ - if (get_linalg_backend() == LinearAlgebraBackend::ScaLAPACK) \ +#define TILEDARRAY_MATH_LINALG_DISPATCH_WO_TTG(FN, MATRIX) \ + TA_MAX_THREADS; \ + if (get_linalg_backend() == LinearAlgebraBackend::TTG) \ + TA_EXCEPTION(TILEDARRAY_MATH_LINALG_DISPATCH_WO_TTG_STRINGIFY( \ + FN) " is not provided by the TTG backend"); \ + if (get_linalg_backend() == LinearAlgebraBackend::ScaLAPACK) \ TA_EXCEPTION("ScaLAPACK linear algebra backend is not available"); \ return non_distributed::FN; #elif !TILEDARRAY_HAS_TTG && TILEDARRAY_HAS_SCALAPACK @@ -315,11 +341,11 @@ inline auto norm2( return scalapack::FN; \ return non_distributed::FN; #else // !TILEDARRAY_HAS_TTG && !TILEDARRAY_HAS_SCALAPACK -#define TILEDARRAY_MATH_LINALG_DISPATCH_WO_TTG(FN, MATRIX) \ - TA_MAX_THREADS; \ - if (get_linalg_backend() == LinearAlgebraBackend::TTG) \ - TA_EXCEPTION("TTG linear algebra backend is not available"); \ - if (get_linalg_backend() == LinearAlgebraBackend::ScaLAPACK) \ +#define TILEDARRAY_MATH_LINALG_DISPATCH_WO_TTG(FN, MATRIX) \ + TA_MAX_THREADS; \ + if (get_linalg_backend() == LinearAlgebraBackend::TTG) \ + TA_EXCEPTION("TTG linear algebra backend is not available"); \ + if (get_linalg_backend() == LinearAlgebraBackend::ScaLAPACK) \ TA_EXCEPTION("ScaLAPACK linear algebra backend is not available"); \ return non_distributed::FN; #endif // !TILEDARRAY_HAS_TTG && !TILEDARRAY_HAS_SCALAPACK diff --git a/src/TiledArray/math/solvers/conjgrad.h b/src/TiledArray/math/solvers/conjgrad.h index cacfd55d63..efabd59af7 100644 --- a/src/TiledArray/math/solvers/conjgrad.h +++ b/src/TiledArray/math/solvers/conjgrad.h @@ -29,6 +29,7 @@ #include #include #include "TiledArray/dist_array.h" +#include "TiledArray/type_traits.h" namespace TiledArray::math { @@ -44,8 +45,8 @@ namespace TiledArray::math { /// stand-alone functions: /// \li std::size_t volume(const D&) (returns the total number of elements) /// \li D clone(const D&) , returns a deep copy -/// \li value_type minabs_value(const D&) -/// \li value_type maxabs_value(const D&) +/// \li value_type abs_min(const D&) +/// \li value_type abs_max(const D&) /// \li void vec_multiply(D& a, const D& b) (element-wise multiply /// of \c a by \c b ) /// \li value_type inner_product(const D& a, const D& b) @@ -60,7 +61,7 @@ namespace TiledArray::math { // clang-format on template struct ConjugateGradientSolver { - typedef typename D::numeric_type value_type; + typedef TiledArray::detail::numeric_t value_type; /// \param a object of type F /// \param b RHS @@ -73,8 +74,8 @@ struct ConjugateGradientSolver { value_type convergence_target = -1.0) { std::size_t n = volume(preconditioner); - const bool use_diis = false; - DIIS diis; + constexpr bool use_diis = false; + std::conditional_t, char> diis{}; // solution vector D XX_i; @@ -120,7 +121,7 @@ struct ConjugateGradientSolver { scale(RR_i, -1.0); axpy(RR_i, 1.0, b); // RR_i = b - a(XX_i) - if (use_diis) diis.extrapolate(XX_i, RR_i, true); + if constexpr (use_diis) diis.extrapolate(XX_i, RR_i, true); // z_0 = D^-1 . r_0 ZZ_i = RR_i; @@ -144,7 +145,7 @@ struct ConjugateGradientSolver { // r_i -= alpha_i Ap_i axpy(RR_i, -alpha_i, APP_i); - if (use_diis) diis.extrapolate(XX_i, RR_i, true); + if constexpr (use_diis) diis.extrapolate(XX_i, RR_i, true); const value_type r_ip1_norm = norm2(RR_i) / rhs_size; if (r_ip1_norm < convergence_target) { diff --git a/src/TiledArray/math/solvers/diis.h b/src/TiledArray/math/solvers/diis.h index 1407ff327e..86f58bfd5a 100644 --- a/src/TiledArray/math/solvers/diis.h +++ b/src/TiledArray/math/solvers/diis.h @@ -29,6 +29,7 @@ #include #include "TiledArray/dist_array.h" #include "TiledArray/external/eigen.h" +#include "TiledArray/type_traits.h" #include #include @@ -82,7 +83,7 @@ namespace TiledArray::math { template class DIIS { public: - typedef typename D::numeric_type value_type; + typedef TiledArray::detail::numeric_t value_type; typedef typename TiledArray::detail::scalar_t scalar_type; typedef Eigen::Matrix diff --git a/src/TiledArray/tile_interface/clone.h b/src/TiledArray/tile_interface/clone.h index 38c9dddf86..9c0221bf92 100644 --- a/src/TiledArray/tile_interface/clone.h +++ b/src/TiledArray/tile_interface/clone.h @@ -36,7 +36,7 @@ namespace TiledArray { /// \tparam Arg The tile argument type /// \param arg The tile argument to be permuted /// \return A (deep) copy of \c arg -template +template ().clone())> inline auto clone(const Arg& arg) { return arg.clone(); } diff --git a/src/TiledArray/tile_op/tile_interface.h b/src/TiledArray/tile_op/tile_interface.h index 6ab18a2384..5bf8219915 100644 --- a/src/TiledArray/tile_op/tile_interface.h +++ b/src/TiledArray/tile_op/tile_interface.h @@ -969,7 +969,8 @@ inline auto min(const Arg& arg) { /// \tparam Arg The tile argument type /// \param arg The argument to find the maximum /// \return A scalar that is equal to abs(max(arg)) -template +template ().abs_max())> inline auto abs_max(const Arg& arg) { return arg.abs_max(); } @@ -979,7 +980,8 @@ inline auto abs_max(const Arg& arg) { /// \tparam Arg The tile argument type /// \param arg The argument to find the minimum /// \return A scalar that is equal to abs(min(arg)) -template +template ().abs_min())> inline auto abs_min(const Arg& arg) { return arg.abs_min(); } @@ -991,7 +993,9 @@ inline auto abs_min(const Arg& arg) { /// \param left The left-hand argument tile /// \param right The right-hand argument tile /// \return A scalar that is equal to sum_i left[i] * right[i] -template +template ().dot( + std::declval()))> inline auto dot(const Left& left, const Right& right) { return left.dot(right); } @@ -1003,7 +1007,9 @@ inline auto dot(const Left& left, const Right& right) { /// \param left The left-hand argument tile /// \param right The right-hand argument tile /// \return A scalar that is equal to sum_i conj(left[i]) * right[i] -template +template ().inner_product( + std::declval()))> inline auto inner_product(const Left& left, const Right& right) { return left.inner_product(right); } diff --git a/tests/solvers.cpp b/tests/solvers.cpp index 15335880a0..ea8c3ce5e5 100644 --- a/tests/solvers.cpp +++ b/tests/solvers.cpp @@ -167,6 +167,49 @@ struct validate> { } }; +// Eigen specializations + +template <> +struct make_Ax { + using T = Eigen::VectorXd; + + struct Ax { + Ax() : A_(3, 3) { A_ << 1, 2, 3, 2, 5, 8, 3, 8, 15; } + void operator()(const T& x, T& result) const { result = A_ * x; } + Eigen::MatrixXd A_; + }; + Ax operator()() const { return Ax{}; } +}; + +template <> +struct make_b { + using T = Eigen::VectorXd; + + T operator()() const { + T result(3); + result << 1, 4, 0; + return result; + } +}; + +template <> +struct make_pc { + using T = Eigen::VectorXd; + + T operator()() const { return T::Ones(3); } +}; + +template <> +struct validate { + using T = Eigen::VectorXd; + + bool operator()(const T& x) const { + T ref(3); + ref << -6.5, 9., -3.5; + return (x - ref).norm() < 1e-11; + } +}; + BOOST_AUTO_TEST_SUITE(solvers) BOOST_AUTO_TEST_CASE_TEMPLATE(conjugate_gradient, Array, array_types) { @@ -178,4 +221,14 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(conjugate_gradient, Array, array_types) { BOOST_CHECK(validate{}(x)); } +BOOST_AUTO_TEST_CASE(conjugate_gradient_eigen) { + using T = Eigen::VectorXd; + auto Ax = make_Ax{}(); + auto b = make_b{}(); + auto pc = make_pc{}(); + T x; + ConjugateGradientSolver{}(Ax, b, x, pc, 1e-11); + BOOST_CHECK(validate{}(x)); +} + BOOST_AUTO_TEST_SUITE_END()