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.