util.h

#ifndef LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_UTIL_H
#define LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_UTIL_H
#include <molpro/Profiler.h>
#include <molpro/linalg/itsolv/ArrayHandlers.h>
#include <molpro/linalg/itsolv/subspace/Matrix.h>
#include <molpro/linalg/itsolv/wrap.h>
 
#include <cmath>
#include <ios>
#include <locale>
#include <map>
#include <ostream>
#include <string>
#include <type_traits>
#include <vector>
 
namespace molpro::linalg::itsolv::util {
 
template <class R>
auto norm(array::ArrayHandler<R, R>& handler, const std::type_identity_t<R>& x) ->
    typename array::ArrayHandler<R, R>::value_type_abs {
  return std::sqrt(std::abs(handler.dot(x, x)));
}
 
template <typename value_type, typename value_type_abs>
void remove_null_vectors(subspace::Matrix<value_type>& lin_trans, std::vector<value_type_abs>& norm, const size_t start,
                         const size_t end, const value_type_abs norm_thresh) {
  for (size_t i = start, j = start; i < end; ++i) {
    if (norm[j] < norm_thresh) {
      norm.erase(begin(norm) + j);
      lin_trans.remove_row(j);
    } else
      ++j;
  }
}
 
template <class ForwardIt, class EndIterator, typename Int>
bool is_iota(ForwardIt it_start, EndIterator it_end, Int value_start) {
  bool result = true;
  for (auto it = it_start; it != it_end && result; ++it, ++value_start)
    result = (*it) == value_start;
  return result;
}
 
template <class R, class Q>
Q construct_zeroed_copy(const R& param, array::ArrayHandler<Q, R>& handler) {
  auto q = handler.copy(param);
  handler.fill(0, q);
  return q;
}
 
template <class R, class Q, class P, typename value_type>
void construct_solutions(const VecRef<R>& params, const std::vector<int>& roots,
                         const subspace::Matrix<value_type>& solutions, const CVecRef<P>& pparams,
                         const CVecRef<Q>& qparams, const CVecRef<Q>& dparams, size_t oP, size_t oQ, size_t oD,
                         array::ArrayHandler<R, R>& handler_rr, array::ArrayHandler<R, P>& handler_rp,
                         array::ArrayHandler<R, Q>& handler_rq) {
  auto prof = molpro::Profiler::single();
  prof->start("itsolv::util::construct_solutions");
  assert(params.size() >= roots.size());
  const auto nroots = roots.size();
  for (size_t i = 0; i < nroots; ++i)
    handler_rr.fill(0, params.at(i));
  // coefficient matrices with the subspace vectors at the given indices in rows and roots in columns, as required by
  // gemm_outer
  auto coefficients = [&](const std::vector<size_t>& indices) {
    auto coeff = subspace::Matrix<value_type>({indices.size(), nroots});
    for (size_t i = 0; i < nroots; ++i)
      for (size_t j = 0; j < indices.size(); ++j)
        coeff(j, i) = solutions(roots[i], indices[j]);
    return coeff;
  };
  auto append_indices = [](std::vector<size_t>& indices, size_t offset, size_t n) {
    for (size_t j = 0; j < n; ++j)
      indices.push_back(offset + j);
  };
  if (nroots > 0) {
    const auto wparams = VecRef<R>(params.begin(), params.begin() + nroots);
    auto p_indices = std::vector<size_t>{};
    append_indices(p_indices, oP, pparams.size());
    handler_rp.gemm_outer(coefficients(p_indices), pparams, wparams);
    // Q and D in one call, so that each solution is passed over only once
    auto qd_indices = std::vector<size_t>{};
    append_indices(qd_indices, oQ, qparams.size());
    append_indices(qd_indices, oD, dparams.size());
    auto qdparams = qparams;
    qdparams.insert(qdparams.end(), dparams.begin(), dparams.end());
    handler_rq.gemm_outer(coefficients(qd_indices), qdparams, wparams);
  }
  prof->stop();
}
 
template <class Container>
void delete_parameters(std::vector<int> indices, Container& params) {
  std::sort(std::begin(indices), std::end(indices), std::greater<int>());
  for (auto i : indices)
    params.erase(begin(params) + i);
}
 
class StringFacet {
public:
  std::string toupper(std::string in) const;
  std::string tolower(std::string in) const;
  bool tobool(const std::string& in) const;
  static void crop_space(std::string& path);
  static std::map<std::string, std::string> parse_keyval_string(std::string s);
 
private:
  const std::ctype<char>& facet = std::use_facet<std::ctype<char>>(std::locale());
};
 
std::size_t tosize(const std::string& key, const std::string& value);
 
class StreamFormatGuard {
public:
  explicit StreamFormatGuard(std::ostream& stream)
      : m_stream(stream), m_flags(stream.flags()), m_precision(stream.precision()) {}
  ~StreamFormatGuard() {
    m_stream.flags(m_flags);
    m_stream.precision(m_precision);
  }
  StreamFormatGuard(const StreamFormatGuard&) = delete;
  StreamFormatGuard& operator=(const StreamFormatGuard&) = delete;
 
private:
  std::ostream& m_stream;
  std::ios_base::fmtflags m_flags;
  std::streamsize m_precision;
};
 
} // namespace molpro::linalg::itsolv::util
#endif // LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_UTIL_H