matrix_free.hpp

Go to the documentation of this file.

#pragma once
 
#include <uni20/core/scalar_concepts.hpp>
#include <uni20/core/math.hpp>
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <concepts>
#include <cstddef>
#include <limits>
#include <optional>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <vector>
 
namespace uni20::krylov
{
 
enum class SpectrumPart
{
  LargestMagnitude,
  SmallestMagnitude,
  LargestAlgebraic,
  SmallestAlgebraic,
  BothEnds,
  LargestReal,
  SmallestReal,
  LargestImaginary,
  SmallestImaginary
};
 
enum class KrylovDiagnosticsLevel
{
  None,
  Summary,
  Full
};
 
enum class RealNonsymmetricPolicy
{
  RequireRealEigenpairs,
  AllowRealSchurPairs,
  PromoteToComplexSuggested
};
 
enum class NonsymmetricStatus
{
  Converged,
  NotConverged,
  ComplexPairEncountered,
  ComplexPromotionRecommended,
  RealSchurPairRequired,
  AmbiguousReality,
  Breakdown,
  InvalidInput
};
 
enum class RitzReality
{
  Real,
  Ambiguous,
  Complex
};
 
enum class SymmetricTransformKind
{
  Regular,
  GeneralizedInverse,
  ShiftInvert,
  Buckling,
  Cayley
};
 
using SymmetricSpectralTransform = SymmetricTransformKind;
 
template <uni20::Real Scalar> struct SymmetricTransformOptions
{
    SymmetricTransformKind transform = SymmetricTransformKind::Regular;
    Scalar sigma = Scalar{};
};
 
template <uni20::Real Scalar> using SymmetricSpectralTransformOptions = SymmetricTransformOptions<Scalar>;
 
template <uni20::Real Scalar> struct NonsymmetricEigenParams
{
    int eigenvalue_count = 1;
    int retained_ritz_count = 0;
    int krylov_dimension = 0;
    int max_iterations = 300;
    Scalar tolerance = Scalar{};
    Scalar complex_pair_tolerance = Scalar{};
    SpectrumPart spectrum = SpectrumPart::LargestMagnitude;
    bool compute_eigenvectors = true;
    RealNonsymmetricPolicy real_policy = RealNonsymmetricPolicy::RequireRealEigenpairs;
    KrylovDiagnosticsLevel diagnostics = KrylovDiagnosticsLevel::None;
};
 
template <uni20::Real Scalar> struct NonsymmetricArnoldiDiagnostics
{
    int restart_count = 0;
    int op_count = 0;
    int final_projected_dimension = 0;
    Scalar final_residual_norm = Scalar{};
    Scalar projected_departure_from_normality = Scalar{};
    std::vector<uni20::complex<Scalar>> final_ritz_values;
    std::vector<Scalar> final_ritz_bounds;
    std::vector<RitzReality> final_ritz_reality;
};
 
template <uni20::Real Scalar, typename Vector> struct NonsymmetricEigenResult
{
    std::vector<uni20::complex<Scalar>> eigenvalues;
    std::vector<Scalar> residual_bounds;
    std::vector<Vector> right_eigenvectors;
    std::vector<RitzReality> reality;
    int converged_count = 0;
    int iteration_count = 0;
    int matvec_count = 0;
    NonsymmetricStatus status = NonsymmetricStatus::Converged;
    std::optional<NonsymmetricArnoldiDiagnostics<Scalar>> diagnostics;
};
 
template <uni20::Real Scalar> struct SymmetricEigenParams
{
    int eigenvalue_count = 1;
    int retained_ritz_count = 0;
    int krylov_dimension = 0;
    int max_iterations = 300;
    Scalar tolerance = Scalar{};
    Scalar breakdown_tolerance = Scalar{};
    SpectrumPart spectrum = SpectrumPart::LargestMagnitude;
    bool compute_eigenvectors = true;
    KrylovDiagnosticsLevel diagnostics = KrylovDiagnosticsLevel::None;
};
 
[[nodiscard]] inline int default_symmetric_krylov_dimension(int eigenvalue_count, std::size_t problem_dimension)
{
  if (eigenvalue_count <= 0)
  {
    throw std::invalid_argument("default symmetric Krylov dimension requires a positive eigenvalue count");
  }
  if (problem_dimension > static_cast<std::size_t>(uni20::numeric_limits<int>::max()))
  {
    throw std::overflow_error("Krylov problem dimension exceeds the supported int range");
  }
  int const dimension = static_cast<int>(problem_dimension);
  int const target = std::max(20, 2 * eigenvalue_count + 1);
  return std::min(dimension, target);
}
 
[[nodiscard]] inline int default_nonsymmetric_krylov_dimension(int eigenvalue_count, std::size_t problem_dimension)
{
  if (eigenvalue_count <= 0)
  {
    throw std::invalid_argument("default nonsymmetric Krylov dimension requires a positive eigenvalue count");
  }
  if (problem_dimension > static_cast<std::size_t>(uni20::numeric_limits<int>::max()))
  {
    throw std::overflow_error("Krylov problem dimension exceeds the supported int range");
  }
  int const dimension = static_cast<int>(problem_dimension);
  int const target = std::max(20, 6 * eigenvalue_count + 8);
  return std::min(dimension, target);
}
 
template <uni20::Real Scalar>
[[nodiscard]] int effective_symmetric_krylov_dimension(SymmetricEigenParams<Scalar> const& params,
                                                       std::size_t problem_dimension)
{
  return params.krylov_dimension > 0 ? params.krylov_dimension
                                     : default_symmetric_krylov_dimension(params.eigenvalue_count, problem_dimension);
}
 
template <uni20::Real Scalar>
[[nodiscard]] int effective_nonsymmetric_krylov_dimension(NonsymmetricEigenParams<Scalar> const& params,
                                                          std::size_t problem_dimension)
{
  return params.krylov_dimension > 0
             ? params.krylov_dimension
             : default_nonsymmetric_krylov_dimension(params.eigenvalue_count, problem_dimension);
}
 
[[nodiscard]] inline int default_symmetric_retained_ritz_count(int eigenvalue_count, int krylov_dimension)
{
  if (eigenvalue_count <= 0)
  {
    throw std::invalid_argument("default symmetric retained Ritz count requires a positive eigenvalue count");
  }
  if (krylov_dimension <= eigenvalue_count)
  {
    return eigenvalue_count;
  }
  return eigenvalue_count;
}
 
[[nodiscard]] inline int default_nonsymmetric_retained_ritz_count(int eigenvalue_count, int krylov_dimension)
{
  if (eigenvalue_count <= 0)
  {
    throw std::invalid_argument("default nonsymmetric retained Ritz count requires a positive eigenvalue count");
  }
  if (krylov_dimension <= eigenvalue_count)
  {
    return eigenvalue_count;
  }
  int const target = std::max(2 * eigenvalue_count, eigenvalue_count + 4);
  return std::min(krylov_dimension - 1, target);
}
 
template <uni20::Real Scalar>
[[nodiscard]] int effective_symmetric_retained_ritz_count(SymmetricEigenParams<Scalar> const& params,
                                                          int krylov_dimension)
{
  int const retained_count = params.retained_ritz_count > 0
                                 ? params.retained_ritz_count
                                 : default_symmetric_retained_ritz_count(params.eigenvalue_count, krylov_dimension);
  if (retained_count < params.eigenvalue_count)
  {
    throw std::invalid_argument("retained Ritz count must be at least the requested eigenvalue count");
  }
  if (retained_count >= krylov_dimension)
  {
    throw std::invalid_argument("retained Ritz count must be smaller than the Krylov dimension");
  }
  return retained_count;
}
 
template <uni20::Real Scalar>
[[nodiscard]] int effective_nonsymmetric_retained_ritz_count(NonsymmetricEigenParams<Scalar> const& params,
                                                             int krylov_dimension)
{
  int const retained_count = params.retained_ritz_count > 0
                                 ? params.retained_ritz_count
                                 : default_nonsymmetric_retained_ritz_count(params.eigenvalue_count, krylov_dimension);
  if (retained_count < params.eigenvalue_count)
  {
    throw std::invalid_argument("nonsymmetric retained Ritz count is below eigenvalue count");
  }
  if (retained_count >= krylov_dimension)
  {
    throw std::invalid_argument("nonsymmetric retained Ritz count must be smaller than the Krylov dimension");
  }
  return retained_count;
}
 
template <uni20::Real Scalar> struct SymmetricRestartCycleDiagnostics
{
    int cycle = 0;
    int converged_count = 0;
    std::size_t projected_dimension = 0;
    std::size_t retained_count = 0;
    std::size_t shift_count = 0;
    std::size_t protected_zero_bound_count = 0;
    bool enlarged_for_partial_convergence = false;
    Scalar residual_norm = Scalar{};
    std::vector<Scalar> ritz_values;
    std::vector<Scalar> ritz_bounds;
    std::vector<Scalar> shifts;
    std::vector<std::size_t> wanted_indices;
    std::vector<std::size_t> shift_indices;
    std::vector<bool> wanted_converged;
};
 
template <uni20::Real Scalar> struct SymmetricLanczosDiagnostics
{
    int restart_count = 0;
    int op_count = 0;
    int final_projected_dimension = 0;
    Scalar final_residual_norm = Scalar{};
    std::vector<Scalar> final_ritz_values;
    std::vector<Scalar> final_ritz_bounds;
    std::vector<std::size_t> final_selected_indices;
    std::vector<SymmetricRestartCycleDiagnostics<Scalar>> restart_cycles;
};
 
template <typename Scalar, typename Vector> struct EigenResult
{
    std::vector<Scalar> eigenvalues;
    std::vector<Scalar> residual_bounds;
    std::vector<Vector> eigenvectors;
    int converged_count = 0;
    int iteration_count = 0;
    int matvec_count = 0;
    int status = 0;
    std::optional<SymmetricLanczosDiagnostics<Scalar>> diagnostics;
};
 
template <typename Ops, typename Vector, typename Scalar>
concept KrylovVectorOps = requires(Ops ops, Vector x, Vector y, Scalar alpha) {
  { ops.problem_dimension() } -> std::convertible_to<std::size_t>;
  { ops.vector_dimension(x) } -> std::convertible_to<std::size_t>;
  { ops.allocate_like(x) } -> std::same_as<Vector>;
  { ops.copy(y, x) } -> std::same_as<void>;
  { ops.axpy(y, alpha, x) } -> std::same_as<void>;
  { ops.scal(y, alpha) } -> std::same_as<void>;
  { ops.set_zero(y) } -> std::same_as<void>;
  { ops.inner_product(x, y) } -> std::convertible_to<Scalar>;
};
 
template <typename Op, typename Vector>
concept KrylovOperator = requires(Op op, Vector x, Vector y) {
  { op.matvec(y, x) } -> std::same_as<void>;
};
 
template <typename Ops, typename Vector, typename Scalar>
concept KrylovMatrixFreeOperator = KrylovVectorOps<Ops, Vector, Scalar> && KrylovOperator<Ops, Vector>;
 
template <typename Ops, typename Vector, typename Scalar>
concept MatrixFreeVectorOps = KrylovMatrixFreeOperator<Ops, Vector, Scalar>;
 
template <typename Ops, typename Vector, typename Scalar>
concept KrylovMetricInnerProduct = requires(Ops ops, Vector x, Vector y) {
  { ops.metric_inner_product(x, y) } -> std::convertible_to<Scalar>;
};
 
template <typename Ops, typename Vector, typename Scalar>
concept MetricInnerProductOps = KrylovMetricInnerProduct<Ops, Vector, Scalar>;
 
template <typename Ops, typename Vector, typename Scalar>
concept KrylovNorm = requires(Ops ops, Vector x) {
  { ops.norm(x) } -> std::convertible_to<uni20::make_real_t<Scalar>>;
};
 
template <typename Scalar, typename Vector, typename Ops>
uni20::make_real_t<Scalar> norm_or_inner_product(Ops& ops, Vector const& x)
{
  using RealScalar = uni20::make_real_t<Scalar>;
  if constexpr (KrylovNorm<Ops, Vector, Scalar>)
  {
    return static_cast<RealScalar>(ops.norm(x));
  }
  else
  {
    Scalar const norm_squared = ops.inner_product(x, x);
    RealScalar const real_norm_squared = static_cast<RealScalar>(uni20::math::real(norm_squared));
    return real_norm_squared > RealScalar{} ? uni20::math::sqrt(real_norm_squared) : RealScalar{};
  }
}
 
template <uni20::Real Scalar>
  requires requires { uni20::numeric_limits<Scalar>::epsilon(); }
Scalar default_complex_pair_tolerance()
{
  return uni20::math::sqrt(uni20::numeric_limits<Scalar>::epsilon());
}
 
template <uni20::Real Scalar> Scalar default_complex_pair_tolerance(Scalar const& exemplar)
{
  return uni20::math::sqrt(uni20::numeric_limits<Scalar>::epsilon(exemplar));
}
 
template <uni20::Complex C>
RitzReality classify_ritz_reality(C theta, uni20::make_real_t<C> tolerance = uni20::make_real_t<C>{},
                                  uni20::make_real_t<C> scale = uni20::make_real_t<C>{1},
                                  uni20::make_real_t<C> ambiguity_factor = uni20::make_real_t<C>{10})
{
  using Scalar = uni20::make_real_t<C>;
  Scalar const effective_tolerance = tolerance > Scalar{} ? tolerance : default_complex_pair_tolerance(theta.real());
  if constexpr (requires {
                  theta.is_exact();
                  theta.at(effective_tolerance.precision());
                })
  {
    if (theta.is_exact()) theta = theta.at(effective_tolerance.precision());
  }
  Scalar const effective_scale = std::max({Scalar{1}, uni20::math::abs(theta), uni20::math::abs(scale)});
  Scalar const threshold = effective_tolerance * effective_scale;
  Scalar const imaginary = uni20::math::abs(theta.imag());
  if (imaginary <= threshold)
  {
    return RitzReality::Real;
  }
  if (imaginary <= ambiguity_factor * threshold)
  {
    return RitzReality::Ambiguous;
  }
  return RitzReality::Complex;
}
 
template <typename Scalar, typename Vector, typename MetricOps, typename BOps, typename InnerOps>
Scalar metric_inner_product_or_apply(MetricOps& metric_ops, BOps& b_ops, InnerOps& inner_ops, Vector const& x,
                                     Vector const& y, Vector& scratch)
{
  if constexpr (MetricInnerProductOps<MetricOps, Vector, Scalar>)
  {
    return metric_ops.metric_inner_product(x, y);
  }
  else
  {
    b_ops.matvec(scratch, y);
    return inner_ops.inner_product(x, scratch);
  }
}
 
template <typename Ops, typename Vector>
std::size_t validate_matrix_free_dimensions(Ops const& ops, Vector const& vector, char const* context)
{
  std::size_t const problem_dimension = ops.problem_dimension();
  std::size_t const vector_dimension = ops.vector_dimension(vector);
  if (problem_dimension == 0)
  {
    throw std::invalid_argument(std::string(context) + " requires a non-empty problem dimension");
  }
  if (vector_dimension != problem_dimension)
  {
    throw std::invalid_argument(std::string(context) + " vector dimension does not match operator dimension");
  }
  return problem_dimension;
}
 
} // namespace uni20::krylov