Evaleev/fix/nonintrusive solver adaptors

src/TiledArray/dist_array.h

@@ -1643,6 +1643,14 @@ auto dot(const DistArray<Tile, Policy>& a, const DistArray<Tile, Policy>& b) {
       .get();
 }
 
+template <typename Tile, typename Policy>
+auto inner_product(const DistArray<Tile, Policy>& a,
+                   const DistArray<Tile, Policy>& b) {
+  return (a(detail::dummy_annotation(rank(a)))
+              .inner_product(b(detail::dummy_annotation(rank(b)))))
+      .get();
+}
+
 template <typename Tile, typename Policy>
 auto squared_norm(const DistArray<Tile, Policy>& a) {
   return a(detail::dummy_annotation(rank(a))).squared_norm();

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 <typename Tile, typename Policy>
 inline void vec_multiply(DistArray<Tile, Policy>& a1,
                          const DistArray<Tile, Policy>& a2) {
@@ -59,4 +62,51 @@ inline void axpy(DistArray<Tile, Policy>& y, S alpha,
 
 }  // namespace TiledArray::math::linalg
 
+namespace Eigen {
+
+// freestanding adaptors for Eigen::MatrixBase needed by solvers like DIIS
+
+template <typename Derived>
+inline void vec_multiply(Eigen::MatrixBase<Derived>& a1,
+                         const Eigen::MatrixBase<Derived>& a2) {
+  a1.array() *= a2.array();
+}
+
+template <typename Derived, typename S>
+inline void scale(Eigen::MatrixBase<Derived>& a, S scaling_factor) {
+  using numeric_type = typename Eigen::MatrixBase<Derived>::value_type;
+  a.array() *= numeric_type(scaling_factor);
+}
+
+template <typename Derived>
+inline void zero(Eigen::MatrixBase<Derived>& a) {
+  a = Derived::Zero(a.rows(), a.cols());
+}
+
+template <typename Derived, typename S>
+inline void axpy(Eigen::MatrixBase<Derived>& y, S alpha,
+                 const Eigen::MatrixBase<Derived>& x) {
+  using numeric_type = typename Eigen::MatrixBase<Derived>::value_type;
+  y.array() += numeric_type(alpha) * x.array();
+}
+
+template <typename Derived>
+inline auto dot(const Eigen::MatrixBase<Derived>& l,
+                const Eigen::MatrixBase<Derived>& r) {
+  return l.adjoint().dot(r);
+}
+
+template <typename Derived>
+inline auto inner_product(const Eigen::MatrixBase<Derived>& l,
+                          const Eigen::MatrixBase<Derived>& r) {
+  return l.dot(r);
+}
+
+template <typename Derived>
+inline auto norm2(const Eigen::MatrixBase<Derived>& m) {
+  return m.template lpNorm<2>();
+}
+
+}  // namespace Eigen
+
 #endif  // TILEDARRAY_MATH_LINALG_BASIC_H__INCLUDED

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 <TiledArray/math/linalg/diis.h>
 #include <TiledArray/math/linalg/basic.h>
+#include <TiledArray/math/linalg/diis.h>
 #include "TiledArray/dist_array.h"
 
 namespace TiledArray::math::linalg {
@@ -48,7 +48,7 @@ namespace TiledArray::math::linalg {
 ///   \li <tt> value_type maxabs_value(const D&) </tt>
 ///   \li <tt> void vec_multiply(D& a, const D& b) </tt> (element-wise multiply
 ///   of \c a by \c b )
-///   \li <tt> value_type dot_product(const D& a, const D& b) </tt>
+///   \li <tt> value_type inner_product(const D& a, const D& b) </tt>
 ///   \li <tt> void scale(D&, value_type) </tt>
 ///   \li <tt> void axpy(D& y,value_type a, const D& x) </tt>
 ///   \li <tt> void assign(D&, const D&) </tt>
@@ -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

src/TiledArray/math/linalg/diis.h

@@ -251,8 +251,10 @@ class DIIS {
 
     // and compute the most recent elements of B, B(i,j) = <ei|ej>
     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)));

src/TiledArray/tensor/complex.h

@@ -40,7 +40,8 @@ namespace detail {
 /// \param r The real scalar
 /// \return `r`
 template <typename R,
-          typename std::enable_if<!is_complex<R>::value>::type* = nullptr>
+          typename std::enable_if<is_numeric_v<R> &&
+                                  !is_complex<R>::value>::type* = nullptr>
 TILEDARRAY_FORCE_INLINE R conj(const R r) {
   return r;
 }
@@ -61,7 +62,8 @@ TILEDARRAY_FORCE_INLINE std::complex<R> conj(const std::complex<R> z) {
 /// \tparam R A numeric type
 /// \return `r`
 template <typename L, typename R,
-          typename std::enable_if<!is_complex<L>::value>::type* = nullptr>
+          typename std::enable_if<is_numeric_v<L> &&
+                                  !is_complex<L>::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 <typename L, typename R,
-          typename std::enable_if<is_complex<L>::value>::type* = nullptr>
+          typename std::enable_if<is_numeric_v<L> &&
+                                  is_complex<L>::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 <typename R,
-          typename std::enable_if<!is_complex<R>::value>::type* = nullptr>
+          typename std::enable_if<is_numeric_v<R> &&
+                                  !is_complex<R>::value>::type* = nullptr>
 TILEDARRAY_FORCE_INLINE R norm(const R r) {
   return r * r;
 }