py-econometrics · schroedk · May 11, 2026 · May 11, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -54,6 +54,12 @@ Modified LSMR is now the sole iterative solver, replacing CG and GMRES.
   `FePreconditioner::additive_schwarz_diagnostics`, and the Python
   `AdditiveSchwarzDiagnostics` class. Scheduling metrics are now private
   to the `Auto` heuristic (closes #34).
+- **BREAKING:** `Operator::apply` / `apply_adjoint` now return
+  `Result<(), SolveError>`; the infallible pair and the `try_apply` /
+  `try_apply_adjoint` defaults are gone. `ApplyError` is removed; its
+  variants moved onto `SolveError`. `PyFePreconditioner.apply` raises
+  `RuntimeError` on local-solver failure instead of returning NaNs
+  (closes #29).
 
 ## [0.1.0] - 2026-03-12
 

diff --git a/crates/schwarz-precond/examples/additive_schwarz.rs b/crates/schwarz-precond/examples/additive_schwarz.rs
@@ -4,8 +4,8 @@
 //! subdomains and diagonal local solvers.
 
 use schwarz_precond::{
-    lsmr, mlsmr, LocalSolveError, LocalSolver, Operator, SchwarzPreconditioner, SubdomainCore,
-    SubdomainEntry,
+    lsmr, mlsmr, LocalSolveError, LocalSolver, Operator, SchwarzPreconditioner, SolveError,
+    SubdomainCore, SubdomainEntry,
 };
 
 // ---------------------------------------------------------------------------
@@ -23,7 +23,7 @@ impl Operator for TridiagOperator {
     fn ncols(&self) -> usize {
         self.n
     }
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
         for i in 0..self.n {
             y[i] = 3.0 * x[i];
             if i > 0 {
@@ -33,12 +33,12 @@ impl Operator for TridiagOperator {
                 y[i] -= x[i + 1];
             }
         }
+        Ok(())
     }
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
-        self.apply(x, y); // symmetric
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
+        self.apply(x, y) // symmetric
     }
 }
-
 // ---------------------------------------------------------------------------
 // Diagonal local solver: y = rhs / diag_val
 // ---------------------------------------------------------------------------

diff --git a/crates/schwarz-precond/examples/custom_local_solver.rs b/crates/schwarz-precond/examples/custom_local_solver.rs
@@ -6,10 +6,9 @@
 use faer::{MatRef, Side};
 
 use schwarz_precond::{
-    mlsmr, LocalSolveError, LocalSolver, Operator, SchwarzPreconditioner, SparseMatrix,
+    mlsmr, LocalSolveError, LocalSolver, Operator, SchwarzPreconditioner, SolveError, SparseMatrix,
     SubdomainCore, SubdomainEntry,
 };
-
 // ---------------------------------------------------------------------------
 // Tridiagonal operator
 // ---------------------------------------------------------------------------
@@ -25,7 +24,7 @@ impl Operator for TridiagOperator {
     fn ncols(&self) -> usize {
         self.n
     }
-    fn apply(&self, x: &[f64], y: &mut [f64]) {
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
         for i in 0..self.n {
             y[i] = 3.0 * x[i];
             if i > 0 {
@@ -35,12 +34,12 @@ impl Operator for TridiagOperator {
                 y[i] -= x[i + 1];
             }
         }
+        Ok(())
     }
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) {
-        self.apply(x, y);
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), SolveError> {
+        self.apply(x, y)
     }
 }
-
 // ---------------------------------------------------------------------------
 // Build the same tridiag as a SparseMatrix (CSR) for submatrix extraction
 // ---------------------------------------------------------------------------

diff --git a/crates/schwarz-precond/src/error.rs b/crates/schwarz-precond/src/error.rs
@@ -1,14 +1,13 @@
 //! Error types for the `schwarz-precond` crate.
 //!
-//! Errors are layered to match the build → apply → solve lifecycle:
+//! Errors are layered to match the build → solve lifecycle:
 //!
 //! - **Build errors** ([`SubdomainCoreBuildError`], [`SubdomainEntryBuildError`],
 //!   [`PreconditionerBuildError`]) — caught during construction, before any
 //!   solve begins.
-//! - **Apply errors** ([`ApplyError`]) — runtime failures during a single
-//!   preconditioner or operator application (e.g. a local solver diverges).
-//! - **Solve errors** ([`SolveError`]) — wraps `ApplyError` for the iterative
-//!   solver layer.
+//! - **Solve errors** ([`SolveError`]) — runtime failures during a solve,
+//!   including operator/preconditioner application (e.g. a local solver
+//!   diverges) and iterative-solver input validation.
 //!
 //! Each level chains to its source via [`Error::source`].
 
@@ -145,51 +144,25 @@ impl Display for LocalSolveError {
 
 impl Error for LocalSolveError {}
 
-/// Runtime failure while applying a Schwarz preconditioner/operator.
+/// Runtime failure while executing a solve.
+///
+/// Covers both operator/preconditioner application failures (e.g. a local
+/// solver diverges) and iterative-solver input validation.
 #[derive(Debug, Clone, PartialEq, Eq)]
-pub enum ApplyError {
-    /// A local subdomain solve failed.
+#[non_exhaustive]
+pub enum SolveError {
+    /// A local subdomain solve failed during a preconditioner apply.
     LocalSolveFailed {
         /// Index of the failing subdomain entry in the preconditioner.
         subdomain: usize,
         /// Local solver error.
         source: LocalSolveError,
     },
-    /// Internal synchronization failed (e.g. poisoned mutex).
+    /// Internal synchronization failed (e.g. poisoned mutex) during an apply.
     Synchronization {
         /// Context string identifying the lock/synchronization site.
         context: &'static str,
     },
-}
-
-impl Display for ApplyError {
-    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
-        match self {
-            Self::LocalSolveFailed { subdomain, source } => {
-                write!(f, "subdomain {subdomain} local solve failed: {source}")
-            }
-            Self::Synchronization { context } => {
-                write!(f, "synchronization failure at {context}")
-            }
-        }
-    }
-}
-
-impl Error for ApplyError {
-    fn source(&self) -> Option<&(dyn Error + 'static)> {
-        match self {
-            Self::LocalSolveFailed { source, .. } => Some(source),
-            Self::Synchronization { .. } => None,
-        }
-    }
-}
-
-/// Runtime failure while executing an iterative solver.
-#[derive(Debug, Clone, PartialEq, Eq)]
-#[non_exhaustive]
-pub enum SolveError {
-    /// Operator/preconditioner apply failed.
-    Apply(ApplyError),
     /// Solver input was invalid before any iteration was attempted.
     InvalidInput {
         /// Context string identifying the validation site.
@@ -202,7 +175,12 @@ pub enum SolveError {
 impl Display for SolveError {
     fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
         match self {
-            Self::Apply(err) => write!(f, "operator apply failed: {err}"),
+            Self::LocalSolveFailed { subdomain, source } => {
+                write!(f, "subdomain {subdomain} local solve failed: {source}")
+            }
+            Self::Synchronization { context } => {
+                write!(f, "synchronization failure at {context}")
+            }
             Self::InvalidInput { context, message } => {
                 write!(f, "invalid solver input at {context}: {message}")
             }
@@ -213,18 +191,13 @@ impl Display for SolveError {
 impl Error for SolveError {
     fn source(&self) -> Option<&(dyn Error + 'static)> {
         match self {
-            Self::Apply(err) => Some(err),
+            Self::LocalSolveFailed { source, .. } => Some(source),
+            Self::Synchronization { .. } => None,
             Self::InvalidInput { .. } => None,
         }
     }
 }
 
-impl From<ApplyError> for SolveError {
-    fn from(value: ApplyError) -> Self {
-        Self::Apply(value)
-    }
-}
-
 pub(crate) fn validate_entries<S: LocalSolver>(
     entries: &[SubdomainEntry<S>],
     n_dofs: usize,

diff --git a/crates/schwarz-precond/src/lib.rs b/crates/schwarz-precond/src/lib.rs
@@ -92,31 +92,22 @@
 ///
 /// Preconditioners are operators too (M^{-1} is a linear map).
 /// All implementors must be Send + Sync to enable Rayon parallelism.
+///
+/// Both apply methods are fallible: implementors that cannot fail in practice
+/// (matrices, identity operators) still return `Result<(), SolveError>` so
+/// callers can use a uniform `?` propagation path. Symmetric operators
+/// should delegate `apply_adjoint` to `apply`.
 pub trait Operator: Send + Sync {
     /// Number of rows in the operator.
     fn nrows(&self) -> usize;
     /// Number of columns in the operator.
     fn ncols(&self) -> usize;
-    /// y = A*x
-    fn apply(&self, x: &[f64], y: &mut [f64]);
+    /// Computes y = A * x. Returns an error if the apply fails at runtime
+    /// (e.g. a local subdomain solver diverges).
+    fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), error::SolveError>;
     /// Computes y = A^T * x. For symmetric operators, this should delegate to `apply`.
-    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]);
-
-    /// Fallible version of [`Operator::apply`].
-    ///
-    /// Implementors with runtime failure modes should override this.
-    fn try_apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), error::ApplyError> {
-        self.apply(x, y);
-        Ok(())
-    }
-
-    /// Fallible version of [`Operator::apply_adjoint`].
-    ///
-    /// Implementors with runtime failure modes should override this.
-    fn try_apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), error::ApplyError> {
-        self.apply_adjoint(x, y);
-        Ok(())
-    }
+    /// Returns an error under the same conditions as `apply`.
+    fn apply_adjoint(&self, x: &[f64], y: &mut [f64]) -> Result<(), error::SolveError>;
 }
 
 // ============================================================================
@@ -133,7 +124,7 @@ mod sparse_matrix;
 
 pub use domain::{PartitionWeights, SubdomainCore};
 pub use error::{
-    ApplyError, LocalSolveError, PreconditionerBuildError, SolveError, SubdomainCoreBuildError,
+    LocalSolveError, PreconditionerBuildError, SolveError, SubdomainCoreBuildError,
     SubdomainEntryBuildError,
 };
 pub use local_solve::{LocalSolver, SubdomainEntry};

diff --git a/crates/schwarz-precond/src/lsmr/bidiag.rs b/crates/schwarz-precond/src/lsmr/bidiag.rs
@@ -278,7 +278,7 @@ pub(super) trait Bidiagonalization {
 impl<A: Operator + ?Sized> Bidiagonalization for GolubKahan<'_, A> {
     fn step(&mut self) -> Result<BidiagStep, SolveError> {
         // β_{k+1} u_{k+1} = A v_k − α_k u_k
-        self.operator.try_apply(&self.bufs.v, &mut self.bufs.av)?;
+        self.operator.apply(&self.bufs.v, &mut self.bufs.av)?;
         let beta_sq = axpy_with_sq_norm(&mut self.bufs.u, &self.bufs.av, -self.alpha);
         let beta = beta_sq.sqrt();
         if beta == 0.0 {
@@ -295,7 +295,7 @@ impl<A: Operator + ?Sized> Bidiagonalization for GolubKahan<'_, A> {
 
         // α_{k+1} v_{k+1} = Aᵀ u_{k+1} − β_{k+1} v_k
         self.operator
-            .try_apply_adjoint(&self.bufs.u, &mut self.bufs.atu)?;
+            .apply_adjoint(&self.bufs.u, &mut self.bufs.atu)?;
         let mut alpha_sq = axpy_with_sq_norm(&mut self.bufs.v, &self.bufs.atu, -beta);
 
         // Windowed MGS in the standard inner product, before normalization.
@@ -330,7 +330,7 @@ impl<A: Operator + ?Sized, M: Operator + ?Sized> Bidiagonalization
         // Phase 1 — ũ_{k+1} unnormalized: scale = −α_k / β_k.
         // Compute ũ_{k+1} = A ṽ_k − (α_k / β_k) ũ_k, then β_{k+1} = ‖ũ_{k+1}‖.
         let scale = -(self.alpha * self.beta_prev_inv);
-        self.operator.try_apply(&self.bufs.v, &mut self.bufs.av)?;
+        self.operator.apply(&self.bufs.v, &mut self.bufs.av)?;
         let beta_sq = axpy_with_sq_norm(&mut self.bufs.u, &self.bufs.av, scale);
         let beta = beta_sq.sqrt();
         if beta == 0.0 {
@@ -353,7 +353,7 @@ impl<A: Operator + ?Sized, M: Operator + ?Sized> Bidiagonalization
         // Precondition: α_k > 0 (the outer loop guard in lsmr_from_bidiag
         // never calls step() once α has collapsed to zero).
         self.operator
-            .try_apply_adjoint(&self.bufs.u, &mut self.bufs.atu)?;
+            .apply_adjoint(&self.bufs.u, &mut self.bufs.atu)?;
         debug_assert!(
             self.alpha > 0.0,
             "self.alpha must be > 0; lsmr_from_bidiag's loop guard prevents step() after alpha=0",
@@ -371,7 +371,7 @@ impl<A: Operator + ?Sized, M: Operator + ?Sized> Bidiagonalization
         // step), and push the normalized v_{k+1} and p̃_{k+1,norm} = p_tilde /
         // α_new into the ring for the next step's MGS.
         self.preconditioner
-            .try_apply(&self.bufs.p_tilde, &mut self.bufs.v)?;
+            .apply(&self.bufs.p_tilde, &mut self.bufs.v)?;
 
         if let Some(reorth) = &self.bufs.local_reorth {
             reorth.reorthogonalize(&mut self.bufs.v, &mut self.bufs.p_tilde);
@@ -464,7 +464,7 @@ impl<'a, A: Operator + ?Sized> GolubKahan<'a, A> {
         }
 
         // α₁ v₁ = Aᵀ u₁
-        operator.try_apply_adjoint(&bufs.u, &mut bufs.v)?;
+        operator.apply_adjoint(&bufs.u, &mut bufs.v)?;
         let alpha = par_dot(&bufs.v, &bufs.v).sqrt();
         if alpha > 0.0 {
             scale_in_place(&mut bufs.v, 1.0 / alpha);
@@ -555,10 +555,10 @@ impl<'a, A: Operator + ?Sized, M: Operator + ?Sized> ModifiedGolubKahan<'a, A, M
         }
 
         // p̃ = Aᵀ u₁
-        operator.try_apply_adjoint(&bufs.u, &mut bufs.p_tilde)?;
+        operator.apply_adjoint(&bufs.u, &mut bufs.p_tilde)?;
 
         // ṽ₁ = M⁻¹ p̃
-        preconditioner.try_apply(&bufs.p_tilde, &mut bufs.v)?;
+        preconditioner.apply(&bufs.p_tilde, &mut bufs.v)?;
 
         // α₁ = √⟨ṽ₁, p̃⟩ via the M-norm dot product trick.
         let vp = par_dot(&bufs.v, &bufs.p_tilde);