There are several libraries to (iteratively) solve large sparse linear equation systems in parallel on a number of CPUs. Our parallel cluster also has attached powerful GPUs, but so far, I did not find any solver that can use CPUs and GPUs in parallel to solve a sparse linear equation system.

MAIN QUESTION: Are there actually no solvers which use CPUs and GPUs or could you recommend me one?

RELATED SIDE QUESTION: Many of the solver libraries can work with an user defined matrix-vector-product which is the basic operation of any iterative solver. So I could use the GPUs to help the CPUs computing the matrix-vector-product, BUT the problem there is: When not specifying the matrix explicitly, but just giving the matrix-vector-product, the (very important) preconditioner routines can not be used. Any ideas how to get over this problem?

  • 1
    $\begingroup$ Some preconditioners (such as incomplete LU factorization) use the actual matrix entries and cannot be used when you only have functions for computing matrix vector products. If you're in this situation, you'd have to come up with a problem specific preconditioner and implement it as a function. If you're willing to discuss your problem in more detail, then perhaps we could help you to come up with a problem specific preconditioner. Without further information, the question really can't be answered. $\endgroup$ May 14, 2014 at 17:35

1 Answer 1


Linear system solvers are generally limited by memory access, so in parallel systems communication becomes the bottleneck. Good scaling can be usually achieved only for very large systems (millions of unknowns), for smaller systems you might be better off with a single-node or single-GPU solver. The libraries with GPU support that I know are listed below, for the "usual" CPU solvers see Libraries for solving sparse linear systems.

As for the custom matrix-vector product, one of the main issues with these libraries is that each uses their own representation of sparse matrices. Most of the time they are based on a standard format like CSR, CSC or Ellpack, so an efficient conversion might be possible, but there is no chance for bitwise compatibility. The problem is that it is not sufficient to optimize the operations, which could be generalized by user-defined operators, but also the memory access pattern for each operation has to be considered. If you can think of a solution that is both general and efficient, please let me know (you should probably write a paper first).

CUDA-only libraries:

  • cuSOLVER provides several direct methods for both dense and sparse systems. It is closed-source, bundled with the CUDA toolkit.
  • cuSPARSE provides incomplete factorizations, which can be used as preconditioners for iterative methods. It is closed-source, bundled with the CUDA toolkit.
  • CUSP provides iterative methods and multiple preconditioners, including a smoothed-aggregation algebraic multigrid. It is free and open-source.

Hybrid or CUDA-accelerated libraries:

  • SuiteSparse has a CUDA-accelerated Cholesky factorization. It is free and open-source.
  • SuperLU_DIST has a CUDA-accelerated LU factorization. It is free and open-source.
  • ViennaCL is a library using CUDA, OpenCL and OpenMP for parallelization, but I don't know if one computation can use multiple backends in parallel. It provides multiple iterative solvers and preconditioners. It is free and open-source.
  • AmgX provides algebraic multigrid and preconditioned iterative methods. It uses CUDA, MPI and OpenMP for parallelization. It is available with a commercial and a free license, the latter is limited to registered developers and non-commercial use.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.