Feature/mrhs solvers

cmake/CPM.cmake

@@ -2,8 +2,8 @@
 #
 # SPDX-FileCopyrightText: Copyright (c) 2019-2023 Lars Melchior and contributors
 
-set(CPM_DOWNLOAD_VERSION 0.38.5)
-set(CPM_HASH_SUM "192aa0ccdc57dfe75bd9e4b176bf7fb5692fd2b3e3f7b09c74856fc39572b31c")
+set(CPM_DOWNLOAD_VERSION 0.40.2)
+set(CPM_HASH_SUM "c8cdc32c03816538ce22781ed72964dc864b2a34a310d3b7104812a5ca2d835d")
 
 if(CPM_SOURCE_CACHE)
   set(CPM_DOWNLOAD_LOCATION "${CPM_SOURCE_CACHE}/cpm/CPM_${CPM_DOWNLOAD_VERSION}.cmake")

include/accelerator.h

@@ -46,7 +46,12 @@ namespace quda
      * @param out Solution vector.
      * @param in Right-hand side.
      */
-    virtual void operator()(ColorSpinorField &out, ColorSpinorField &in)
+    virtual void operator()(cvector_ref<ColorSpinorField> &out, cvector_ref<const ColorSpinorField> &in) override
+    {
+      for (auto i = 0u; i < in.size(); i++) operator()(out[i], in[i]);
+    }
+
+    void operator()(ColorSpinorField &out, const ColorSpinorField &in)
     {
       if (transformer.trained) {
         transformer.apply(*base_solver, out, in);
@@ -64,7 +69,7 @@ namespace quda
      * @param null Solver to solve for null vectors.
      * @param in meta color spinor field.
      */
-    virtual void train_param(Solver &null, ColorSpinorField &in)
+    virtual void train_param(Solver &null, const ColorSpinorField &in) override
     {
       if (!active_training && !transformer.trained) {
         active_training = true;

include/blas_quda.h

@@ -30,6 +30,7 @@ namespace quda {
 
     inline void copy(cvector_ref<ColorSpinorField> &dst, cvector_ref<const ColorSpinorField> &src)
     {
+      if (dst.size() != src.size()) errorQuda("Mismatched vector sets %lu != %lu", dst.size(), src.size());
       for (auto i = 0u; i < src.size(); i++) { dst[i].copy(src[i]); }
     }
 
@@ -293,7 +294,7 @@ namespace quda {
 
     inline array<double, 2> max_deviation(const ColorSpinorField &x, const ColorSpinorField &y)
     {
-      return max_deviation(cvector_ref<const ColorSpinorField>(x), cvector_ref<const ColorSpinorField>(y));
+      return max_deviation(cvector_ref<const ColorSpinorField>(x), cvector_ref<const ColorSpinorField>(y))[0];
     }
 
     /**
@@ -302,13 +303,15 @@ namespace quda {
     */
     cvector<double> norm1(cvector_ref<const ColorSpinorField> &x);
 
+    inline double norm1(const ColorSpinorField &x) { return norm1(cvector_ref<const ColorSpinorField> {x})[0]; }
+
     /**
        @brief Compute the L2 norm (||x||^2) of a field
        @param[in] x The field we are reducing
     */
     cvector<double> norm2(cvector_ref<const ColorSpinorField> &x);
 
-    inline double norm2(const ColorSpinorField &x) { return norm2(cvector_ref<const ColorSpinorField> {x}); }
+    inline double norm2(const ColorSpinorField &x) { return norm2(cvector_ref<const ColorSpinorField> {x})[0]; }
 
     /**
        @brief Compute y += a * x and then (x, y)
@@ -319,6 +322,11 @@ namespace quda {
     cvector<double> axpyReDot(cvector<double> &a, cvector_ref<const ColorSpinorField> &x,
                               cvector_ref<ColorSpinorField> &y);
 
+    inline double axpyReDot(double a, const ColorSpinorField &x, ColorSpinorField &y)
+    {
+      return axpyReDot(cvector<double>(a), cvector_ref<const ColorSpinorField>(x), y)[0];
+    }
+
     /**
        @brief Compute the real-valued inner product (x, y)
        @param[in] x input vector
@@ -328,7 +336,7 @@ namespace quda {
 
     inline double reDotProduct(const ColorSpinorField &x, const ColorSpinorField &y)
     {
-      return reDotProduct(cvector_ref<const ColorSpinorField> {x}, cvector_ref<const ColorSpinorField> {y});
+      return reDotProduct(cvector_ref<const ColorSpinorField> {x}, cvector_ref<const ColorSpinorField> {y})[0];
     }
 
     /**
@@ -342,6 +350,12 @@ namespace quda {
     cvector<double> axpbyzNorm(cvector<double> &a, cvector_ref<const ColorSpinorField> &x, cvector<double> &b,
                                cvector_ref<const ColorSpinorField> &y, cvector_ref<ColorSpinorField> &z);
 
+    inline double axpbyzNorm(double a, const ColorSpinorField &x, double b, const ColorSpinorField &y,
+                             ColorSpinorField &z)
+    {
+      return axpbyzNorm(cvector<double>(a), cvector_ref<const ColorSpinorField>(x), cvector<double>(b), y, z)[0];
+    }
+
     /**
        @brief Compute y += a * x and then ||y||^2
        @param[in] a scalar multiplier
@@ -354,6 +368,11 @@ namespace quda {
       return axpbyzNorm(a, x, 1.0, y, y);
     }
 
+    inline double axpyNorm(double a, const ColorSpinorField &x, ColorSpinorField &y)
+    {
+      return axpyNorm(a, cvector_ref<const ColorSpinorField> {x}, cvector_ref<ColorSpinorField> {y})[0];
+    }
+
     /**
        @brief Compute the complex-valued inner product (x, y)
        @param[in] x input vector
@@ -363,7 +382,7 @@ namespace quda {
 
     inline Complex cDotProduct(const ColorSpinorField &x, const ColorSpinorField &y)
     {
-      return cDotProduct(cvector_ref<const ColorSpinorField> {x}, cvector_ref<const ColorSpinorField> {y});
+      return cDotProduct(cvector_ref<const ColorSpinorField> {x}, cvector_ref<const ColorSpinorField> {y})[0];
     }
 
     /**
@@ -373,6 +392,11 @@ namespace quda {
     */
     cvector<double4> cDotProductNormAB(cvector_ref<const ColorSpinorField> &x, cvector_ref<const ColorSpinorField> &y);
 
+    inline double4 cDotProductNormAB(const ColorSpinorField &x, const ColorSpinorField &y)
+    {
+      return cDotProductNormAB(cvector_ref<const ColorSpinorField> {x}, cvector_ref<const ColorSpinorField> {y})[0];
+    }
+
     /**
        @brief Return complex-valued inner product (x,y) and ||x||^2
        @param[in] x input vector
@@ -387,6 +411,11 @@ namespace quda {
       return a;
     }
 
+    inline double3 cDotProductNormA(const ColorSpinorField &x, const ColorSpinorField &y)
+    {
+      return cDotProductNormA(cvector_ref<const ColorSpinorField> {x}, cvector_ref<const ColorSpinorField> {y})[0];
+    }
+
     /**
        @brief Return complex-valued inner product (x,y) and ||y||^2
        @param[in] x input vector
@@ -401,6 +430,11 @@ namespace quda {
       return a;
     }
 
+    inline double3 cDotProductNormB(const ColorSpinorField &x, const ColorSpinorField &y)
+    {
+      return cDotProductNormB(cvector_ref<const ColorSpinorField> {x}, cvector_ref<const ColorSpinorField> {y})[0];
+    }
+
     /**
        @brief Apply the operation z += a * x + b * y, y -= b * w,
        compute complex-valued inner product (u, y) and ||y||^2
@@ -418,6 +452,14 @@ namespace quda {
                                                     cvector_ref<const ColorSpinorField> &w,
                                                     cvector_ref<const ColorSpinorField> &u);
 
+    inline double3 caxpbypzYmbwcDotProductUYNormY(const Complex &a, const ColorSpinorField &x, const Complex &b,
+                                                  ColorSpinorField &y, ColorSpinorField &z, const ColorSpinorField &w,
+                                                  const ColorSpinorField &u)
+    {
+      return caxpbypzYmbwcDotProductUYNormY(cvector<Complex>(a), cvector_ref<const ColorSpinorField>(x), b, y, z, w,
+                                            u)[0];
+    }
+
     /**
        @brief Compute y = a * x + b * y and then ||y||^2
        @param[in] a scalar multiplier
@@ -440,6 +482,11 @@ namespace quda {
       return caxpbyNorm(a, x, 1.0, y);
     }
 
+    inline double caxpyNorm(const Complex &a, const ColorSpinorField &x, ColorSpinorField &y)
+    {
+      return caxpyNorm(a, cvector_ref<const ColorSpinorField> {x}, cvector_ref<ColorSpinorField> {y})[0];
+    }
+
     /**
        @brief Compute y -= x and then ||y||^2
        @param[in] x input vector
@@ -450,6 +497,11 @@ namespace quda {
       return caxpbyNorm(1.0, x, -1.0, y);
     }
 
+    inline double xmyNorm(const ColorSpinorField &x, ColorSpinorField &y)
+    {
+      return xmyNorm(cvector_ref<const ColorSpinorField> {x}, cvector_ref<ColorSpinorField> {y})[0];
+    }
+
     /**
        @brief Compute z = a * b * x + y, x = a * x, and then ||z||^2
        @param[in] a scalar multiplier
@@ -461,6 +513,12 @@ namespace quda {
     cvector<double> cabxpyzAxNorm(cvector<double> &a, cvector<Complex> &b, cvector_ref<ColorSpinorField> &x,
                                   cvector_ref<const ColorSpinorField> &y, cvector_ref<ColorSpinorField> &z);
 
+    inline double cabxpyzAxNorm(double a, const Complex &b, ColorSpinorField &x, const ColorSpinorField &y,
+                                ColorSpinorField &z)
+    {
+      return cabxpyzAxNorm(cvector<double>(a), cvector<Complex>(b), cvector_ref<ColorSpinorField>(x), y, z)[0];
+    }
+
     /**
        @brief Compute y += a * x and the resulting complex-valued inner product (z, y)
        @param[in] a scalar multiplier
@@ -471,6 +529,11 @@ namespace quda {
     cvector<Complex> caxpyDotzy(cvector<Complex> &a, cvector_ref<const ColorSpinorField> &x,
                                 cvector_ref<ColorSpinorField> &y, cvector_ref<const ColorSpinorField> &z);
 
+    inline Complex caxpyDotzy(const Complex &a, const ColorSpinorField &x, ColorSpinorField &y, const ColorSpinorField &z)
+    {
+      return caxpyDotzy(cvector<Complex>(a), cvector_ref<const ColorSpinorField>(x), y, z)[0];
+    }
+
     /**
        @brief Compute y += a * x and then compute ||y||^2 and
        real-valued inner product (y_out, y_out-y_in)
@@ -481,6 +544,11 @@ namespace quda {
     cvector<double2> axpyCGNorm(cvector<double> &a, cvector_ref<const ColorSpinorField> &x,
                                 cvector_ref<ColorSpinorField> &y);
 
+    inline double2 axpyCGNorm(double a, const ColorSpinorField &x, ColorSpinorField &y)
+    {
+      return axpyCGNorm(cvector<double>(a), cvector_ref<const ColorSpinorField>(x), y)[0];
+    }
+
     /**
        @brief Computes ||x||^2, ||r||^2 and the MILC/FNAL heavy quark
        residual norm
@@ -492,7 +560,7 @@ namespace quda {
 
     inline double3 HeavyQuarkResidualNorm(const ColorSpinorField &x, const ColorSpinorField &r)
     {
-      return HeavyQuarkResidualNorm(cvector_ref<const ColorSpinorField>(x), cvector_ref<const ColorSpinorField>(r));
+      return HeavyQuarkResidualNorm(cvector_ref<const ColorSpinorField>(x), cvector_ref<const ColorSpinorField>(r))[0];
     }
 
     /**
@@ -510,7 +578,7 @@ namespace quda {
                                              const ColorSpinorField &r)
     {
       return xpyHeavyQuarkResidualNorm(cvector_ref<const ColorSpinorField>(x), cvector_ref<const ColorSpinorField>(y),
-                                       cvector_ref<const ColorSpinorField>(r));
+                                       cvector_ref<const ColorSpinorField>(r))[0];
     }
 
     /**
@@ -522,6 +590,11 @@ namespace quda {
     cvector<double3> tripleCGReduction(cvector_ref<const ColorSpinorField> &x, cvector_ref<const ColorSpinorField> &y,
                                        cvector_ref<const ColorSpinorField> &z);
 
+    inline double3 tripleCGReduction(const ColorSpinorField &x, const ColorSpinorField &y, const ColorSpinorField &z)
+    {
+      return tripleCGReduction(cvector_ref<const ColorSpinorField>(x), y, z)[0];
+    }
+
     /**
        @brief Computes ||x||^2, ||y||^2, the real-valued inner product (y, z), and ||z||^2
        @param[in] x input vector
@@ -531,6 +604,11 @@ namespace quda {
     cvector<double4> quadrupleCGReduction(cvector_ref<const ColorSpinorField> &x, cvector_ref<const ColorSpinorField> &y,
                                           cvector_ref<const ColorSpinorField> &z);
 
+    inline double4 quadrupleCGReduction(const ColorSpinorField &x, const ColorSpinorField &y, const ColorSpinorField &z)
+    {
+      return quadrupleCGReduction(cvector_ref<const ColorSpinorField>(x), y, z)[0];
+    }
+
     /**
        @brief Computes z = x, w = y, x += a * y, y -= a * v and ||y||^2
        @param[in] a scalar multiplier
@@ -544,6 +622,12 @@ namespace quda {
                                          cvector_ref<ColorSpinorField> &y, cvector_ref<ColorSpinorField> &z,
                                          cvector_ref<ColorSpinorField> &w, cvector_ref<const ColorSpinorField> &v);
 
+    inline double quadrupleCG3InitNorm(double a, ColorSpinorField &x, ColorSpinorField &y, ColorSpinorField &z,
+                                       ColorSpinorField &w, const ColorSpinorField &v)
+    {
+      return quadrupleCG3InitNorm(cvector<double>(a), cvector_ref<ColorSpinorField>(x), y, z, w, v)[0];
+    }
+
     /**
        @brief Computes x = b * (x + a * y) + ( 1 - b) * z,
        y = b * (y + a * v) + (1 - b) * w, z = x_in, w = y_in, and
@@ -560,6 +644,12 @@ namespace quda {
                                            cvector_ref<ColorSpinorField> &y, cvector_ref<ColorSpinorField> &z,
                                            cvector_ref<ColorSpinorField> &w, cvector_ref<const ColorSpinorField> &v);
 
+    inline double quadrupleCG3UpdateNorm(double a, double b, ColorSpinorField &x, ColorSpinorField &y,
+                                         ColorSpinorField &z, ColorSpinorField &w, const ColorSpinorField &v)
+    {
+      return quadrupleCG3UpdateNorm(cvector<double>(a), b, cvector_ref<ColorSpinorField>(x), y, z, w, v)[0];
+    }
+
     namespace block
     {
 

include/color_spinor_field.h

@@ -943,6 +943,36 @@ namespace quda
   void resize(std::vector<ColorSpinorField> &v, size_t new_size, QudaFieldCreate create,
               const ColorSpinorField &src = ColorSpinorField());
 
+  /**
+     @brief Create a vector of fields that aliases another vector of
+     fields' storage.  The alias field can use a different precision
+     than this field, though it cannot be greater.  This
+     functionality is useful for the case where we have multiple
+     temporaries in different precisions, but do not need them
+     simultaneously.  Use this functionality with caution.
+     @param[out] alias The vector of aliased fields
+     @param[in] v The vector of fields to alias
+     @param[in] param Parameters for the alias field
+  */
+  void create_alias(cvector_ref<ColorSpinorField> &alias, cvector_ref<const ColorSpinorField> &v,
+                    const ColorSpinorParam &param = ColorSpinorParam());
+
+  /**
+     @brief Create a vector of fields that aliases another vector of
+     fields' storage.  The alias field can use a different precision
+     than this field, though it cannot be greater.  This functionality
+     is useful for the case where we have multiple temporaries in
+     different precisions, but do not need them simultaneously.  This
+     variant is used with std::vector as opposed to vector_ref, and
+     allows for correct resizing.  Use this functionality with
+     caution.
+     @param[out] alias The vector of aliased fields
+     @param[in] v The vector of fields to alias
+     @param[in] param Parameters for the alias field
+  */
+  void create_alias(std::vector<ColorSpinorField> &alias, cvector_ref<const ColorSpinorField> &v,
+                    const ColorSpinorParam &param = ColorSpinorParam());
+
   void copyGenericColorSpinor(ColorSpinorField &dst, const ColorSpinorField &src, QudaFieldLocation location,
                               void *Dst = nullptr, const void *Src = nullptr);
 

include/dirac_quda.h

@@ -542,7 +542,6 @@ namespace quda {
     DiracWilson(const DiracWilson &dirac);
     DiracWilson(const DiracParam &param, const int nDims); // to correctly adjust face for DW and non-deg twisted mass
 
-    virtual ~DiracWilson();
     DiracWilson& operator=(const DiracWilson &dirac);
 
     virtual void Dslash(cvector_ref<ColorSpinorField> &out, cvector_ref<const ColorSpinorField> &in,
@@ -592,7 +591,6 @@ namespace quda {
   public:
     DiracWilsonPC(const DiracParam &param);
     DiracWilsonPC(const DiracWilsonPC &dirac);
-    virtual ~DiracWilsonPC();
     DiracWilsonPC& operator=(const DiracWilsonPC &dirac);
 
     void M(cvector_ref<ColorSpinorField> &out, cvector_ref<const ColorSpinorField> &in) const override;
@@ -1962,7 +1960,6 @@ namespace quda {
        @param[in] param Parameters defining this operator
      */
     DiracCoarse(const DiracCoarse &dirac, const DiracParam &param);
-    virtual ~DiracCoarse();
 
     virtual bool isCoarse() const override { return true; }
 
@@ -2108,8 +2105,6 @@ namespace quda {
      */
     DiracCoarsePC(const DiracCoarse &dirac, const DiracParam &param);
 
-    virtual ~DiracCoarsePC();
-
     /**
        @brief Apply preconditioned Dslash out = (D * in)
        @param[out] out Output field