From 6e9a98b3254977d81091101e6241d3b2bb6fe9e1 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Mon, 5 Apr 2021 13:59:25 -0400 Subject: [PATCH 1/3] solver adaptors for Eigen data --- src/TiledArray/math/linalg/basic.h | 50 ++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/src/TiledArray/math/linalg/basic.h b/src/TiledArray/math/linalg/basic.h index 32d7cbb94d..88f995f8ea 100644 --- a/src/TiledArray/math/linalg/basic.h +++ b/src/TiledArray/math/linalg/basic.h @@ -27,9 +27,12 @@ #define TILEDARRAY_MATH_LINALG_BASIC_H__INCLUDED #include "TiledArray/dist_array.h" +#include "TiledArray/external/eigen.h" namespace TiledArray::math::linalg { +// freestanding adaptors for DistArray needed by solvers like DIIS + template inline void vec_multiply(DistArray& a1, const DistArray& a2) { @@ -59,4 +62,51 @@ inline void axpy(DistArray& y, S alpha, } // namespace TiledArray::math::linalg +namespace Eigen { + +// freestanding adaptors for Eigen::MatrixBase needed by solvers like DIIS + +template +inline void vec_multiply(Eigen::MatrixBase& a1, + const Eigen::MatrixBase& a2) { + a1.array() *= a2.array(); +} + +template +inline void scale(Eigen::MatrixBase& a, S scaling_factor) { + using numeric_type = typename Eigen::MatrixBase::value_type; + a.array() *= numeric_type(scaling_factor); +} + +template +inline void zero(Eigen::MatrixBase& a) { + a = Derived::Zero(a.rows(), a.cols()); +} + +template +inline void axpy(Eigen::MatrixBase& y, S alpha, + const Eigen::MatrixBase& x) { + using numeric_type = typename Eigen::MatrixBase::value_type; + y.array() += numeric_type(alpha) * x.array(); +} + +template +inline auto dot(const Eigen::MatrixBase& l, + const Eigen::MatrixBase& r) { + return l.adjoint().dot(r); +} + +template +inline auto inner_product(const Eigen::MatrixBase& l, + const Eigen::MatrixBase& r) { + return l.dot(r); +} + +template +inline auto norm2(const Eigen::MatrixBase& m) { + return m.template lpNorm<2>(); +} + +} // namespace Eigen + #endif // TILEDARRAY_MATH_LINALG_BASIC_H__INCLUDED From 7a07752f3e6e3874b7321845ff78f414beaa912a Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Mon, 5 Apr 2021 14:00:31 -0400 Subject: [PATCH 2/3] extra sfinae sprinkles --- src/TiledArray/tensor/complex.h | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/TiledArray/tensor/complex.h b/src/TiledArray/tensor/complex.h index 0aba127c9f..66854b3f43 100644 --- a/src/TiledArray/tensor/complex.h +++ b/src/TiledArray/tensor/complex.h @@ -40,7 +40,8 @@ namespace detail { /// \param r The real scalar /// \return `r` template ::value>::type* = nullptr> + typename std::enable_if && + !is_complex::value>::type* = nullptr> TILEDARRAY_FORCE_INLINE R conj(const R r) { return r; } @@ -61,7 +62,8 @@ TILEDARRAY_FORCE_INLINE std::complex conj(const std::complex z) { /// \tparam R A numeric type /// \return `r` template ::value>::type* = nullptr> + typename std::enable_if && + !is_complex::value>::type* = nullptr> TILEDARRAY_FORCE_INLINE auto inner_product(const L l, const R r) { return l * r; } @@ -72,7 +74,8 @@ TILEDARRAY_FORCE_INLINE auto inner_product(const L l, const R r) { /// \tparam R A numeric type /// \return `r` template ::value>::type* = nullptr> + typename std::enable_if && + is_complex::value>::type* = nullptr> TILEDARRAY_FORCE_INLINE auto inner_product(const L l, const R r) { return TiledArray::detail::conj(l) * r; } @@ -85,7 +88,8 @@ TILEDARRAY_FORCE_INLINE auto inner_product(const L l, const R r) { /// \param r The real scalar /// \return `r` template ::value>::type* = nullptr> + typename std::enable_if && + !is_complex::value>::type* = nullptr> TILEDARRAY_FORCE_INLINE R norm(const R r) { return r * r; } From 15658a2da5b5e025595401c2c9ca390fe708ad3c Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Mon, 5 Apr 2021 14:01:11 -0400 Subject: [PATCH 3/3] DIIS + CG both use inner_product rather than dot --- src/TiledArray/dist_array.h | 8 ++++++++ src/TiledArray/math/linalg/conjgrad.h | 13 ++++++------- src/TiledArray/math/linalg/diis.h | 6 ++++-- 3 files changed, 18 insertions(+), 9 deletions(-) diff --git a/src/TiledArray/dist_array.h b/src/TiledArray/dist_array.h index 0378c01f43..7a1d90393a 100644 --- a/src/TiledArray/dist_array.h +++ b/src/TiledArray/dist_array.h @@ -1643,6 +1643,14 @@ auto dot(const DistArray& a, const DistArray& b) { .get(); } +template +auto inner_product(const DistArray& a, + const DistArray& b) { + return (a(detail::dummy_annotation(rank(a))) + .inner_product(b(detail::dummy_annotation(rank(b))))) + .get(); +} + template auto squared_norm(const DistArray& a) { return a(detail::dummy_annotation(rank(a))).squared_norm(); diff --git a/src/TiledArray/math/linalg/conjgrad.h b/src/TiledArray/math/linalg/conjgrad.h index 8ee0ffffe4..25ea53aff3 100644 --- a/src/TiledArray/math/linalg/conjgrad.h +++ b/src/TiledArray/math/linalg/conjgrad.h @@ -26,8 +26,8 @@ #ifndef TILEDARRAY_MATH_LINALG_CONJGRAD_H__INCLUDED #define TILEDARRAY_MATH_LINALG_CONJGRAD_H__INCLUDED -#include #include +#include #include "TiledArray/dist_array.h" namespace TiledArray::math::linalg { @@ -48,7 +48,7 @@ namespace TiledArray::math::linalg { /// \li value_type maxabs_value(const D&) /// \li void vec_multiply(D& a, const D& b) (element-wise multiply /// of \c a by \c b ) -/// \li value_type dot_product(const D& a, const D& b) +/// \li value_type inner_product(const D& a, const D& b) /// \li void scale(D&, value_type) /// \li void axpy(D& y,value_type a, const D& x) /// \li void assign(D&, const D&) @@ -71,7 +71,6 @@ struct ConjugateGradientSolver { /// elements in the residual. value_type operator()(F& a, const D& b, D& x, const D& preconditioner, value_type convergence_target = -1.0) { - std::size_t n = volume(preconditioner); const bool use_diis = false; @@ -133,10 +132,10 @@ struct ConjugateGradientSolver { unsigned int iter = 0; while (not converged) { // alpha_i = (r_i . z_i) / (p_i . A . p_i) - value_type rz_norm2 = dot(RR_i, ZZ_i); + value_type rz_norm2 = inner_product(RR_i, ZZ_i); a(PP_i, APP_i); - const value_type pAp_i = dot(PP_i, APP_i); + const value_type pAp_i = inner_product(PP_i, APP_i); const value_type alpha_i = rz_norm2 / pAp_i; // x_i += alpha_i p_i @@ -157,7 +156,7 @@ struct ConjugateGradientSolver { ZZ_i = RR_i; vec_multiply(ZZ_i, preconditioner); - const value_type rz_ip1_norm2 = dot(ZZ_i, RR_i); + const value_type rz_ip1_norm2 = inner_product(ZZ_i, RR_i); const value_type beta_i = rz_ip1_norm2 / rz_norm2; @@ -186,7 +185,7 @@ struct ConjugateGradientSolver { } // namespace TiledArray::math::linalg namespace TiledArray { - using TiledArray::math::linalg::ConjugateGradientSolver; +using TiledArray::math::linalg::ConjugateGradientSolver; } #endif // TILEDARRAY_MATH_LINALG_CONJGRAD_H__INCLUDED diff --git a/src/TiledArray/math/linalg/diis.h b/src/TiledArray/math/linalg/diis.h index b1335e3666..951fc6a11a 100644 --- a/src/TiledArray/math/linalg/diis.h +++ b/src/TiledArray/math/linalg/diis.h @@ -251,8 +251,10 @@ class DIIS { // and compute the most recent elements of B, B(i,j) = for (unsigned int i = 0; i < nvec - 1; i++) - B_(i, nvec - 1) = B_(nvec - 1, i) = dot(errors_[i], errors_[nvec - 1]); - B_(nvec - 1, nvec - 1) = dot(errors_[nvec - 1], errors_[nvec - 1]); + B_(i, nvec - 1) = B_(nvec - 1, i) = + inner_product(errors_[i], errors_[nvec - 1]); + B_(nvec - 1, nvec - 1) = + inner_product(errors_[nvec - 1], errors_[nvec - 1]); using std::abs; using std::sqrt; const auto current_error_2norm = sqrt(abs(B_(nvec - 1, nvec - 1)));