KrylovScope

krylovscope estimates extremal Ritz values of large symmetric operators with a GPU-resident Lanczos process. It is intended for solver setup, spectral conditioning diagnostics, preconditioner diagnostics, domain-decomposition setup, and applications that need spectral endpoint estimates without assembling or solving a dense eigensystem.

Version 0.1.0 provides the validated CUDA path.

日本語の導入・運用・診断手順は GUIDE.ja.md を参照してください。

CUDA algorithm

• resident q / q_prev / q_next / z: f64
• basis storage: f32
• deterministic radix-16 scalar reductions: f64
• tile-2 mixed Gram: f32 basis × f64 vector -> f64
• one-pass CGS reorthogonalization
• device-side breakdown/non-finite latch
• scale-relative breakdown detection and scaled tridiagonal extraction
• zero scalar readbacks during recurrence
• one packed coefficient readback at the end
• optional free-DOF mask for constrained systems
• explicit prepare() to isolate first-use CubeCL/CUDA compilation

The design was developed through the E.36-E.38 validation series and has been exercised with a real matrix-free Tet10 elasticity operator at 363,216 DOF.

Important numerical semantics

RitzBounds contains the extremal Ritz values of the finite Krylov projection. Despite the type name, these values are estimates, not guaranteed enclosing bounds for the spectrum of the original operator.

In particular, the lower endpoint can converge more slowly and can be more seed-dependent than the upper endpoint. For difficult spectra, use multiple seeds with estimate_ensemble_prepared() and inspect min_relative_spread, max_relative_spread, tail_beta, and the number of completed steps.

RitzBounds::condition_estimate() is therefore an estimated condition number, not a rigorous condition-number bound. If a Ritz endpoint is used to choose a stability-critical polynomial or time-step interval, apply an application-level safety margin or an independent enclosure method.

Using a CSR matrix

use krylovscope::cubecl::cuda::{CudaDevice, CudaRuntime};
use krylovscope::{SparseSession, SparseSessionConfig, SparseView};

# fn demo(row_ptr: &[usize], col_idx: &[usize], values: &[f64]) -> Result<(), String> {
let device = CudaDevice::default();
let mut spectral = SparseSession::<CudaRuntime>::new(
    &device,
    SparseView::CsrF64Usize {
        row_ptr,
        col_idx,
        values,
    },
    SparseSessionConfig::default(),
)?;

spectral.prepare(0x1234)?;
let estimate = spectral.estimate_prepared(0x5678)?;
let bounds = estimate.require_bounds()?;
println!("lambda = [{}, {}]", bounds.lambda_min, bounds.lambda_max);
# Ok(()) }

Enable CUDA:

[dependencies]
krylovscope = { version = "0.1", features = ["cuda"] }

Matrix-free application operators

Applications can implement ResidentF64Operator<R> and use SpectralEstimator<R, O> directly. The application operator launches resident y = A*x; Krylov vectors, basis storage, reductions, preparation, reorthogonalization, and final Ritz extraction remain owned by krylovscope. apply() must overwrite every element of y; the output buffer is not promised to be zero-initialized.

For constrained FEM, construct the estimator with SpectralEstimator::new_masked. A nonzero mask entry denotes a free DOF. The library projects both the start vector and every operator output so the Krylov recurrence remains in the selected subspace even when the underlying operator does not apply the constraint mask itself.

Prepared execution

The first CubeCL/CUDA use can include substantial compilation and initialization cost. Production callers should normally call prepare(seed) once, then use estimate_prepared(seed) for subsequent estimates. prepare() is idempotent.

Current limitations

• CUDA is the production backend in 0.1.0.
• The validated public estimator is capped at 32 Lanczos steps (MAX_STEPS).
• The operator is expected to represent a symmetric problem; sparse sessions can validate structure/symmetry before upload.
• Ritz endpoints are finite-step estimates, not certified spectral enclosures.
• Basis storage is f32 and reorthogonalization is one-pass CGS; difficult clustered spectra should be checked with multiple seeds and residual diagnostics.

Validation

The release line has been validated against CPU references, deterministic breakdown/non-finite tests, CSR/PDSM operator paths, prepared/reuse semantics, free-DOF masking, and a production-size matrix-free Tet10 structural operator. Scale-invariance, arbitrary start normalization, masked starts, CSR safety validation, and symmetry-policy edge cases are covered by regression tests.

The final E.38E integration gate used 121,072 nodes, 78,004 Tet10 elements, 363,216 DOF, and 5,151 constrained DOF. The external-crate integration, operator comparison, masked subspace, two 32-step spectral estimates, and resident synchronization checks all passed.

License

MIT. See LICENSE.