Linear Solvers

Linear Solvers

Goal

Provide a universal workflow to select a solver, assess conditioning, and diagnose convergence for linear systems arising in numerical simulations.

Requirements

• Python 3.10+
• NumPy, SciPy (for matrix operations)
• See individual scripts for dependencies

Inputs to Gather

| Input | Description | Example | |-------|-------------|---------| | Matrix size | Dimension of system | n = 1000000 | | Sparsity | Fraction of nonzeros | 0.01% | | Symmetry | Is A = Aᵀ? | yes | | Definiteness | Is A positive definite? | yes (SPD) | | Conditioning | Estimated condition number | 10⁶ |

Decision Guidance

Solver Selection Flowchart

Is matrix dense and small enough to factor in memory (dense float64
storage n²·8 bytes  0.95),
because stagnation is a tail property — early fast drops can hide a flat tail.
Read `asymptotic_rate` when judging the regime below:

| Asymptotic rate | Meaning | Action |
|-----------------|---------|--------|
|  0.95 | Stagnation | Change solver or preconditioner |

### Stagnation Diagnosis

| Pattern | Likely Cause | Fix |
|---------|--------------|-----|
| Flat residual | Poor preconditioner | Improve preconditioner |
| Oscillating | Near-singular or indefinite | Check matrix, try different solver |
| Very slow decay | Ill-conditioned | Apply scaling, use AMG |

## Verification checklist

Do not trust a solve until each of these is satisfied with a recorded value, not a "looks fine":

- [ ] Recorded `asymptotic_rate` from `convergence_diagnostics.py` and confirmed it is below the 0.95 stagnation threshold (and ideally < 0.5); a low whole-history `rate` alone does not rule out a flat tail.
- [ ] Checked the relative residual from `residual_norms.py` against the physics-scaled `--rel-tol` (default 1e-6), not just the absolute norm; for unscaled RHS use `--require-both` so an undersized `rhs` cannot fake convergence.
- [ ] Confirmed `solver_selector.py` `recommended` matches the actual matrix properties recorded from `sparsity_stats.py` (`symmetry`, and definiteness if known) — e.g. CG only when symmetric AND positive-definite, MINRES for symmetric-indefinite, GMRES/BiCGSTAB for nonsymmetric.
- [ ] For systems flagged ill-conditioned, ran `scaling_equilibration.py` and recorded `row_scale_max/row_scale_min` and `col_scale_max/col_scale_min`; for symmetric matrices used `--symmetric` (D A D) so symmetry is preserved, and applied row_scale THEN col_scale for nonsymmetric two-sided scaling.
- [ ] Reviewed `sparsity_stats.py` `notes`/`zero_rows`/`zero_cols` from `scaling_equilibration.py` — any zero row or column means the system is structurally singular and the scale-of-1 fallback is not a fix.
- [ ] Confirmed the preconditioner from `preconditioner_advisor.py` is admissible for the chosen Krylov method — in particular a MINRES preconditioner must be SPD (an indefinite incomplete LDLᵀ is invalid).

## Common pitfalls & rationalizations

| Tempting shortcut | Why it's wrong / what to do |
|-------------------|-----------------------------|
| "The mean `rate` is low, so it converged." | `rate` is the whole-history mean and is dominated by early fast drops; stagnation is a tail property. Read `asymptotic_rate` and confirm it is below 0.95. |
| "The absolute residual is tiny, so we're done." | A small absolute norm can be meaningless if the RHS is large or unscaled. Check the `relative_norms` / `relative_value` against a physics-scaled `--rel-tol`. |
| "It's symmetric, so just use CG." | CG requires symmetric AND positive-definite. A symmetric-indefinite matrix needs MINRES (with an SPD preconditioner); using CG can break down or stall. Confirm definiteness before selecting. |
| "Large system, so factor it directly." | `solver_selector.py` gates dense direct solvers on dense float64 storage (n²·8 bytes < ~2 GB, n ≈ 16384); above that a dense Cholesky/LU is infeasible and you must route to an iterative method. |
| "Scaling is just dividing each row by its max." | One-sided row scaling does not equilibrate. For nonsymmetric matrices derive `col_scale` from the row-scaled matrix and apply both; for symmetric matrices use the symmetric D A D scale or you destroy symmetry. |
| "GMRES stagnates, so add more iterations." | A flat tail means the preconditioner or restart length is the problem, not iteration count. Strengthen the preconditioner (higher ILU fill / AMG), increase the restart parameter, or switch methods. |

## Security

### Input Validation
- All numeric inputs (residuals, tolerances, matrix entries) are validated as finite numbers
- Comma-separated residual/vector inputs are capped at 100,000 entries
- The `solver_selector.py` `--size` parameter is bounded at 10 billion
- `--matrix-type` is validated against a fixed allowlist (`spd`, `symmetric-indefinite`, `nonsymmetric`)
- Boolean flags (`--symmetric`, `--positive-definite`, `--sparse`, `--ill-conditioned`) are type-safe argparse flags

### File Access
- `sparsity_stats.py` and `scaling_equilibration.py` read a single matrix file (`.npy` format) specified by `--matrix`
- `np.load()` is called with `allow_pickle=False` to prevent arbitrary code execution via crafted `.npy` files
- Matrix files are rejected if they exceed 500 MB before any parsing occurs
- Matrix dimension limits (100,000 per dimension) prevent memory exhaustion
- All other scripts read no external files; inputs are provided via CLI arguments

### Tool Restrictions
- **Read**: Used to inspect script source, references, and matrix files
- **Write**: Used to save analysis results or solver recommendations; writes are scoped to the user's working directory
- **Grep/Glob**: Used to locate relevant files and search references
- The skill's `allowed-tools` excludes `Bash` to prevent the agent from executing arbitrary commands when processing untrusted matrix files or numeric inputs

### Safety Measures
- No `eval()`, `exec()`, or dynamic code generation
- All subprocess calls use explicit argument lists (no `shell=True`)
- Reduced tool surface (no Bash) limits the agent to read/write operations only
- JSON output mode produces structured, parseable results without shell-interpretable content

## Limitations

- **Large dense matrices**: Direct solvers may run out of memory
- **Highly indefinite**: Standard preconditioners may fail
- **Saddle-point**: Requires specialized block preconditioners

## References

- `references/solver_decision_tree.md` - Selection logic
- `references/preconditioner_catalog.md` - Preconditioner options
- `references/convergence_patterns.md` - Diagnosing failures
- `references/scaling_guidelines.md` - Equilibration guidance

## Version History

- **v1.2.0** (2026-06-23): Fixed asymptotic stagnation detection, dense-feasibility solver gating, saddle-point/small-dense direct-solver routing, equilibrating two-sided scaling, CG iteration-bound table, and doc/eval consistency
- **v1.1.0** (2024-12-24): Enhanced documentation, decision guidance, examples
- **v1.0.0**: Initial release with 6 solver analysis scripts

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [HeshamFS](https://github.com/HeshamFS)
- **Source:** [HeshamFS/materials-simulation-skills](https://github.com/HeshamFS/materials-simulation-skills)
- **License:** Apache-2.0

Install and usage instructions live in the source repository linked above.