subspace.hpp File Reference

Subspace construction and orthogonalization kernels. (namespace num::kernel::subspace) More...

#include "core/matrix.hpp"
#include "core/policy.hpp"
#include "core/types.hpp"
#include "core/vector.hpp"
#include "linalg/sparse/sparse.hpp"
#include <memory>
#include <vector>

Go to the source code of this file.

Classes
struct	num::kernel::subspace::LinearOp
	Abstract matrix-free linear operator: y = A*x. More...
struct	num::kernel::subspace::DenseOp
	Wrap a dense Matrix as a LinearOp. More...
struct	num::kernel::subspace::SparseOp
	Wrap a SparseMatrix as a LinearOp. More...
struct	num::kernel::subspace::CallableOp< F >
	Wrap any callable void(const Vector&, Vector&) as a LinearOp. More...

Namespaces
namespace	num
namespace	num::kernel
namespace	num::kernel::subspace

Functions
template<typename F >
CallableOp< F >	num::kernel::subspace::make_op (F f, idx n)
	Factory: wrap a callable as a stack-allocated CallableOp<F>.
template<typename F >
std::unique_ptr< LinearOp >	num::kernel::subspace::make_op_ptr (F f, idx n)
	Factory: wrap a callable as a heap-allocated LinearOp.
real	num::kernel::subspace::mgs_orthogonalize (const std::vector< Vector > &basis, Vector &v, std::vector< real > &h, idx k)
	Modified Gram-Schmidt: orthogonalize v against basis[0..k-1].
real	num::kernel::subspace::mgs_orthogonalize (const Matrix &basis, idx k, Vector &v)
	Modified Gram-Schmidt: orthogonalize v against columns 0..k-1 of a column-stored basis matrix.
real	num::kernel::subspace::arnoldi_step (const LinearOp &A, std::vector< Vector > &basis, std::vector< real > &h, idx k, Vector &scratch, real breakdown_tol=real(1e-14))
	One Arnoldi step: expand the orthonormal basis by one vector.

Detailed Description

Subspace construction and orthogonalization kernels. (namespace num::kernel::subspace)

These are the inner loops that any algorithm building or maintaining an orthonormal basis reduces to – Krylov solvers (GMRES, MINRES), Lanczos, rational Krylov, exponential integrators (expv), randomized SVD, POD, and model-order reduction all share the same core operations.

Linear operator interface: LinearOp – abstract base: apply(x, y), rows(), cols() DenseOp – wraps Matrix (Backend-dispatched matvec) SparseOp – wraps SparseMatrix CallableOp<F> – wraps any callable void(const Vector&, Vector&) make_op(f, n) – factory returning CallableOp<F> (stack, typed) make_op_ptr(f, n) – factory returning unique_ptr<LinearOp> (heap, erased)

Orthogonalization kernels: mgs_orthogonalize – Modified Gram-Schmidt: orthogonalize one vector against an existing basis, sequential and stable. arnoldi_step – One Arnoldi step: apply operator + MGS + normalize. The inner loop that GMRES, Lanczos, and expv share.

Key design choice – scratch buffer: arnoldi_step takes a pre-allocated scratch Vector for the operator apply. Callers allocate it once before the loop; the kernel reuses it at every step. This avoids one heap allocation per Arnoldi step.

Key design choice – no OMP inside subspace kernels: mgs_orthogonalize and arnoldi_step call num::dot / num::axpy / num::norm with default_backend (which is BLAS > SIMD > blocked). The outer basis- expansion loop is inherently serial (each step depends on the previous). OMP at this level would only cause nested-parallelism issues.

Definition in file subspace.hpp.