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.8+
• NumPy, SciPy (for matrix operations)
• See individual scripts for dependencies

Inputs to Gather

Decision Guidance

Solver Selection Flowchart

Is matrix small (n < 5000) and dense?
├── YES → Use direct solver (LU, Cholesky)
└── NO → Is matrix symmetric?
    ├── YES → Is it positive definite?
    │   ├── YES → Use CG with AMG/IC preconditioner
    │   └── NO → Use MINRES
    └── NO → Is it nearly symmetric?
        ├── YES → Use BiCGSTAB
        └── NO → Use GMRES with ILU/AMG

Quick Reference

Script Outputs (JSON Fields)

Workflow

1. Characterize matrix - symmetry, definiteness, sparsity 2. Analyze sparsity - Run scripts/sparsity_stats.py 3. Select solver - Run scripts/solver_selector.py 4. Choose preconditioner - Run scripts/preconditioner_advisor.py 5. Apply scaling - If ill-conditioned, use scripts/scaling_equilibration.py 6. Monitor convergence - Use scripts/convergence_diagnostics.py 7. Diagnose issues - Check residual history with scripts/residual_norms.py

Conversational Workflow Example

User: My GMRES solver is stagnating after 50 iterations. The residual drops to 1e-3 then stops improving.

Agent workflow: 1. Diagnose convergence:

   python3 scripts/convergence_diagnostics.py --residuals 1,0.1,0.01,0.005,0.003,0.002,0.002,0.002 --json

2. Check for preconditioning advice:

   python3 scripts/preconditioner_advisor.py --matrix-type nonsymmetric --sparse --stagnation --json

3. Recommend: Increase restart parameter, try ILU(k) with higher k, or switch to AMG.

Pre-Solve Checklist

• [ ] Confirm matrix symmetry/definiteness
• [ ] Decide direct vs iterative based on size and sparsity
• [ ] Set residual tolerance relative to physics scale
• [ ] Choose preconditioner appropriate to matrix structure
• [ ] Apply scaling/equilibration if needed
• [ ] Track convergence and adjust if stagnation occurs

CLI Examples

# Analyze sparsity pattern
python3 scripts/sparsity_stats.py --matrix A.npy --json

# Select solver for SPD sparse system
python3 scripts/solver_selector.py --symmetric --positive-definite --sparse --size 1000000 --json

# Get preconditioner recommendation
python3 scripts/preconditioner_advisor.py --matrix-type spd --sparse --json

# Diagnose convergence from residual history
python3 scripts/convergence_diagnostics.py --residuals 1,0.2,0.05,0.01 --json

# Apply scaling
python3 scripts/scaling_equilibration.py --matrix A.npy --symmetric --json

# Compute residual norms
python3 scripts/residual_norms.py --residual 1,0.1,0.01 --rhs 1,0,0 --json

Error Handling

Interpretation Guidance

Convergence Rate

Stagnation Diagnosis

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.1.0 (2024-12-24): Enhanced documentation, decision guidance, examples
• v1.0.0: Initial release with 6 solver analysis scripts

Convergence Patterns

Comprehensive guide for diagnosing and improving iterative solver convergence.

Convergence Monitoring

Key Metrics

Convergence Rate

Asymptotic convergence rate ρ = lim_{k→∞} ||r_k|| / ||r_{k-1}||

Linear convergence: ||r_k|| ≤ ρ^k ||r_0||, ρ < 1
Superlinear: ρ_k → 0 as k → ∞

Typical Convergence Patterns

Pattern 1: Fast Monotonic Decay

Residual:
1e0   ****
1e-2       ****
1e-4            ****
1e-6                 ****
      |----|----|----|----|
      0    10   20   30   40  iterations

Characteristics:

• Smooth, consistent reduction
• Constant or decreasing rate
• Reaches tolerance quickly

Indicates:

• Good preconditioner match
• Well-conditioned problem
• Appropriate solver choice

Action: None needed, this is ideal.

Pattern 2: Initial Stall Then Decay

Residual:
1e0   ****----
1e-2          \____
1e-4               \____
1e-6                    \****
      |----|----|----|----|
      0    10   20   30   40  iterations

Characteristics:

• Slow start, then acceleration
• Superlinear convergence later
• Common with GMRES

Indicates:

• Krylov subspace building useful information
• Eventually finds good direction

Action:

• Be patient
• Increase max iterations
• Don't restart GMRES too early

Pattern 3: Stagnation

Residual:
1e0   ****
1e-2       ****--------------------
1e-4
1e-6
      |----|----|----|----|----|----|
      0    20   40   60   80   100  iterations

Characteristics:

• Residual stops decreasing
• Convergence rate ≈ 1
• May plateau at various levels

Indicates:

• Preconditioner inadequate
• Matrix too ill-conditioned
• Reached limits of floating-point precision

Actions:

Pattern 4: Oscillation

Residual:
1e0   *   *   *
1e-2   * * * * *
1e-4    *   *   *
1e-6
      |----|----|----|----|
      0    10   20   30   40  iterations

Characteristics:

• Residual bounces up and down
• Net progress may be slow
• Common with BiCGSTAB

Indicates:

• Near-singular or indefinite matrix
• Eigenvalue close to origin
• Loss of orthogonality (GMRES)

Actions:

• Try different solver (GMRES if using BiCGSTAB)
• Increase GMRES restart parameter
• Better preconditioner
• Check matrix for issues

Pattern 5: Divergence

Residual:
1e0   ****
1e2        ****
1e4             ****
NaN                  ****
      |----|----|----|----|
      0    10   20   30   40  iterations

Characteristics:

• Residual increases each iteration
• Eventually overflow/NaN

Indicates:

• Matrix singular or nearly so
• Instability in method
• Severe ill-conditioning

Actions:

• Check matrix (is it singular?)
• Apply scaling
• Try direct solver
• Regularize if appropriate

Pattern 6: Plateau with Breakthrough

Residual:
1e0   ****
1e-2       ****--------****
1e-4                        \____
1e-6                             ****
      |----|----|----|----|----|----|
      0    20   40   60   80   100  iterations

Characteristics:

• Long plateau, then sudden progress
• Multiple plateaus possible

Indicates:

• Multiple scales in problem
• Different eigenvalue clusters

Action:

• Increase max iterations
• Consider multigrid or multilevel preconditioner

Solver-Specific Patterns

CG (Conjugate Gradient)

Expected behavior:

• Monotonic residual decrease
• At most n iterations for exact arithmetic
• Affected by eigenvalue distribution

Warning signs:

GMRES

Expected behavior:

• Monotonically decreasing residual (in exact arithmetic)
• May stall then accelerate
• Memory grows with iterations

Restart effects:

BiCGSTAB

Expected behavior:

• May oscillate
• Two matrix-vector products per iteration
• Residual can temporarily increase

Warning signs:

Diagnostic Procedures

Step 1: Basic Health Check

1. Verify matrix properties:
   - Is it supposed to be SPD? Check.
   - Approximate condition number?
   - Check for zero rows/columns.

2. Check RHS:
   - Is ||b|| reasonable?
   - No NaN or Inf?

3. Verify preconditioner:
   - Did it build successfully?
   - Is it appropriate for this matrix type?

Step 2: Convergence Analysis

1. Plot log(||r_k||) vs k
2. Compute average convergence rate
3. Identify pattern (decay, stagnation, oscillation)
4. Compare to expected behavior for solver type

Step 3: Condition Number Estimation

1. Run CG/GMRES without preconditioner
2. Estimate κ from convergence rate:
   ρ ≈ (√κ - 1)/(√κ + 1) for CG
3. If κ > 10⁶, scaling/preconditioning essential

Step 4: Preconditioner Quality Check

1. Compare iterations with/without preconditioner
2. If minimal improvement:
   - Preconditioner too weak
   - Wrong type for this matrix
   - Implementation bug

3. Target: 10-50× fewer iterations with preconditioner

Convergence Criteria

Relative Residual

Stop when: ||r_k|| / ||b|| < tol
Typical tol: 1e-6 to 1e-10

Caution: Can be fooled by ill-conditioned systems.

Absolute Residual

Stop when: ||r_k|| < tol
Use when: ||b|| is very small or known scale

True Error (if available)

Stop when: ||x - x_k|| / ||x|| < tol
Best criterion but rarely computable

Practical Multi-Criteria

def converged(r_k, r_0, b, x_k, tol_rel=1e-6, tol_abs=1e-10):
    rel_resid = np.linalg.norm(r_k) / np.linalg.norm(b)
    abs_resid = np.linalg.norm(r_k)
    reduction = np.linalg.norm(r_k) / np.linalg.norm(r_0)

    # Any criterion met
    return (rel_resid < tol_rel or
            abs_resid < tol_abs or
            reduction < tol_rel)

Improving Convergence

General Strategies

Preconditioner Strengthening Ladder

Level 0: No preconditioner
  ↓ slow →
Level 1: Jacobi (diagonal)
  ↓ slow →
Level 2: ILU(0) or IC(0)
  ↓ slow →
Level 3: ILU(1) or IC(1)
  ↓ slow →
Level 4: ILUT with moderate fill
  ↓ slow →
Level 5: AMG or high-fill ILU
  ↓ slow →
Level 6: Direct solver

When to Switch Solvers

Implementation: Convergence Logger

class ConvergenceMonitor:
    """Track and analyze iterative solver convergence."""

    def __init__(self, b_norm):
        self.residuals = []
        self.b_norm = b_norm

    def log(self, r_norm):
        self.residuals.append(r_norm)

    def relative_residual(self):
        return [r / self.b_norm for r in self.residuals]

    def convergence_rate(self, window=5):
        """Average convergence rate over last window iterations."""
        if len(self.residuals) < window + 1:
            return None
        recent = self.residuals[-window-1:]
        rates = [recent[i+1]/recent[i] for i in range(window)]
        return np.mean(rates)

    def diagnose(self):
        if len(self.residuals) < 2:
            return "Insufficient data"

        rate = self.convergence_rate()
        if rate is None:
            return "Need more iterations"

        rel_final = self.residuals[-1] / self.b_norm

        if rel_final < 1e-10:
            return "Converged well"
        elif rate > 0.99:
            return "Stagnated - strengthen preconditioner"
        elif rate > 0.9:
            return "Slow - consider better preconditioner"
        elif rate < 0:
            return "Oscillating - try different solver"
        elif rate > 1:
            return "Diverging - check matrix properties"
        else:
            return "Good progress - continue"

    def plot(self):
        """Plot convergence history."""
        import matplotlib.pyplot as plt
        plt.semilogy(self.residuals, 'b-o')
        plt.xlabel('Iteration')
        plt.ylabel('Residual norm')
        plt.title(f'Convergence: final rate = {self.convergence_rate():.3f}')
        plt.grid(True)
        plt.show()

Quick Troubleshooting Guide

Preconditioner Catalog

Comprehensive reference for preconditioners in iterative linear solvers.

Preconditioner Fundamentals

Purpose

Transform Ax = b into M⁻¹Ax = M⁻¹b where M⁻¹A has better conditioning.

Goals:

• Cluster eigenvalues away from zero
• Reduce condition number κ(M⁻¹A)
• Make M⁻¹ cheap to apply

Application Modes

Key Trade-offs

Incomplete Factorization Family

Incomplete Cholesky (IC)

For SPD matrices: A ≈ LLᵀ where L is sparse.

IC(0) - Zero fill-in:

• Same sparsity pattern as lower triangle of A
• Cheap, often effective
• May fail for indefinite or poorly scaled

IC(k) - Level-k fill:

• Allow fill-in up to k levels from original pattern
• More robust, higher cost
• k = 1 or 2 often sufficient

Modified IC (MIC):

• Add dropped entries to diagonal
• Better for M-matrices (e.g., Laplacian)
• Preserves row sums

Incomplete LU (ILU)

For general matrices: A ≈ LU where L, U are sparse.

ILU(0) - Zero fill-in:

Pattern(L + U) = Pattern(A)
Cheap, first try for nonsymmetric

ILU(k) - Level-k fill:

Allow fill paths up to k edges
k = 1: moderate fill
k = 2: substantial fill

ILUT - Threshold-based:

Parameters: τ (drop tolerance), p (max fill per row)
Drop if |entry| < τ × ||row||
Keep at most p entries per row

Typical values:

τ = 1e-4 to 1e-2
p = 10 to 50 (or 2× to 5× original nnz/row)

Choosing IC vs ILU Parameters

Algebraic Multigrid (AMG)

When to Use

AMG Components

Coarsening:

• Classical (Ruge-Stüben): Strength-based C/F splitting
• Aggregation: Group nodes into aggregates
• Smoothed aggregation: SA-AMG, good for elasticity

Interpolation:

• Direct: Use strong connections
• Standard: Include weak connections
• Extended+i: For harder problems

Smoothing:

• Jacobi: Simple, parallelizable
• Gauss-Seidel: Better smoothing, less parallel
• Polynomial: Good for GPU

AMG Tuning

AMG for Nonsymmetric

Some AMG can handle mildly nonsymmetric:

• Use symmetric part for coarsening
• May need more smoothing
• Verify convergence experimentally

Specialized Preconditioners

Jacobi and Block Jacobi

Point Jacobi: M = diag(A)

Simple, parallel, weak
Use as smoother in multigrid

Block Jacobi: M = block_diag(A)

Blocks from natural structure (elements, nodes)
Stronger than point Jacobi
Embarrassingly parallel

Gauss-Seidel

Forward/Backward GS:

M = L + D (forward) or D + U (backward)
Stronger than Jacobi
Sequential, hard to parallelize

Symmetric GS (SSOR):

Apply forward then backward
Good smoother for symmetric problems
Parameter ω (relaxation): typically 1.0

SSOR (Symmetric SOR)

For SPD systems:

M = (D/ω + L) × (D/ω)⁻¹ × (D/ω + U)
ω: overrelaxation parameter (0 < ω < 2)
ω = 1: Symmetric Gauss-Seidel
ω optimal ≈ 2/(1 + sin(πh)) for model problem

Polynomial Preconditioners

Approximate M⁻¹ ≈ p(A) for some polynomial p.

Neumann series:

M⁻¹ ≈ I + (I - A) + (I - A)² + ...
Requires ρ(I - A) < 1

Chebyshev:

Optimal polynomial for given eigenvalue bounds
Requires [λ_min, λ_max] estimates
Good for GPU (only matrix-vector products)

Block Preconditioners

For Saddle-Point Systems

K = [A   B ]    x = [u]    b = [f]
    [Bᵀ  -C]        [p]        [g]

Block Diagonal:

P = [Â   0 ]
    [0   Ŝ ]

 ≈ A (e.g., AMG for A)
Ŝ ≈ S = C + BᵀA⁻¹B (Schur complement)

Block Triangular:

P = [Â   0 ]    or    P = [Â   B]
    [Bᵀ  -Ŝ]              [0  -Ŝ]

Schur Complement Approximations

Field-Split Preconditioners

For multi-physics with fields u₁, u₂, ...:

Multiplicative: Solve u₁, then u₂ using updated u₁
Additive: Solve each independently, sum corrections

Preconditioner Selection Guide

By Matrix Type

By Problem Size

By Available Resources

Troubleshooting

Preconditioner Fails to Build

Poor Convergence Despite Preconditioner

Parameter Tuning Strategy

1. Start simple: ILU(0) or IC(0) 2. Monitor: iteration count, residual curve 3. Adjust: If slow, strengthen preconditioner 4. Balance: Setup time vs iteration time 5. Validate: Check solution accuracy

Implementation Notes

Setup vs Apply Cost

Reusing Preconditioners

When matrix changes slightly:

Same pattern: May reuse structure, update values
Small changes: Lag preconditioner (update every k solves)
Large changes: Rebuild preconditioner

Quality Metrics

Scaling and Equilibration

Comprehensive guide for matrix scaling to improve conditioning and solver performance.

When Scaling is Needed

Indicators of Poor Scaling

Matrix Value Ranges

Scaling Methods

Row Scaling

Multiply each row by a scalar: D_r × A

For each row i:
    scale_i = 1 / max_j |A_ij|  (max scaling)
    or
    scale_i = 1 / ||row_i||₂    (norm scaling)
    or
    scale_i = 1 / ||row_i||₁    (sum scaling)

Effect: All row maxima (or norms) become 1.

Preserves: Column space structure.

Column Scaling

Multiply each column by a scalar: A × D_c

For each column j:
    scale_j = 1 / max_i |A_ij|

Effect: All column maxima become 1.

Preserves: Row space structure.

Row and Column Scaling

Apply both: D_r × A × D_c

Balanced approach:

Scaled system: (D_r × A × D_c) × (D_c⁻¹ × x) = D_r × b
Solve for: y = D_c⁻¹ × x
Recover: x = D_c × y

Symmetric Scaling

For symmetric matrices, use same scaling for rows and columns: D × A × D

Preserves symmetry: (DAD)ᵀ = DAD
Essential for CG, MINRES

Computing symmetric scaling:

d_i = 1 / sqrt(|A_ii|)  (diagonal scaling)
or
d_i = 1 / sqrt(max_j |A_ij|)  (max scaling)

Equilibration Algorithms

Ruiz Equilibration

Iterative algorithm to make row and column norms approximately equal:

Repeat until converged:
    For each row i: r_i = ||row_i||_∞
    For each col j: c_j = ||col_j||_∞
    D_r = diag(1/sqrt(r_i))
    D_c = diag(1/sqrt(c_j))
    A = D_r × A × D_c

Properties:

• Converges to doubly stochastic-like scaling
• Usually 5-10 iterations sufficient
• Works for nonsymmetric matrices

Sinkhorn-Knopp (Doubly Stochastic)

For nonnegative matrices, scale to doubly stochastic:

All row sums = 1
All column sums = 1

Algorithm:

Repeat:
    Normalize rows to sum to 1
    Normalize columns to sum to 1
Until converged

Geometric Mean Scaling

Scale by geometric mean of max and min absolute values:

For row i:
    r_max = max_j |A_ij|
    r_min = min_{j: A_ij≠0} |A_ij|
    scale_i = 1 / sqrt(r_max × r_min)

Good for: Matrices with entries spanning many orders of magnitude.

Special Considerations

Preserving Symmetry

CRITICAL: For symmetric matrices, use symmetric scaling only!

Wrong: D_r × A × D_c with D_r ≠ D_c (destroys symmetry)
Right: D × A × D with same D for rows and columns

Check after scaling:

||A_scaled - A_scaled^T||_F / ||A_scaled||_F < ε

Preserving Positive Definiteness

Symmetric scaling preserves definiteness:

If A is SPD and D is nonsingular diagonal:
Then D × A × D is also SPD

Check after scaling (if needed):

Try Cholesky factorization
Or check smallest eigenvalue > 0

Scaling the RHS

When scaling Ax = b → (D_r A D_c) y = D_r b:

Solve: (D_r A D_c) y = D_r b
Recover: x = D_c y
Residual: r = b - Ax = b - A(D_c y)

Important: Scale b with the same D_r used for rows!

Practical Scaling Strategies

Safe Default Strategy

def scale_matrix(A, symmetric=False):
    # Row scaling
    row_norms = np.max(np.abs(A), axis=1)
    D_r = 1.0 / np.maximum(row_norms, 1e-15)

    if symmetric:
        # Use same scaling for rows and columns
        D_c = D_r
    else:
        # Scale by rows first, then columns
        A_scaled = np.diag(D_r) @ A
        col_norms = np.max(np.abs(A_scaled), axis=0)
        D_c = 1.0 / np.maximum(col_norms, 1e-15)

    return D_r, D_c

When to Apply Scaling

Scaling Order

1. Analyze original matrix (symmetry, definiteness) 2. Choose appropriate scaling method 3. Scale matrix and RHS 4. Build preconditioner on scaled system 5. Solve scaled system 6. Unscale solution

Diagnosing Scaling Issues

Before Scaling

Check: ratio of max to min nonzero absolute values
       ||A||_∞ / min_{ij: A_ij≠0} |A_ij|

If ratio > 10⁶: scaling strongly recommended

After Scaling

Verify:
- Row norms approximately equal
- Column norms approximately equal
- No overflow/underflow
- Symmetry preserved (if was symmetric)

Monitoring Convergence

Common Mistakes

Mistake 1: Destroying Symmetry

Problem: Used different row and column scaling on symmetric matrix
Result: CG fails (requires symmetric)
Fix: Use symmetric scaling D × A × D

Mistake 2: Forgetting to Unscale

Problem: Returned y instead of x = D_c × y
Result: Wrong solution
Fix: Always unscale: x = D_c × y

Mistake 3: Not Scaling RHS

Problem: Scaled A but not b
Result: Solving wrong system
Fix: b_scaled = D_r × b

Mistake 4: Scaling Zero Rows

Problem: Row of zeros gets 1/0 = inf scaling
Result: NaN in solution
Fix: Handle zero rows specially (or they indicate singular matrix)

Advanced Topics

Iterative Refinement with Scaling

For very ill-conditioned systems:

1. Scale: Ã = D_r A D_c, b̃ = D_r b
2. Solve: Ã ỹ = b̃ (in lower precision if available)
3. Compute residual: r = b - A(D_c ỹ) (in high precision)
4. Solve: Ã δỹ = D_r r
5. Update: ỹ = ỹ + δỹ
6. Repeat until converged
7. Unscale: x = D_c ỹ

Scaling for Eigenvalue Problems

When solving Ax = λx after scaling D_r A D_c:

Eigenvalues: Same (scaling is similarity transform)
Eigenvectors: v_original = D_c × v_scaled

Block Scaling

For block-structured matrices:

Scale each block independently, or
Scale entire rows/columns of blocks

Useful for: Multiphysics with different variable scales.

Implementation Example

def equilibrate_matrix(A, tol=1e-6, max_iter=20):
    """Ruiz equilibration for general matrices."""
    n, m = A.shape
    D_r = np.ones(n)
    D_c = np.ones(m)

    for _ in range(max_iter):
        # Row scaling
        row_norms = np.max(np.abs(A), axis=1)
        row_norms = np.maximum(row_norms, 1e-15)
        d_r = 1.0 / np.sqrt(row_norms)
        A = np.diag(d_r) @ A
        D_r *= d_r

        # Column scaling
        col_norms = np.max(np.abs(A), axis=0)
        col_norms = np.maximum(col_norms, 1e-15)
        d_c = 1.0 / np.sqrt(col_norms)
        A = A @ np.diag(d_c)
        D_c *= d_c

        # Check convergence
        if np.max(np.abs(row_norms - 1)) < tol:
            if np.max(np.abs(col_norms - 1)) < tol:
                break

    return A, D_r, D_c

def solve_scaled(A, b, solver, symmetric=False):
    """Solve with scaling and unscaling."""
    A_scaled, D_r, D_c = equilibrate_matrix(A)
    b_scaled = D_r * b

    if symmetric:
        # Ensure symmetric scaling
        D = np.sqrt(D_r * D_c)
        A_scaled = np.diag(D) @ A @ np.diag(D)
        b_scaled = D * b
        y = solver(A_scaled, b_scaled)
        x = D * y
    else:
        y = solver(A_scaled, b_scaled)
        x = D_c * y

    return x

Quick Reference

Solver Decision Tree

Comprehensive decision guide for selecting linear solvers for Ax = b.

Matrix Classification

Key Properties to Determine

Quick Classification

Check symmetry: ||A - Aᵀ||_F / ||A||_F < ε
Check SPD: try Cholesky, if succeeds → SPD
Check sparsity: nnz / n² gives density
Check conditioning: estimate κ using power iteration or Lanczos

Primary Decision Tree

START: Need to solve Ax = b
│
├─ Is n < 5000 and matrix dense?
│   └── YES → Use DIRECT solver
│       ├── Symmetric PD → Cholesky (LLᵀ)
│       ├── Symmetric indefinite → LDLᵀ (Bunch-Kaufman)
│       ├── Nonsymmetric → LU with pivoting
│       └── Least squares → QR or SVD
│
└── NO → Use ITERATIVE solver
    │
    ├─ Is matrix symmetric?
    │   │
    │   ├── YES → Is it positive definite?
    │   │   │
    │   │   ├── YES (SPD) → CG (Conjugate Gradient)
    │   │   │   ├── Elliptic PDE → AMG preconditioner
    │   │   │   ├── Banded structure → IC(0) or IC(k)
    │   │   │   └── General sparse → AMG or IC
    │   │   │
    │   │   └── NO/Unknown (symmetric indefinite)
    │   │       ├── Eigenvalues both signs → MINRES
    │   │       ├── Near-singular → SYMMLQ
    │   │       └── Preconditioner: ILU or block diagonal
    │   │
    │   └── NO (nonsymmetric) → Is it nearly symmetric?
    │       │
    │       ├── YES (A ≈ Aᵀ) → BiCGSTAB
    │       │   ├── Smooth convergence needed → BiCGSTAB(ℓ)
    │       │   └── Preconditioner: ILUT or AMG
    │       │
    │       └── NO (strongly nonsymmetric) → GMRES
    │           ├── Memory limited → GMRES(m) restarted
    │           ├── Very nonsymmetric → Full GMRES if affordable
    │           └── Preconditioner: ILU(k), ILUT, or AMG

Direct Solver Selection

When to Use Direct

Direct Solver Types

Memory Considerations

Iterative Solver Details

Conjugate Gradient (CG)

Requirements: A must be SPD

Convergence:

||e_k||_A ≤ 2 × ((√κ - 1)/(√κ + 1))^k × ||e_0||_A

Breakdown: CG is breakdown-free for SPD matrices.

GMRES

Requirements: None (works for any nonsingular matrix)

Properties:

• Minimizes residual over Krylov subspace
• Optimal for nonsymmetric systems
• Memory: O(m × n) for m iterations

Restart Strategy:

GMRES(m): Restart after m iterations
├── m = 20-30: Typical for memory-limited
├── m = 50-100: Better convergence
└── m = n: Full GMRES (no restart)

When GMRES stalls:

• Increase restart parameter m
• Improve preconditioner
• Check for near-singularity

BiCGSTAB

Requirements: Works for nonsymmetric, best if A ≈ Aᵀ

Properties:

• Lower memory than GMRES
• May have irregular convergence
• Two matrix-vector products per iteration

Variants:

MINRES

Requirements: A symmetric (not necessarily PD)

Properties:

• Works for indefinite symmetric systems
• Minimizes residual (like GMRES for symmetric)
• Short recurrences (memory efficient)

Use for:

• Saddle-point systems
• Systems from variational problems
• Eigenvalue problems (shift-invert)

Special Matrix Structures

Saddle-Point Systems

[A   B ] [x]   [f]
[Bᵀ  -C] [y] = [g]

Solvers:

• Schur complement method
• Uzawa iteration
• Block preconditioned GMRES/MINRES

Preconditioners:

Block diagonal: [A⁻¹   0  ]
                [0    S⁻¹]

Block triangular: [A⁻¹  0 ]
                  [Bᵀ   S⁻¹]

Where S = C + BᵀA⁻¹B is the Schur complement.

Banded Systems

For bandwidth b << n:

Block Systems

When A has natural block structure:

Block Jacobi: Solve diagonal blocks independently
Block Gauss-Seidel: Sequential block updates
Block ILU: Factor with block structure preserved

Solver Selection by Application

Elliptic PDEs (Laplacian-like)

Parabolic PDEs (Heat equation)

Hyperbolic PDEs (Advection-dominated)

Navier-Stokes

Phase-Field

Failure Modes and Remedies

Common Problems

When All Else Fails

1. Check the matrix: Is it actually nonsingular? 2. Scale rows/columns: Equilibrate to unit row norms 3. Try direct solver: Even for large sparse, may work 4. Regularize: Add small diagonal for near-singular 5. Reformulate problem: Different discretization

Quick Reference Table