gram_schmidt.h

#ifndef LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_SUBSPACE_GRAM_SCHMIDT_H
#define LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_SUBSPACE_GRAM_SCHMIDT_H
#include <molpro/linalg/array/ArrayHandler.h>
#include <molpro/linalg/itsolv/helper.h>
#include <molpro/linalg/itsolv/subspace/Matrix.h>
#include <molpro/linalg/itsolv/wrap.h>
 
namespace molpro::linalg::itsolv::subspace::util {
template <typename T>
std::vector<T> gram_schmidt(const Matrix<T>& s, Matrix<T>& l, T norm_thresh = precision_scaled<T>(1e-14)) {
  assert(s.rows() == s.cols());
  auto n = s.rows();
  l.fill(0);
  l.resize({n, n});
  auto norm = std::vector<T>(n, 0);
  auto w = std::vector<T>{};
  for (size_t i = 0; i < n; ++i) {
    w.assign(i, 0.);
    for (size_t j = 0; j < i; ++j) {
      for (size_t k = 0; k <= j; ++k) {
        w[j] += s(i, k) * l(j, k);
      }
    }
    for (size_t j = 0; j < i; ++j) {
      if (norm[j] > norm_thresh) {
        for (size_t k = 0; k <= j; ++k) {
          l(i, k) -= w[j] / norm[j] * l(j, k);
        }
      }
    }
    l(i, i) = 1.;
    for (size_t j = 0; j <= i; ++j) {
      for (size_t k = 0; k < j; ++k) {
        norm[i] += 2 * l(i, j) * l(i, k) * s(j, k);
      }
      norm[i] += l(i, j) * l(i, j) * s(j, j);
    }
  }
  std::transform(begin(norm), end(norm), begin(norm), [](auto el) { return std::sqrt(std::abs(el)); });
  return norm;
}
 
// FIXME This is implicitly constructed in Gram Schmidt and we can make it an opitonal return variable
template <typename value_type, typename value_type_abs>
Matrix<value_type> construct_lin_trans_in_orthogonal_set(const Matrix<value_type>& overlap,
                                                         const Matrix<value_type>& lin_trans,
                                                         const std::vector<value_type_abs>& norm) {
  const auto nrows = lin_trans.rows();
  const auto ncols = lin_trans.cols();
  auto t = Matrix<value_type>({nrows, ncols});
  for (size_t i = 0; i < nrows; ++i) {
    for (size_t j = 0; j < i; ++j) {
      for (size_t k = 0; k <= j; ++k) {
        t(i, j) -= lin_trans(j, k) * overlap(i, k) / std::pow(norm[j], 2);
      }
    }
  }
  return t;
}
 
template <class R, typename value_type_abs>
auto modified_gram_schmidt(VecRef<R>& params, array::ArrayHandler<R, R>& handler, value_type_abs null_thresh) {
  auto null_param_indices = std::vector<size_t>{};
  const size_t n = params.size();
  for (size_t i = 0; i < n; ++i) {
    auto norm = handler.dot(params[i], params[i]);
    norm = std::sqrt(std::abs(norm));
    if (norm > null_thresh) {
      handler.scal(1. / norm, params[i]);
      for (size_t j = i + 1; j < n; ++j) {
        auto ov = handler.dot(params[i], params[j]);
        handler.axpy(-ov, params[i], params[j]);
      }
    } else {
      null_param_indices.emplace_back(i);
    }
  }
  return null_param_indices;
}
 
} // namespace molpro::linalg::itsolv::subspace::util
 
#endif // LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_SUBSPACE_GRAM_SCHMIDT_H