num::kernel::subspace Namespace Reference

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

Functions
template<typename F >
CallableOp< F >	make_op (F f, idx n)
	Factory: wrap a callable as a stack-allocated CallableOp<F>.
template<typename F >
std::unique_ptr< LinearOp >	make_op_ptr (F f, idx n)
	Factory: wrap a callable as a heap-allocated LinearOp.
real	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	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	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.

Function Documentation

arnoldi_step()

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.

Algorithm:

1. w = A * basis[k] (operator apply via scratch buffer)
2. h[0..k] = mgs_orthogonalize(basis[0..k], w)
3. h[k+1] = ||w|| (subdiagonal Hessenberg element)
4. if h[k+1] > breakdown_tol: w /= h[k+1]; basis.push_back(w) (new orthonormal basis vector)

Parameters

A	Linear operator
basis	Current orthonormal basis (basis[0..k] present on entry; basis[k+1] appended on exit if no breakdown)
h	Hessenberg column k; must be pre-allocated with size >= k+2
k	Current step index (0-based)
scratch	Pre-allocated Vector of size n; reused across steps to avoid per-step heap allocation
breakdown_tol	Lucky breakdown threshold (default 1e-14)

Returns

h[k+1]: the subdiagonal element. Near-zero signals lucky breakdown (A-invariant subspace found; solution is exact).

Definition at line 73 of file subspace.cpp.

References num::kernel::subspace::LinearOp::apply(), num::beta(), mgs_orthogonalize(), and num::scale().

Referenced by num::gmres().

make_op()

template<typename F >

CallableOp< F > num::kernel::subspace::make_op	(	F	f,
		idx	n
	)

Factory: wrap a callable as a stack-allocated CallableOp<F>.

Use when the operator type is known at the call site (zero overhead; compiler can inline apply()).

Definition at line 128 of file subspace.hpp.

Referenced by num::gmres().

make_op_ptr()

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.

Use when you need type-erased storage (unique_ptr<LinearOp>).

Definition at line 136 of file subspace.hpp.

mgs_orthogonalize() [1/2]

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.

Overload for algorithms (e.g. Lanczos) that store their Krylov basis column-by-column in a dense Matrix rather than a std::vector<Vector>. The Matrix is stored row-major, so column l is accessed with stride basis.cols() between successive elements.

When NUMERICS_HAS_BLAS is defined, uses cblas_ddot / cblas_daxpy with the correct stride so the BLAS handles the strided access natively. Otherwise falls back to a plain loop via basis(i, l).

Parameters

basis	Row-major Matrix of shape (n, max_steps); columns are the orthonormal basis vectors.
k	Number of columns to project out (columns 0..k-1).
v	Vector to orthogonalize in-place (length n).

Returns

||v||₂ after all projections are removed.

Definition at line 47 of file subspace.cpp.

References num::Matrix::cols(), num::Matrix::data(), num::BasicVector< T >::data(), num::norm(), and num::Matrix::rows().

mgs_orthogonalize() [2/2]

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].

For each i in [0, k), computes h[i] = <v, basis[i]> then subtracts h[i] * basis[i] from v. Updates v in-place.

h must be pre-allocated with size >= k.

Returns ||v||₂ after all projections are removed (= the next subdiagonal element in the Hessenberg matrix when used inside arnoldi_step).

MGS performs the projections sequentially; each one modifies v before the next dot product, giving better numerical stability than classical GS at the cost of k sequential passes over v.

Uses num::dot and num::axpy with default_backend (BLAS > SIMD > seq).

Definition at line 32 of file subspace.cpp.

References num::axpy(), num::dot(), and num::norm().

Referenced by arnoldi_step(), num::expv(), and num::lanczos().