diff --git a/fvgp/fvgp.py b/fvgp/fvgp.py index a1b2073..beb9896 100755 --- a/fvgp/fvgp.py +++ b/fvgp/fvgp.py @@ -182,8 +182,10 @@ class fvGP(GP): * ``"sparseSolve"`` — direct sparse solve via scipy. * ``"sparseCGpre"`` — preconditioned conjugate-gradient. The preconditioner type is selected by ``args["sparse_preconditioner_type"]`` (default ``"ilu"``; - also ``"ic"``/``"incomplete_cholesky"``, ``"block_jacobi"``, - ``"schwarz"``/``"additive_schwarz"``, or ``"amg"`` (requires pyamg)). + also compiled incomplete Cholesky ``"ichol"``/``"ic"``/``"incomplete_cholesky"``, + zero-fill ``"ichol0"``, legacy Python IC(0) ``"native_ic"``/``"native_ichol"``, + ``"block_jacobi"``, ``"schwarz"``/``"additive_schwarz"``, or ``"amg"`` + (requires pyamg). The compiled IC options require the optional ``ilupp`` package. * ``"sparseMINRESpre"`` — preconditioned MINRES; same preconditioner choices. * ``"sparseCGpre_"`` / ``"sparseMINRESpre_"`` — shortcut that sets ``args["sparse_preconditioner_type"]`` to ```` (e.g. ``"sparseCGpre_amg"``). @@ -239,8 +241,8 @@ class fvGP(GP): - "sparse_krylov_warm_start" : True/False; default = False — feed the previous training iteration's ``KVinvY`` as ``x0`` to the next solve - "sparse_preconditioner_type" : str; default = "ilu". One of "ilu", - "ic"/"ichol"/"incomplete_cholesky", "block_jacobi", "schwarz"/ - "additive_schwarz", "amg" (requires pyamg) + "ichol"/"ic"/"incomplete_cholesky", "ichol0", "native_ic"/"native_ichol", + "block_jacobi", "schwarz"/"additive_schwarz", "amg" (requires pyamg) - "sparse_preconditioner_refresh_interval" : int; default = 1 — reuse the cached preconditioner for up to N consecutive solves before rebuilding. ``set_KV`` always force-refreshes. @@ -250,12 +252,31 @@ class fvGP(GP): additive Schwarz - "sparse_preconditioner_drop_tol" / "sparse_preconditioner_fill_factor" — forwarded to scipy ``spilu`` for "ilu" + - "sparse_preconditioner_ichol_fill_in" / "sparse_preconditioner_ichol_threshold" + — forwarded to ``ilupp`` thresholded IC for "ichol" - "sparse_preconditioner_amg_*" — forwarded to pyamg (``max_levels``, ``max_coarse``, ``strength``, ``cycle``, etc.) - "sparse_preconditioner_shift" / "_growth" / "_attempts" — diagonal - shift retry knobs for "ic" / "block_jacobi" / "additive_schwarz" when + shift retry knobs for legacy "native_ic"/"native_ichol" / "block_jacobi" / "additive_schwarz" when a local Cholesky encounters a non-PD block + Practical sparse-solver guidance: + + - Keep compact-support covariance matrices genuinely sparse before using global + factor-based preconditioners. If the kernel support is too broad, ILU/IC + failures are usually memory/fill failures, not proof that Krylov solvers are bad. + - Prefer ``sparseCGpre_*`` as the primary path for covariance systems, which are + symmetric positive definite in the usual case. ``MINRES`` can be useful as a + comparison, but it may need a stricter ``sparse_minres_tol`` to satisfy a raw + relative residual check. + - Sweep preconditioner parameters at the target scale. For ILU, ``drop_tol`` and + ``fill_factor`` control a memory/solve-time tradeoff; a slightly more expressive + factor can cost only marginally more to build but reduce solve time substantially. + - For repeated nearby K+V updates, enable ``sparse_krylov_warm_start`` and reuse + preconditioners with ``sparse_preconditioner_refresh_interval``. The best refresh + interval is problem-dependent because preconditioner build cost can be comparable + to an unconditioned solve. + Cholesky compute-device routing: - "Chol_factor_compute_device" : str; default = "cpu"/"gpu" diff --git a/fvgp/gp.py b/fvgp/gp.py index fe35f4c..110db08 100755 --- a/fvgp/gp.py +++ b/fvgp/gp.py @@ -183,8 +183,10 @@ class GP: * ``"sparseSolve"`` — direct sparse solve via scipy. * ``"sparseCGpre"`` — preconditioned conjugate-gradient. The preconditioner type is selected by ``args["sparse_preconditioner_type"]`` (default ``"ilu"``; - also ``"ic"``/``"incomplete_cholesky"``, ``"block_jacobi"``, - ``"schwarz"``/``"additive_schwarz"``, or ``"amg"`` (requires pyamg)). + also compiled incomplete Cholesky ``"ichol"``/``"ic"``/``"incomplete_cholesky"``, + zero-fill ``"ichol0"``, legacy Python IC(0) ``"native_ic"``/``"native_ichol"``, + ``"block_jacobi"``, ``"schwarz"``/``"additive_schwarz"``, or ``"amg"`` + (requires pyamg). The compiled IC options require the optional ``ilupp`` package. * ``"sparseMINRESpre"`` — preconditioned MINRES; same preconditioner choices. * ``"sparseCGpre_"`` / ``"sparseMINRESpre_"`` — shortcut that sets ``args["sparse_preconditioner_type"]`` to ```` (e.g. ``"sparseCGpre_amg"``). @@ -240,8 +242,8 @@ class GP: - "sparse_krylov_warm_start" : True/False; default = False — feed the previous training iteration's ``KVinvY`` as ``x0`` to the next solve - "sparse_preconditioner_type" : str; default = "ilu". One of "ilu", - "ic"/"ichol"/"incomplete_cholesky", "block_jacobi", "schwarz"/ - "additive_schwarz", "amg" (requires pyamg) + "ichol"/"ic"/"incomplete_cholesky", "ichol0", "native_ic"/"native_ichol", + "block_jacobi", "schwarz"/"additive_schwarz", "amg" (requires pyamg) - "sparse_preconditioner_refresh_interval" : int; default = 1 — reuse the cached preconditioner for up to N consecutive solves before rebuilding. ``set_KV`` always force-refreshes. @@ -251,12 +253,31 @@ class GP: additive Schwarz - "sparse_preconditioner_drop_tol" / "sparse_preconditioner_fill_factor" — forwarded to scipy ``spilu`` for "ilu" + - "sparse_preconditioner_ichol_fill_in" / "sparse_preconditioner_ichol_threshold" + — forwarded to ``ilupp`` thresholded IC for "ichol" - "sparse_preconditioner_amg_*" — forwarded to pyamg (``max_levels``, ``max_coarse``, ``strength``, ``cycle``, etc.) - "sparse_preconditioner_shift" / "_growth" / "_attempts" — diagonal - shift retry knobs for "ic" / "block_jacobi" / "additive_schwarz" when + shift retry knobs for legacy "native_ic"/"native_ichol" / "block_jacobi" / "additive_schwarz" when a local Cholesky encounters a non-PD block + Practical sparse-solver guidance: + + - Keep compact-support covariance matrices genuinely sparse before using global + factor-based preconditioners. If the kernel support is too broad, ILU/IC + failures are usually memory/fill failures, not proof that Krylov solvers are bad. + - Prefer ``sparseCGpre_*`` as the primary path for covariance systems, which are + symmetric positive definite in the usual case. ``MINRES`` can be useful as a + comparison, but it may need a stricter ``sparse_minres_tol`` to satisfy a raw + relative residual check. + - Sweep preconditioner parameters at the target scale. For ILU, ``drop_tol`` and + ``fill_factor`` control a memory/solve-time tradeoff; a slightly more expressive + factor can cost only marginally more to build but reduce solve time substantially. + - For repeated nearby K+V updates, enable ``sparse_krylov_warm_start`` and reuse + preconditioners with ``sparse_preconditioner_refresh_interval``. The best refresh + interval is problem-dependent because preconditioner build cost can be comparable + to an unconditioned solve. + Cholesky compute-device routing: - "Chol_factor_compute_device" : str; default = "cpu"/"gpu" diff --git a/fvgp/gp_kv.py b/fvgp/gp_kv.py index 392c791..2755737 100755 --- a/fvgp/gp_kv.py +++ b/fvgp/gp_kv.py @@ -49,6 +49,7 @@ def __init__(self, self.Preconditioner_signature = None self.Preconditioner_KV_shape = None self.Preconditioner_reuse_counter = 0 + self.Last_preconditioner_error = None self.allowed_modes = ["Chol", "CholInv", "Inv", "sparseMINRES", "sparseCG", "sparseLU", "sparseMINRESpre", "sparseCGpre", "sparseMINRESpre_", "sparseCGpre_", @@ -121,6 +122,7 @@ def _reset_sparse_preconditioner(self): self.Preconditioner_signature = None self.Preconditioner_KV_shape = None self.Preconditioner_reuse_counter = 0 + self.Last_preconditioner_error = None def _can_reuse_sparse_preconditioner(self, KV): if self.mode not in self._PRECONDITIONED_MODES: @@ -140,9 +142,12 @@ def _build_sparse_preconditioner_or_none(self, KV): try: factor, operator = calculate_sparse_preconditioner(KV, args=self.args) except Exception as exc: + self.Last_preconditioner_error = f"{type(exc).__name__}: {exc}" warnings.warn( f"Failed to build sparse preconditioner for mode {self.mode}; " - "falling back to the unpreconditioned iterative solve." + f"falling back to the unpreconditioned iterative solve. " + f"Reason: {self.Last_preconditioner_error}. " + f"{sparse_preconditioner_failure_guidance(self.args)}" ) logger.warning("Sparse preconditioner construction failed for {}: {}", self.mode, exc) return None, None @@ -172,6 +177,7 @@ def _get_or_refresh_preconditioner(self, KV, force_refresh=False): self.Preconditioner_signature = self._preconditioner_signature() self.Preconditioner_KV_shape = KV.shape self.Preconditioner_reuse_counter = 0 + self.Last_preconditioner_error = None return operator ################################################################## @@ -505,6 +511,7 @@ def __getstate__(self): Preconditioner_signature=self.Preconditioner_signature, Preconditioner_KV_shape=self.Preconditioner_KV_shape, Preconditioner_reuse_counter=self.Preconditioner_reuse_counter, + Last_preconditioner_error=self.Last_preconditioner_error, custom_obj=self.custom_obj, allowed_modes=self.allowed_modes, logdet_KV=self.logdet_KV @@ -520,6 +527,7 @@ def __setstate__(self, state): ("Preconditioner_signature", None), ("Preconditioner_KV_shape", None), ("Preconditioner_reuse_counter", 0), + ("Last_preconditioner_error", None), ): if attr not in self.__dict__: setattr(self, attr, default) diff --git a/fvgp/gp_lin_alg.py b/fvgp/gp_lin_alg.py index e1fa28c..189fd87 100755 --- a/fvgp/gp_lin_alg.py +++ b/fvgp/gp_lin_alg.py @@ -6,7 +6,7 @@ fine-grained options such as ``"GPU_engine"`` and ``"GPU_device"``. Also exposes a sparse preconditioner framework (``calculate_sparse_preconditioner``) -with ILU, incomplete Cholesky (IC(0)), block Jacobi, additive Schwarz, and AMG +with ILU, incomplete Cholesky, block Jacobi, additive Schwarz, and AMG backends, along with a block conjugate-gradient solver for multi-RHS systems. """ import importlib @@ -19,7 +19,7 @@ from scipy import sparse from scipy.linalg import cho_factor, cho_solve, solve_triangular from scipy.sparse import identity -from scipy.sparse.linalg import cg, minres, onenormest, splu, spsolve +from scipy.sparse.linalg import cg, minres, onenormest, splu, spsolve, spsolve_triangular warnings.simplefilter("once", UserWarning) @@ -345,9 +345,16 @@ def normalize_sparse_preconditioner_type(preconditioner_type): preconditioner_type = str(preconditioner_type).lower() aliases = { "ilu": "ilu", - "ic": "incomplete_cholesky", - "ichol": "incomplete_cholesky", - "incomplete_cholesky": "incomplete_cholesky", + "ic": "ichol", + "ichol": "ichol", + "incomplete_cholesky": "ichol", + "ichol0": "ichol0", + "native_ic": "native_incomplete_cholesky", + "native_ichol": "native_incomplete_cholesky", + "legacy_ic": "native_incomplete_cholesky", + "legacy_ichol": "native_incomplete_cholesky", + "native_incomplete_cholesky": "native_incomplete_cholesky", + "legacy_incomplete_cholesky": "native_incomplete_cholesky", "block_jacobi": "block_jacobi", "blockjacobi": "block_jacobi", "schwarz": "additive_schwarz", @@ -358,12 +365,58 @@ def normalize_sparse_preconditioner_type(preconditioner_type): raise ValueError( "Unknown sparse preconditioner type " f"{preconditioner_type!r}. Expected one of " - "{'ilu', 'ic', 'incomplete_cholesky', 'block_jacobi', 'blockjacobi', " + "{'ilu', 'ichol', 'ic', 'ichol0', 'incomplete_cholesky', " + "'native_ic', 'native_ichol', 'legacy_ic', 'legacy_ichol', " + "'block_jacobi', 'blockjacobi', " "'schwarz', 'additive_schwarz', 'amg'}." ) return aliases[preconditioner_type] +def _raise_missing_ilupp(preconditioner_type, exc): + raise ImportError( + "The sparse incomplete-Cholesky preconditioners (`ichol`, `ic`, " + "`incomplete_cholesky`, and `ichol0`) require the optional `ilupp` " + "package. Install it in the Python environment running fvGP with:\n\n" + " pip install ilupp\n\n" + f"Requested sparse preconditioner resolved to backend={preconditioner_type!r}." + ) from exc + + +def sparse_preconditioner_failure_guidance(args=None): + args = _normalize_args(args) + preconditioner_type = args.get("sparse_preconditioner_type") + try: + preconditioner_type = normalize_sparse_preconditioner_type(preconditioner_type) + except Exception: + preconditioner_type = str(preconditioner_type) + + guidance = [ + "Practical guidance: preconditioner failures often mean the covariance graph is too dense or the factor is too expressive for available memory.", + "First check the compact-support kernel length scale/support radius and keep matrix density low before tuning solver parameters.", + "Run a small preconditioner build sweep before a full solve run; a buildable preconditioner can still be slow to apply.", + ] + if preconditioner_type == "ilu": + guidance.append( + "For ILU, sweep `sparse_preconditioner_drop_tol` and `sparse_preconditioner_fill_factor`; looser drop tolerances and smaller fill factors are more likely to fit, while stronger factors may reduce solve time." + ) + elif preconditioner_type in {"ichol", "ichol0"}: + guidance.append( + "For IC/IChol, install the optional backend with `pip install ilupp`; if thresholded IC does not fit, try softer fill/threshold settings or `ichol0`, then verify actual solve time." + ) + elif preconditioner_type in {"block_jacobi", "additive_schwarz"}: + guidance.append( + "For block/local preconditioners, sweep block size and overlap; these may fit easily but can be weaker than ILU on large covariance systems." + ) + guidance.append( + "For repeated nearby K+V updates, `sparse_krylov_warm_start=True` and a nontrivial `sparse_preconditioner_refresh_interval` can avoid rebuilding every solve." + ) + guidance.append( + "If MINRES returns with a poor raw residual, try a stricter `sparse_minres_tol` before judging the method." + ) + return " ".join(guidance) + + def resolve_gp2scale_linalg_mode(mode, args=None): """Split a mode string like ``"sparseCGpre_amg"`` into ``("sparseCGpre", args)``. @@ -546,7 +599,7 @@ def _build_additive_schwarz_preconditioner(KV, args=None): def _build_ic0_factor(KV, args=None): """Pure-Python IC(0) preconditioner. - Correct but slow; intended for moderate problem sizes. Falls back to + Correct but slow; intended only as a legacy/debugging path. Falls back to increasing diagonal shifts if a non-positive pivot is encountered. """ args = _normalize_args(args) @@ -604,38 +657,41 @@ def _build_ic0_factor(KV, args=None): diag[i] = np.sqrt(diag_sq) rows.append(computed_entries) - columns = [[] for _ in range(n)] + indptr = [0] + indices = [] + data = [] for i, row in enumerate(rows): for j, value in row.items(): - columns[j].append((i, value)) + indices.append(int(j)) + data.append(float(value)) + indices.append(i) + data.append(float(diag[i])) + indptr.append(len(indices)) + + L = sparse.csr_matrix( + ( + np.asarray(data, dtype=np.float64), + np.asarray(indices, dtype=np.int32), + np.asarray(indptr, dtype=np.int32), + ), + shape=A.shape, + ) + LT = L.transpose().tocsr() def solve(vector): vector = np.asarray(vector, dtype=np.float64) - y = np.zeros_like(vector) - for i, row in enumerate(rows): - total = vector[i] - for j, value in row.items(): - total -= value * y[j] - y[i] = total / diag[i] - - z = np.zeros_like(y) - for i in range(n - 1, -1, -1): - total = y[i] - for row_index, value in columns[i]: - if row_index > i: - total -= value * z[row_index] - z[i] = total / diag[i] - return z + y = spsolve_triangular(L, vector, lower=True) + return spsolve_triangular(LT, y, lower=False) factor = { - "type": "incomplete_cholesky", + "type": "native_incomplete_cholesky", "diag": diag, "rows": rows, - "columns": columns, - "solve": solve, + "L": L, + "LT": LT, "shift": shift, } - operator = sparse.linalg.LinearOperator(A.shape, matvec=solve, rmatvec=solve, dtype=A.dtype) + operator = _build_dtype_adapted_operator(A.shape, solve, factor_dtype=np.float64) return factor, operator except Exception as exc: last_exc = exc @@ -644,6 +700,30 @@ def solve(vector): raise np.linalg.LinAlgError(f"IC(0) preconditioner construction failed after shifted retries: {last_exc}") +def _build_dtype_adapted_operator(shape, solve, factor_dtype, operator_dtype=np.float64): + factor_dtype = np.dtype(factor_dtype) + operator_dtype = np.dtype(operator_dtype) + + def _apply(vec): + arr = np.asarray(vec, dtype=operator_dtype) + if arr.ndim == 1: + solved = solve(np.asarray(arr, dtype=factor_dtype)) + return np.asarray(solved, dtype=operator_dtype) + columns = [ + np.asarray(solve(np.asarray(arr[:, i], dtype=factor_dtype)), dtype=operator_dtype) + for i in range(arr.shape[1]) + ] + return np.column_stack(columns) + + return sparse.linalg.LinearOperator( + shape, + matvec=_apply, + rmatvec=_apply, + matmat=_apply, + dtype=operator_dtype, + ) + + def _build_ilu_preconditioner(KV, args=None): args = _normalize_args(args) A = KV.tocsc() @@ -659,12 +739,34 @@ def _build_ilu_preconditioner(KV, args=None): spilu_kwargs["diag_pivot_thresh"] = args["sparse_preconditioner_diag_pivot_thresh"] factor = sparse.linalg.spilu(A, **spilu_kwargs) - operator = sparse.linalg.LinearOperator( - A.shape, - matvec=factor.solve, - rmatvec=factor.solve, - dtype=A.dtype, - ) + operator = _build_dtype_adapted_operator(A.shape, factor.solve, factor_dtype=A.dtype) + return factor, operator + + +def _build_ichol0_preconditioner(KV, args=None): + try: + import ilupp + except ImportError as exc: + _raise_missing_ilupp("ichol0", exc) + + A = _as_symmetric_csr(KV).astype(np.float64) + factor = ilupp.IChol0Preconditioner(A) + operator = _build_dtype_adapted_operator(A.shape, factor.dot, factor_dtype=np.float64) + return factor, operator + + +def _build_ichol_preconditioner(KV, args=None): + args = _normalize_args(args) + try: + import ilupp + except ImportError as exc: + _raise_missing_ilupp("ichol", exc) + + A = _as_symmetric_csr(KV).astype(np.float64) + add_fill_in = int(args.get("sparse_preconditioner_ichol_fill_in", 16)) + threshold = float(args.get("sparse_preconditioner_ichol_threshold", 1e-4)) + factor = ilupp.ICholTPreconditioner(A, add_fill_in=add_fill_in, threshold=threshold) + operator = _build_dtype_adapted_operator(A.shape, factor.dot, factor_dtype=np.float64) return factor, operator @@ -712,7 +814,9 @@ def calculate_sparse_preconditioner(KV, args=None): builders = { "ilu": _build_ilu_preconditioner, - "incomplete_cholesky": _build_ic0_factor, + "native_incomplete_cholesky": _build_ic0_factor, + "ichol0": _build_ichol0_preconditioner, + "ichol": _build_ichol_preconditioner, "block_jacobi": _build_block_jacobi_preconditioner, "additive_schwarz": _build_additive_schwarz_preconditioner, "amg": _build_amg_preconditioner, diff --git a/tests/test_fvgp.py b/tests/test_fvgp.py index 536eb1a..2064202 100755 --- a/tests/test_fvgp.py +++ b/tests/test_fvgp.py @@ -385,15 +385,20 @@ def test_linalg_modes(): # the GP must function end-to-end (posterior + data updates) with that # preconditioner backing the iterative solve. canonical_to_aliases = { - "sparseCGpre": ["ilu", "ic", "block_jacobi", "schwarz"], - "sparseMINRESpre": ["ilu", "ic", "block_jacobi", "schwarz"], + "sparseCGpre": ["ilu", "native_ic", "block_jacobi", "schwarz"], + "sparseMINRESpre": ["ilu", "native_ic", "block_jacobi", "schwarz"], } + if _il.util.find_spec("ilupp") is not None: + canonical_to_aliases["sparseCGpre"].extend(["ichol", "ichol0"]) + canonical_to_aliases["sparseMINRESpre"].extend(["ichol", "ichol0"]) if _il.util.find_spec("pyamg") is not None: canonical_to_aliases["sparseCGpre"].append("amg") canonical_to_aliases["sparseMINRESpre"].append("amg") canonical_to_type = { "ilu": "ilu", - "ic": "incomplete_cholesky", + "ichol": "ichol", + "ichol0": "ichol0", + "native_ic": "native_incomplete_cholesky", "block_jacobi": "block_jacobi", "schwarz": "additive_schwarz", "amg": "amg", @@ -1356,8 +1361,11 @@ def _make_test_spd_sparse(n=40, seed=0): def test_normalize_sparse_preconditioner_type(): assert normalize_sparse_preconditioner_type("ILU") == "ilu" - assert normalize_sparse_preconditioner_type("ic") == "incomplete_cholesky" - assert normalize_sparse_preconditioner_type("ichol") == "incomplete_cholesky" + assert normalize_sparse_preconditioner_type("ic") == "ichol" + assert normalize_sparse_preconditioner_type("ichol") == "ichol" + assert normalize_sparse_preconditioner_type("ichol0") == "ichol0" + assert normalize_sparse_preconditioner_type("native_ic") == "native_incomplete_cholesky" + assert normalize_sparse_preconditioner_type("native_ichol") == "native_incomplete_cholesky" assert normalize_sparse_preconditioner_type("BlockJacobi") == "block_jacobi" assert normalize_sparse_preconditioner_type("schwarz") == "additive_schwarz" assert normalize_sparse_preconditioner_type("AMG") == "amg" @@ -1374,8 +1382,17 @@ def test_resolve_gp2scale_linalg_mode(): mode, args = resolve_gp2scale_linalg_mode("sparseCGpre_amg") assert mode == "sparseCGpre" and args["sparse_preconditioner_type"] == "amg" + mode, args = resolve_gp2scale_linalg_mode("sparseMINRESpre_ichol") + assert mode == "sparseMINRESpre" and args["sparse_preconditioner_type"] == "ichol" + mode, args = resolve_gp2scale_linalg_mode("sparseMINRESpre_ic") - assert mode == "sparseMINRESpre" and args["sparse_preconditioner_type"] == "incomplete_cholesky" + assert mode == "sparseMINRESpre" and args["sparse_preconditioner_type"] == "ichol" + + mode, args = resolve_gp2scale_linalg_mode("sparseMINRESpre_native_ic") + assert mode == "sparseMINRESpre" and args["sparse_preconditioner_type"] == "native_incomplete_cholesky" + + mode, args = resolve_gp2scale_linalg_mode("sparseCGpre_native_ichol") + assert mode == "sparseCGpre" and args["sparse_preconditioner_type"] == "native_incomplete_cholesky" # Consistent explicit type is allowed mode, args = resolve_gp2scale_linalg_mode( @@ -1401,16 +1418,36 @@ def test_calculate_sparse_preconditioner_ilu(): assert res < 1e-6 -def test_calculate_sparse_preconditioner_ic0(): +def test_calculate_sparse_preconditioner_native_ic0(): A = _make_test_spd_sparse(n=30) - factor, op = calculate_sparse_preconditioner(A, args={"sparse_preconditioner_type": "incomplete_cholesky"}) - assert factor["type"] == "incomplete_cholesky" + factor, op = calculate_sparse_preconditioner(A, args={"sparse_preconditioner_type": "native_ic"}) + assert factor["type"] == "native_incomplete_cholesky" b = np.random.rand(A.shape[0]) x = calculate_sparse_conj_grad(A, b, M=op, args={"sparse_cg_tol": 1e-8}) res = np.linalg.norm(A @ x[:, 0] - b) / np.linalg.norm(b) assert res < 1e-6 +def test_missing_ilupp_message_for_ic_aliases(monkeypatch): + import builtins + + real_import = builtins.__import__ + + def fake_import(name, *args, **kwargs): + if name == "ilupp": + raise ImportError("simulated missing ilupp") + return real_import(name, *args, **kwargs) + + monkeypatch.setattr(builtins, "__import__", fake_import) + A = _make_test_spd_sparse(n=10) + + with pytest.raises(ImportError, match="pip install ilupp"): + calculate_sparse_preconditioner(A, args={"sparse_preconditioner_type": "ic"}) + + with pytest.raises(ImportError, match="pip install ilupp"): + calculate_sparse_preconditioner(A, args={"sparse_preconditioner_type": "ichol0"}) + + def test_calculate_sparse_preconditioner_block_jacobi(): A = _make_test_spd_sparse(n=30) factor, op = calculate_sparse_preconditioner( @@ -1761,7 +1798,7 @@ def test_kv_preconditioner_signature_invalidates_cache(): assert op0 is not None # Mutate args to flip the preconditioner type - gp.data.args["sparse_preconditioner_type"] = "ic" + gp.data.args["sparse_preconditioner_type"] = "native_ic" KV = kv.addKV(kv.K, kv.V) kv.update_KV(KV) assert kv.Preconditioner_operator is not None @@ -1864,5 +1901,3 @@ def test_preconditioner_build_failure_falls_back(): # And the build-failure warning fired msgs = [str(w.message) for w in caught if issubclass(w.category, UserWarning)] assert any("Failed to build sparse preconditioner" in m for m in msgs) - -