Problems solving linear system with multiple right-hand side

I am trying to use cusolverSpZcsrlsvqr with multiple right hand side (~20 000) but it is very inefficient. Is there another culib solution of that kind of problems?

void GpuMatrix::solve(const d_sparse_matrix &A, const d_dense_matrix &d_b, d_dense_matrix &d_x) {
    auto m = A.m_;
    auto n = d_b.n_;
    auto nnz = A.nnz_;
    int singularity;
    int reorder = 0;
    cusparseMatDescr_t descr;
    cusparseSetMatIndexBase(descr, CUSPARSE_INDEX_BASE_ZERO);
    cusparseSetMatType(descr, CUSPARSE_MATRIX_TYPE_GENERAL);
    std::cout << "rhs= " << n << std::endl;
    for (int i = 0; i < n; ++i) {
        int offset = i*m;
        CHECK_CUSOLVER(cusolverSpZcsrlsvqr(sv_handle_, m, nnz, descr, reinterpret_cast<const cuDoubleComplex *>(A.vm),
                                 ,,reinterpret_cast<const cuDoubleComplex *>(d_b.vm+offset),
                                           tol_,  reorder,reinterpret_cast<cuDoubleComplex *>(d_x.vm+offset), &singularity));