BLIT — Block Iterative Linear Solver Toolbox

% 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;
}