From 9cab78f4937f6f9eb0f6bf1d128566c8589f4ac7 Mon Sep 17 00:00:00 2001 From: Eloy Romero Date: Tue, 5 Apr 2022 14:03:28 -0400 Subject: [PATCH 001/260] Towards the first version of mgproton (in progress); New class to deal with options (doesn't compile); New functions to fill the tensors. --- ...inline_baryon_matelem_colorvec_superb_w.cc | 2 +- lib/util/ferm/mgproton.h | 414 +++++++++++ lib/util/ferm/superb_contractions.h | 87 ++- lib/util/ferm/superb_options.h | 684 ++++++++++++++++++ 4 files changed, 1182 insertions(+), 5 deletions(-) create mode 100644 lib/util/ferm/mgproton.h create mode 100644 lib/util/ferm/superb_options.h diff --git a/lib/meas/inline/hadron/inline_baryon_matelem_colorvec_superb_w.cc b/lib/meas/inline/hadron/inline_baryon_matelem_colorvec_superb_w.cc index dabdd52f1..8867e2dc3 100644 --- a/lib/meas/inline/hadron/inline_baryon_matelem_colorvec_superb_w.cc +++ b/lib/meas/inline/hadron/inline_baryon_matelem_colorvec_superb_w.cc @@ -396,7 +396,7 @@ namespace Chroma } //---------------------------------------------------------------------------- - //! Normalize just one displacement array and return std::vector + //! Return multi1d from std::vector multi1d tomulti1d(const std::vector& orig) { multi1d r(orig.size()); diff --git a/lib/util/ferm/mgproton.h b/lib/util/ferm/mgproton.h new file mode 100644 index 000000000..83f5be383 --- /dev/null +++ b/lib/util/ferm/mgproton.h @@ -0,0 +1,414 @@ +// -*- C++ -*- +/*! \file + * \brief Multigrid prototype next + * + * Hadron spectrum calculations utilities + */ + +#ifndef __INCLUDE_MGPROTON__ +#define __INCLUDE_MGPROTON__ + +#include "chromabase.h" +#include "util/ferm/superb_contractions.h" + +#ifdef BUILD_SB +namespace Chroma +{ + + namespace SB + { + + /// Verbosity level for solvers + enum Verbosity { + NoOutput, ///< Print nothing + JustSummary, ///< Print summary of the performance, eg., number of iterations and final residual + Detailed, ///< Print progress + VeryDetailed ///< Print whatever + }; + + /// Return the map from string to Verbosity values + inline const std::map& getVerbosityMap() + { + static const std::map m{{"NoOutput", NoOutput}, + {"JustSummary", JustSummary}, + {"Detailed", Detailed}, + {"VeryDetailed", VeryDetailed}, + {"false", NoOutput}}; + return m; + } + + /// Representation of an operator, function of type tensor -> tensor where the input and the + /// output tensors have the same dimensions + + template + using OperatorFun = std::function(const Tensor&)>; + + /// Representation of an operator together with a map to convert from domain labels (columns) to + /// image labels (rows) + + template + struct Operator : public OperatorFun { + remap map; + Operator(const OperatorFun& op, const remap& map) + : OperatorFun(op), map(map) + { + } + }; + + namespace detail + { + /// Return an order with all domain labels + /// \param map: from domain labels to image labels + /// \param exclude: do not return the given labels + + std::string getDomain(const remap& map, const std::string& exclude = "") + { + std::string o; + for (const auto& it : map) + if (std::find(exclude.begin(), exclude.end(), it.first) == exclude.end()) + o += it.first; + return o; + } + + /// Return an order with all domain labels + /// \param map: from domain labels to image labels + /// \param exclude: do not return the given labels + + std::string getImage(const remap& map, const std::string& exclude = "") + { + std::string o; + for (const auto& it : map) + if (std::find(exclude.begin(), exclude.end(), it.second) == exclude.end()) + o += it.second; + return o; + } + + /// Return the inverse map + /// \param map: from domain labels to image labels + + remap reverse(const remap& map) + { + remap o; + for (const auto& it : map) + o.insert({it.second, it.first}); + return o; + } + } + + namespace detail + { + + } + + /// Orthonormalize W against V, W_out <- (I-V*V')*W_in*R, with R such that W_out'*W_out = I, + /// for each combination of values of the order_t dimensions, or once if it is empty. + /// \param V: orthogonal matrix to orthogonalize against + /// \param W: matrix to orthogonalize + /// \param order_t: dimension labels that do not participate in the orthogonalization + /// \param order_rows: dimension labels that are the rows of the matrices V and W + /// \param order_cols: dimension labels that are the columns of the matrices V and W + + template + void ortho(const Tensor& V, Tensor& W, Maybe order_t, + Maybe order_rows, Maybe order_cols, + unsigned int max_its = 3) + { + // Check that V and W are compatible, excepting the column dimensions + if (V) + detail::check_compatible(V, W, none, order_cols); + + // Find the ordering + std::string torder, rorder, Wcorder, Vcorder; + detail::get_mat_ordering(W, order_t, order_rows, order_cols, torder, rorder, Wcorder); + if (V) + detail::get_mat_ordering(V, order_t, order_rows, order_cols, torder, rorder, Vcorder); + + // Create an alternative view of W with different labels for the column + remap Wac = getNewLabels(Wcorder, toString(V.kvdim()) + toString(W.kvdim())); + std::string Wacorder = detail::update_order(Wcorder, Wac); + auto Wa = W.rename_dims(Wac); + + for (unsigned int i = 0; i < max_its; ++i) + { + // W = W - V*(V'*W) + if (V) + contract(V.scale(-1), contract(V.conj(), W, rorder), Vcorder, AddTo, W); + + // W = W/chol(Wa'*W) + contract(W, + inverse(chol(contract(Wa.conj(), W, rorder), torder, Wacorder, Wcorder), torder, + Wacorder, Wcorder), + CopyTo, W); + } + + // For each column in W do classic Gram-Schmidt + } + + /// Solve iteratively op * y = x using FGMRES + /// \param op: problem matrix + /// \param prec: preconditioner + /// \param x: input right-hand-sides + /// \param y: the solution vectors + /// \param order_op_t: op and prec dimension labels that do not participate in the solver + /// \param order_op_rows: dimension labels that are the rows of op and y, and the columns of prec + /// \param order_op_cols: dimension labels that are the columns of op and the rows of y and prec + /// \param order_xy_cols: dimension labels that are the columns of x and y + /// \param max_basis_size: maximum rank of the search subspace basis per t + /// \param tol: maximum tolerance + /// \param max_its: maximum number of iterations + /// \param error_if_not_converged: throw an error if the tolerance was not satisfied + /// \param max_simultaneous_rhs: maximum number of right-hand-sides solved simultaneously; all by default + /// \param verb: verbosity level + /// \param prefix: prefix printed before every line + + template + void fgmres(const Operator& op, Maybe&> prec, + const Tensor& x, Tensor& y, + const std::string& order_t, const std::string& order_rows, + const std::string& order_cols, unsigned int max_basis_size, double tol, + unsigned int max_its = 0, bool error_if_not_converged = true, + unsigned int max_simultaneous_rhs = 0, Verbosity verb = NoOutput, + std::string prefix = "") + { + // Find volume of columns in W + + // For each column in W do classic Gram-Schmidt + } + + template + Operator getSolver(const Operator& op, const Options& ops); + + namespace detail + { + /// Function that gets an operator and options and returns an operator that approximates the inverse + template + using Solver = + std::function(const Operator& op, const Options& ops)>; + + /// Returns an approximate inverse operator given an operator, options, and a list of solvers + /// \param solvers: map of solvers + /// \param op: operator to make the inverse of + /// \param ops: options to select the solver from `solvers` and influence the solver construction + + template + Operator getSolver(const std::map>& solvers, + const Operator& op, const Options& ops) + { + std::string type = ops.getValue("type", StringOption{""}).getString(); + if (type.size() == 0) + throw std::runtime_error("Error getting a solver from option `" + ops.prefix + + "': missing tag `type'"); + if (solvers.count(type) == 0) + { + std::string supported_tags; + for (const auto& it : solvers) + supported_tags += it->first + " "; + throw std::runtime_error("Error getting a solver from option `" + ops.prefix + + "': unsupported tag value `" + type + + "'; supported values: " + supported_tags); + } + + return solvers[type](op, ops); + } + + /// Returns a FGMRES solver + /// \param op: operator to make the inverse of + /// \param ops: options to select the solver from `solvers` and influence the solver construction + + template + Operator getFGMRESSolver(const Operator& op, const Options& ops) + { + // Get preconditioner + Maybe> prec = none; + const Options& precOps = ops.getValue("prec", NoneOption{}, Dictionary); + if (precOps) prec = Maybe>(getSolver(op, precOps)); + + // Get the remainder options + std::string order_t = getOption(ops, "order_t", ""); + std::string order_rows = + getOption(ops, "order_rows", getImage(op.map, order_t)); + std::string order_cols = + getOption(ops, "order_cols", getDomain(op.map, order_t)); + unsigned int max_basis_size = getOption(ops, "max_basis_size", 0); + double tol = getOption(ops, "tol", 0.0); + unsigned int max_its = getOption(ops, "max_its", 0); + bool error_if_not_converged = getOption(ops, "max_basis_size", true); + unsigned int max_simultaneous_rhs = getOption(ops, "max_simultaneous_rhs", 0); + Verbosity verb = getOption(ops, "verbosity", getVerbosityMap(), NoOutput); + std::string prefix = getOption(ops, "prefix", ""); + + // Return the solver + const Operator op_ = op; // copy by value of the operator + return {[=](const Tensor& x) { + auto y = x.like_this(); + fgmres(op_, prec, x, y, order_t, order_rows, order_cols, max_basis_size, tol, + max_its, error_if_not_converged, max_simultaneous_rhs, verb, prefix); + return y; + }, + reverse(op.map)}; + } + + /// Returns the \gamma_5 or a projection for a given number of spins + /// \param ns: number of spins + + template + Tensor<2, COMPLEX> getGamma5(int ns, DeviceHost dev = OnDefaultDevice) + { + if (ns == 1) + { + Tensor<2, COMPLEX> r("ij", {1, 1}, OnHost, OnEveryoneReplicated); + r.set({{0, 0}}, COMPLEX{1}); + return r.make_sure(none, dev); + } + else if (ns == Ns) + { + return SB::Gamma(Ns * Ns - 1); + } + else if (ns == 2) + { + Tensor<1, COMPLEX> ones("j", {Ns}, OnHost, OnEveryoneReplicated); + for (int i = 0; i < Ns; ++i) + ones.set({{i}}, COMPLEX{1}); + Tensor<1, COMPLEX> p("i", {Ns}, OnHost, OnEveryoneReplicated); + p.contract(getGamma5(Ns), {}, NotConjugate, ones, {}, NotConjugate); + Tensor<2, COMPLEX> r("ij", {2, 2}, OnHost, OnEveryoneReplicated); + r.set_zero(); + r.set({{0, 0}}, COMPLEX{1}); + r.set({{1, 1}}, COMPLEX{-1}); + return r.make_sure(none, dev); + } + else + { + throw std::runtime_error("Error in getGamma5: Unsupported spin number"); + } + } + + /// Returns the prolongator constructed from + /// \param solvers: map of solvers + /// \param op: operator to make the inverse of + /// \param ops: options to select the solver from `solvers` and influence the solver construction + + template + Operator getMGProlongator(const Operator& op, + unsigned int num_null_vecs, + const std::map& blocking, + const Operator& null_solver) + { + // Create num_null_vecs random vectors and solve them + auto b = op.domain_eg.like_this(none, {'n', num_null_vecs}); + urand(b, -1, 1); + auto nv = null_solver(std::move(b)); + + // Do chirality splitting nv2 = [ nv * gpos, nv * gneg ] + // TODO: experiment without chirality splitting + int ns = nv.kvdims()['s']; + if (ns != 2 && ns != Ns) + throw std::runtime_error("Error in getMGProlongator: Unsupported spin number"); + + auto nv2 = nv.like_this(none, {'n', num_null_vecs * 2}); + auto g5 = getGamma5(ns, OnHost), gpos = g5.clone(), gneg = g5.clone(); + for (int i = 0; i < Ns; ++i) // make diagonal entries of gpos all positive or zero + gpos.set({{i, i}}, g5.get({{i, i}} + 1)); + for (int i = 0; i < Ns; ++i) // make diagonal entries of gneg all negative or zero + gneg.set({{i, i}}, g5.get({{i, i}} - 1)); + nv2.kvslice_from_size({}, {{'n', num_null_vecs}}) + .contract(g5pos, {}, NotConjugate, nv, {}, NotConjugate); + nv2.kvslice_from_size({{'n', num_null_vecs}}, {{'n', num_null_vecs}}) + .contract(g5neg, {}, NotConjugate, nv, {}, NotConjugate); + + // Do the blocking + auto nv_blk = nv.split_dimension('x', "Wx", blocking['x']) + .split_dimension('y', "Yy", blocking['y']) + .split_dimension('z', "Zz", blocking['z']) + .split_dimension('t', "Tt", blocking['t']) + .split_dimension('n', "Nn", num_null_vecs); + + // Do the orthogonalization on each block and chirality + ortho(Tensor<1, COMPLEX>(), nv_blk, "xyztn", "WYZTXsc", "N"); + } + + + /// Returns a MG preconditioner. + /// + /// It returns an approximation of Op^{-1} = Op^{-1}*Q + Op^{-1}(I-Q), where Q is a projector + /// on the left singular space of Op. (MG can be derived also from P*Op^{-1} + (I-P)*Op^{-1}, + /// where P is on the right singular space, producing pre-smoothers instead.) + /// The approximation is constructed using an oblique projector and doing the inversions + /// approximately: + /// 1) Q = Op*V*(W'*Op*V)^{-1}*W', where W and V are left and right singular spaces. + /// 2) [Op^{-1}*Q + Op^{-1}(I-Q)]*x \approx + /// V*solver(W'*Op*V, W'*x) + solver(Op, x - Op*V*solver(W'*Op*V, W'*x)) + /// Note that if Op is \gamma_5-Hermitian, and W=\gamma_5*V, and \gamma_5 commutes with V + /// (\gamma_5*V = V*\gamma_5^coarse), then the projector reduces to Q=Op*V(V'*Op*V)^{-1}*V', + /// and, this coarse operator is easier to solve than V'\gamma_5*Op*V. + /// + /// \param op: operator to make the inverse of + /// \param ops: options to select the solver from `solvers` and influence the solver construction + + template + Operator getMGPrec(const Operator& op, const Options& ops) + { + // Get prolongator, V + unsigned int num_null_vecs = getOption(ops, "num_null_vecs"); + std::vector blocking_v = + getOption>(ops, "blocking"); + if (blocking_v.size() != Nd) + throw std::runtime_error("Error generating MG preconditioner for prefix `" + ops.prefix + + "': the blocking should be a vector with four elements"); + std::map blocking{ + {'x', blocking_v[0]}, {'y', blocking_v[1]}, {'z', blocking_v[2]}, {'t', blocking_v[3]}}; + const Options& nvInvOps = ops.getValue("prec_null_vecs", none, Dictionary); + const Operator nullSolver = getSolver(op, nvInvOps); + auto V = getMGProlongator(op, num_null_vecs, blocking, nullSolver); + + // Compute the coarse operator + const Operator op_c = getCoarseOperator(op, V); + + // Get the solver for the projector + const Options& coarseSolverOps = ops.getValue("solver_coarse", none, Dictionary); + const Operator coarseSolver = getSolver(op_c, coarseSolverOps); + + // Get the solver for the smoother + const Options& opSolverOps = ops.getValue("solver_smoother", none, Dictionary); + const Operator opSolver = getSolver(op, opSolverOps); + + // Return the solver + const Operator op_ = op; // copy by value of the operator + return {[=](Tensor x) { + // y0 = V*solver(V'*Op*V, V'x) + auto y0 = V.contract(coarseSolver(V.conj_contract(x))); + + // x1 = x - op*y0 + auto x1 = op_(y0.scale(-1)); + x1 += std::move(x); + + // y = y0 + solver(Op, x1) + auto y = opSolver(std::move(x1)); + y += y0; + + return y; + }, + reverse(op.map)}; + } + } + + /// Returns an operator that approximate the inverse of a given operator + /// \param op: operator to make the inverse of + /// \param ops: options to select the solver from `solvers` and influence the solver construction + + template + Operator getSolver(const Operator& op, const Options& ops) + { + static const std::map> solvers{ + {"fgmres", detail::getFGMRESSolver}, // flexible GMRES + {"mg", detail::getMGPrec} // Multigrid + }; + + return detail::getSolver(solvers, op, ops); + } + } +} + +#endif // BUILD_SB + +#endif // __INCLUDE_MGPROTON__ diff --git a/lib/util/ferm/superb_contractions.h b/lib/util/ferm/superb_contractions.h index e2183e5aa..b68789cde 100644 --- a/lib/util/ferm/superb_contractions.h +++ b/lib/util/ferm/superb_contractions.h @@ -32,6 +32,7 @@ # include # include # include +# include # include # include # include @@ -120,6 +121,12 @@ namespace Chroma return opt_val.first; } + /// Return whether it has been initialized with a value + explicit operator bool() const noexcept + { + return hasSome(); + } + /// Return the value if it has been initialized with some T getSome() const { @@ -1115,7 +1122,7 @@ namespace Chroma /// Return whether the tensor is not empty explicit operator bool() const noexcept { - return superbblas::detail::volume(size) > 0; + return volume() > 0; } // Return whether from != 0 or size != dim @@ -1683,7 +1690,7 @@ namespace Chroma if (dist != OnMaster && dist != Local) return false; - if (superbblas::detail::volume(size) > 0 && N > 1) + if (volume() > 0 && N > 1) { bool non_full_dim = false; // some dimension is not full for (unsigned int i = 0; i < N - 1; ++i) @@ -1758,7 +1765,7 @@ namespace Chroma throw std::runtime_error("Not allowed for tensor with a scale not being one"); // Only on primary node read the data - std::size_t vol = superbblas::detail::volume(size); + std::size_t vol = volume(); std::size_t disp = detail::coor2index(from, dim, strides); std::size_t word_size = sizeof(typename detail::WordType::type); bin.readArrayPrimaryNode((char*)&data.get()[disp], word_size, sizeof(T) / word_size * vol); @@ -1799,7 +1806,7 @@ namespace Chroma ss << "% " << repr(data.get()) << std::endl; ss << "% dist=" << p->p << std::endl; ss << name << "=reshape(["; - std::size_t vol = superbblas::detail::volume(size); + std::size_t vol = volume(); for (std::size_t i = 0; i < vol; ++i) { //using detail::repr::operator<<; @@ -2235,6 +2242,78 @@ namespace Chroma } }; + /// Modify the content of a tensor with the result of a function + /// \param t: tensor template + /// \param func: function () -> COMPLEX + /// \param threaded: whether to run threaded + + template + void fillWithCPUFuncNoArgs(const Tensor& t, Func func, bool threaded = true) + { + auto l = t.getLocal().like_this(none, OnHost); + std::size_t vol = l.volume(); + COMPLEX* ptr = r.data.get(); + + if (threaded) + { +# ifdef _OPENMP +# pragma omp parallel for schedule(static) +# endif + for (std::size_t i = 0; i < vol; ++i) + ptr[i] = func(); + } + else + { + for (std::size_t i = 0; i < vol; ++i) + ptr[i] = func(); + } + + l.copyTo(t.getLocal()); + } + + namespace detail + { + inline std::mt19937_64& getSeed() + { + // This is quick and dirty and nonreproducible if the lattice is distributed + // among the processes is different ways. + static std::mt19937_64 twister_engine(10 + QMP_get_node_number()); + return twister_engine; + } + } + + /// Modify with complex random uniformly distributed numbers with the real and the imaginary part between [a,b] + /// \param t: tensor to fill with random numbers + /// \param a: minimum random value + /// \param b: maximum random value + + template ::value, bool>::type = true> + void urand(const Tensor& t, T::value_type a=0, T::value_type b=1) + { + std::uniform_real_distribution d(a, b); + fillWithCPUFuncNoArgs( + t, + []() { + return {d(detail::getSeed()), d(detail::getSeed())}; + }, + false); + } + + /// Modify with random uniformly distributed numbers between [a,b] + /// \param t: tensor to fill with random numbers + /// \param a: minimum random value + /// \param b: maximum random value + + template ::value, bool>::type = true> + void urand(const Tensor& t, T::value_type a=0, T::value_type b=1) + { + std::uniform_real_distribution d(a, b); + fillWithCPUFuncNoArgs( + t, []() { return d(detail::getSeed()); }, false); + } + /// Return a tensor filled with the value of the function applied to each element /// \param order: dimension labels, they should start with "xyztX" /// \param size: length of each dimension diff --git a/lib/util/ferm/superb_options.h b/lib/util/ferm/superb_options.h new file mode 100644 index 000000000..21d0706dc --- /dev/null +++ b/lib/util/ferm/superb_options.h @@ -0,0 +1,684 @@ +// -*- C++ -*- +/*! \file + * \brief Multigrid prototype next + * + * Hadron spectrum calculations utilities + */ + +#ifndef __INCLUDE_SUPERB_OPTIONS__ +#define __INCLUDE_SUPERB_OPTIONS__ + +#include "chromabase.h" +#include "util/ferm/superb_contractions.h" + +#ifdef BUILD_SB +namespace Chroma +{ + + namespace SB + { + + namespace detail { + /// Return the previous lines up this one + /// \param file: content + /// \param char_num: character index of the last line to print + /// \param num_prev_lines: maximum number of previous lines to print + /// \param prefix: string to print previous to each line + + void get_prev_lines(const std::string& file, unsigned int char_num, + unsigned int num_prev_lines = 0, const std::string prefix = "") const + { + if (char_num > file.size()) + throw std::runtime_error("Erroneous character number"); + std::size_t line_num = std::count(file.begin(), file.begin() + char_num, '\n') + 1; + auto p = file.begin(); + for (unsigned int l = 1; l <= line_num && p != file.end(); ++l) + { + auto p1 = std::find(p, file.end(), '\n'); + if (l + num_prev_lines >= line_num) + QDPIO::cout << prefix << std::string(p, p1) << std::endl; + p = p1 + (p1 != file.end() ? 1 : 0); + } + } + } + + /// Class for storing options + struct Options { + /// Get track of the path of this option on a set of options + std::string prefix; + /// Content of the file where the option comes from + std::shared_ptr file; + /// First line number in `file` associated to the option + unsigned int char_num; + /// Whether this option has been checked + mutable bool visited; + + protected: + Options() : char_num(0), visited(false) + { + } + + /// Copy `prefix`, `file` and `char_num` from a given option + /// \param op: option to copy the info + + void copyFileInfo(const Options& op) + { + prefix = op.prefix; + file = op.file; + char_num = op.char_num; + } + + /// Throw an error + /// \param s: error message + + void throw_error(const std::string& err_msg) const + { + if (file) + throw std::runtime_error("Error at prefix `" + prefix + "'(l. " + + std::to_string(line_num) + "): " + err_msg + "\n" + + detail::get_prev_lines(*file, char_num, 5, "| ")); + else + throw std::runtime_error("Error at prefix `" + prefix + ": " + err_msg); + } + + public: + /// Type of the option + enum Type { None, String, Double, Vector, Dictionary }; + + /// Return the type of this option + virtual Type getType() const = 0; + + /// Return if this options isn't None + explicit operator bool() const noexcept + { + return getType() != None; + } + + /// Return the string content of the option + virtual std::string getString() const + { + throw_error("expected the value to be a string"); + } + + /// Return the double content of the option + virtual double getDouble() const + { + throw_error("expected the value to be a double"); + } + + /// Return the integer content of the option + virtual int getInt() const + { + throw_error("expected the value to be an integer"); + } + + /// Return the unsigned integer content of the option + virtual unsigned int getUInt() const + { + throw_error("expected the value to be an unsigned integer"); + } + + /// Return the unsigned integer content of the option + virtual bool getBool() const + { + throw_error("expected the value to be a boolean"); + } + + /// Return the vector content of the vector + virtual const std::vector& getVector() const + { + throw_error("expected the value to be a vector"); + } + + /// Return the vector content of the vector + virtual std::vector& getVector() + { + throw_error("expected the value to be a vector"); + } + + /// Return the map content of the vector + virtual const std::map& getDictionary() const + { + throw_error("expected the value to be a dictionary"); + } + + /// Return the map content of the vector + virtual std::map& getDictionary() + { + throw_error("expected the value to be a dictionary"); + } + + /// Return the option content on a path + Option getValue(const std::string& path, Maybe