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