sparse-matrix-krylov

Answer

The implementation provides a complete parallel GMRES solver for large sparse linear systems with three core components:

1. CSR Sparse Matrix Format

The CSRMatrix class stores matrices in Compressed Sparse Row format with values, col_indices, and row_ptr arrays. It provides efficient matvec() (sparse matrix-vector product) and a threaded matvec_parallel() for parallel execution.

2. ILU(0) Preconditioner

The ILU0Preconditioner computes an incomplete LU factorization with zero fill-in — the sparsity pattern of L and U matches A's pattern exactly. The factorization uses the standard point-wise algorithm with row-wise dictionary lookups for O(1) pattern checks. Forward and backward substitution implement the solve.

3. Restarted GMRES with Right Preconditioning

The GMRES class implements restarted GMRES(m) using: - Right preconditioning: Solve A M^{-1} u = b where x = M^{-1} u, with the ILU(0) factor as M - Arnoldi iteration: Modified Gram-Schmidt for stable orthogonalization, building an orthonormal basis of the Krylov subspace K_m(A M^{-1}, r_0) - Givens rotations: For stable least-squares solution of the Hessenberg system without explicitly forming Q - Restart: After m Krylov vectors, the cycle repeats with the new residual

# Core GMRES cycle with right preconditioning
def _gmres_cycle(self, x, b, m):
    r = b - self.A.matvec(x)
    beta = np.linalg.norm(r)
    V = np.zeros((n, m + 1)); V[:, 0] = r / beta

    for k in range(m):
        # Right-preconditioned Arnoldi: w = A * (M^{-1} * v_k)
        z = self.M.solve(V[:, k]) if self.M else V[:, k].copy()
        w = self.A.matvec(z)

        # Modified Gram-Schmidt orthogonalization
        for j in range(k + 1):
            H[j, k] = np.dot(V[:, j], w)
            w -= H[j, k] * V[:, j]

        H[k + 1, k] = np.linalg.norm(w)
        V[:, k + 1] = w / H[k + 1, k]

        # Apply Givens rotations to track residual norm cheaply
        # ... (see full implementation)

    # Solve upper triangular system: R*y = g
    # Update: x += M^{-1} * (V * y)
    temp = V[:, :k_used] @ y
    delta = self.M.solve(temp) if self.M else temp
    return x + delta, k_used, converged

4. Utility Functions

• laplacian_1d(n) / laplacian_2d(n) — generate standard test matrices (tridiagonal and 5-point stencil)
• random_sparse(n, density) — random diagonally dominant SPD matrices
• gmres(...) — convenience wrapper