Solver library for matrix-free linear equation system

Question

I will have to solve a large linear system. I'm now looking for a solver that works "matrix-free" (So that I just have to specify a matrix-vector product, but not the matrix). As far as I understand (I'm not an expert.) LAPACK does not provide this option.

Can you recommend any libraries? Or is my approach to implement this matrix-free not reasonable at all?

Remark: I plan to use C++ and the matrix comes up from a FEM for a system of PDEs, but is not symmetric. The matrix is sparse, but as I plan to provide the matrix-vector product by myself, and the resulting vector will be dense, sparsity should not play a role here - as far as I understand.

Regards, Michael

Answer 1

Using a matrix-free approach for a sufficiently large linear system sounds reasonable to me.

You probably want to take a look at PETSc or Trilinos; those are the two main libraries that serve as abstraction layers for a number of different numerical methods (preconditioners, sparse linear solvers, scalable nonlinear equation solvers, time steppers, optimization). They also build in abstractions that make it easier to incorporate parallelism into an application than using bare MPI calls.

You could also look at higher levels of abstraction, if that is of interest; there are a number of finite element libraries that build on PETSc or Trilinos.