Showing revision 2.8

BLIT — Block Iterative Linear Solvers

High-performance sparse solvers for symmetric FEM matrices

Real Symmetric Complex Symmetric (A = Aᵀ) Sparse & Dense
Python
MATLAB
Fortran
C++
Get Started GitHub

Why Block QMR?

🧮

Block Solver

Multiple RHS simultaneously, sharing Krylov subspace for 4× fewer iterations

∑

Complex Symmetry

Quasi inner product ⟨x,y⟩ = Σxₖyₖ — natural for A = Aᵀ systems

🛡️

Preconditioned

Built-in ILU, Jacobi, and split preconditioners

⚡

Up to 65× Faster

Outperforms direct solvers for large systems

💾

Memory Efficient

Short three-term recurrences, not full Krylov basis

📈

Smooth Convergence

Monotonic quasi-minimal residual reduction

Applications

🔬 Diffuse Optical Tomography
📡 Frequency-Domain EM
🔊 Acoustic Scattering
🩺 Microwave Imaging
🌊 Helmholtz Equations
⚙️ Structural Dynamics

Simple API

# pip install blocksolver
from blocksolver import blqmr
from scipy.sparse import diags

# Complex symmetric FEM matrix
n = 5000
A = diags([-1, 4+0.1j, -1], [-1, 0, 1], shape=(n, n), dtype=complex)

# Block solve: 16 RHS at once
B = np.random.randn(n, 16) + 0j
result = blqmr(A, B, tol=1e-10, precond_type='diag')

print(f"Converged: {result.converged}")
% Complex symmetric FEM matrix
n = 5000;
K = spdiags([-ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
A = K + 0.1i * speye(n);

% Block solve: 16 RHS at once
B = rand(n, 16) + 1i*rand(n, 16);
opt.precond = 'diag';
[X, flag] = blqmr(A, B, 1e-10, 1000, [], [], [], opt);
program solve_fem
    use blit_blqmr_real
    implicit none
    
    type(BLQMRSolver) :: qmr
    integer :: n, nnz, nrhs
    integer, allocatable :: Ap(:), Ai(:)
    real(8), allocatable :: Ax(:), B(:,:), X(:,:)
    
    n = 5000; nnz = 14998; nrhs = 16
    ! ... allocate and fill CSC arrays ...
    
    call BLQMRCreate(qmr, n)
    qmr%maxit = 1000
    qmr%qtol = 1.0d-10
    qmr%pcond_type = 3  ! Jacobi
    
    call BLQMRPrep(qmr, Ap, Ai, Ax, nnz)
    call BLQMRSolve(qmr, Ap, Ai, Ax, nnz, X, B, nrhs)
    call BLQMRDestroy(qmr)
end program
#include "blit_solvers.h"

int main() {
    const int n = 5000, nnz = 14998, nrhs = 16;
    std::vector<int> Ap(n+1), Ai(nnz);
    std::vector<double> Ax(nnz), b(n*nrhs), x(n*nrhs);
    
    BlitBLQMR<double> solver(n, nrhs, 1000, 1e-3);
    solver.Prepare(Ap.data(), Ai.data(), Ax.data(), nnz);
    solver.Solve(x.data(), b.data(), nrhs);
    return 0;
}

Performance

65×
Max Speedup vs Direct
16-32
Optimal Block Size
<0.25
Super-linear Iter Ratio
3K+
Crossover Nodes

📊 Speedup vs Direct Solver (Complex, 4 RHS)

Finding: BLQMR achieves 14× speedup at 125K nodes. For single RHS, speedup reaches 65×!

⏱️ Solve Time Comparison

Scaling: Direct solver time grows as O(n²), while BLQMR scales nearly linearly.

🔢 Block Size Optimization (64 RHS, 20³)

Optimal: Block size 16-32 balances iteration reduction with per-iteration cost.

📉 Iteration Efficiency

Super-linear: Actual ratio beats theoretical 1/m — Krylov info sharing works!

🎯 Crossover Analysis: BLQMR vs Direct Solver

Speedup = Direct_time / BLQMR_time. Green (>1.0) = BLQMR wins.

Grid1 RHS2 RHS4 RHS8 RHS16 RHS32 RHS64 RHS128 RHS
10³ (1K)0.650.450.440.290.190.160.150.13
15³ (3K)4.41.30.960.610.340.260.220.20
20³ (8K)6.53.32.21.30.690.530.420.35
25³ (16K)10.34.42.71.50.940.730.570.51
30³ (27K)16.39.15.43.21.91.31.00.77
BLQMR Wins
Break-even
Direct Wins
Key Insight: For n ≥ 3,375 nodes, BLQMR wins with few RHS. Larger systems benefit even with many RHS. The crossover shifts right as problem size decreases.

Get Started

1

Python

pip install blocksolver

2

Fortran Backend (Optional)

Linux: apt install gfortran libsuitesparse-dev

3

MATLAB/Octave

addpath('/path/to/blit/matlab')

References

Microwave image reconstruction from 3-D fields coupled to 2-D parameter estimation
Fang Q, Meaney PM, Geimer SD, et al.
IEEE Transactions on Medical Imaging, 23(4):475-484, 2004
A block QMR method for computing multiple simultaneous solutions
Boyse WE, Seidl AA
SIAM Journal on Scientific Computing, 17(1):263-274, 1996
Powered by Habitat