def _cusp_solver(M, parameters):
    cache_key = lambda t, p: (t,
                              p['ksp_type'],
                              p['pc_type'],
                              p['ksp_rtol'],
                              p['ksp_atol'],
                              p['ksp_max_it'],
                              p['ksp_gmres_restart'],
                              p['ksp_monitor'])
    module = _cusp_cache.get(cache_key(M.ctype, parameters))
    if module:
        return module

    import codepy.toolchain
    from cgen import FunctionBody, FunctionDeclaration
    from cgen import Block, Statement, Include, Value
    from codepy.bpl import BoostPythonModule
    from codepy.cuda import CudaModule
    gcc_toolchain = codepy.toolchain.guess_toolchain()
    nvcc_toolchain = codepy.toolchain.guess_nvcc_toolchain()
    if 'CUSP_HOME' in os.environ:
        nvcc_toolchain.add_library('cusp', [os.environ['CUSP_HOME']], [], [])
    host_mod = BoostPythonModule()
    nvcc_mod = CudaModule(host_mod)
    nvcc_includes = ['thrust/device_vector.h',
                     'thrust/fill.h',
                     'cusp/csr_matrix.h',
                     'cusp/krylov/cg.h',
                     'cusp/krylov/bicgstab.h',
                     'cusp/krylov/gmres.h',
                     'cusp/precond/diagonal.h',
                     'cusp/precond/smoothed_aggregation.h',
                     'cusp/precond/ainv.h',
                     'string']
    nvcc_mod.add_to_preamble([Include(s) for s in nvcc_includes])
    nvcc_mod.add_to_preamble([Statement('using namespace std')])

    # We're translating PETSc preconditioner types to CUSP
    diag = Statement('cusp::precond::diagonal< ValueType, cusp::device_memory >M(A)')
    ainv = Statement(
        'cusp::precond::scaled_bridson_ainv< ValueType, cusp::device_memory >M(A)')
    amg = Statement(
        'cusp::precond::smoothed_aggregation< IndexType, ValueType, cusp::device_memory >M(A)')
    none = Statement(
        'cusp::identity_operator< ValueType, cusp::device_memory >M(nrows, ncols)')
    preconditioners = {
        'diagonal': diag,
        'jacobi': diag,
        'ainv': ainv,
        'ainvcusp': ainv,
        'amg': amg,
        'hypre': amg,
        'none': none,
        None: none
    }
    try:
        precond_call = preconditioners[parameters['pc_type']]
    except __HOLE__:
        raise RuntimeError("Cusp does not support preconditioner type %s" %
                           parameters['pc_type'])
    solvers = {
        'cg': Statement('cusp::krylov::cg(A, x, b, monitor, M)'),
        'bicgstab': Statement('cusp::krylov::bicgstab(A, x, b, monitor, M)'),
        'gmres': Statement('cusp::krylov::gmres(A, x, b, %(ksp_gmres_restart)d, monitor, M)' % parameters)
    }
    try:
        solve_call = solvers[parameters['ksp_type']]
    except KeyError:
        raise RuntimeError("Cusp does not support solver type %s" %
                           parameters['ksp_type'])
    monitor = 'monitor(b, %(ksp_max_it)d, %(ksp_rtol)g, %(ksp_atol)g)' % parameters

    nvcc_function = FunctionBody(
        FunctionDeclaration(Value('void', '__cusp_solve'),
                            [Value('CUdeviceptr', '_rowptr'),
                             Value('CUdeviceptr', '_colidx'),
                             Value('CUdeviceptr', '_csrdata'),
                             Value('CUdeviceptr', '_b'),
                             Value('CUdeviceptr', '_x'),
                             Value('int', 'nrows'),
                             Value('int', 'ncols'),
                             Value('int', 'nnz')]),
        Block([
            Statement('typedef int IndexType'),
            Statement('typedef %s ValueType' % M.ctype),
            Statement(
                'typedef typename cusp::array1d_view< thrust::device_ptr<IndexType> > indices'),
            Statement(
                'typedef typename cusp::array1d_view< thrust::device_ptr<ValueType> > values'),
            Statement(
                'typedef cusp::csr_matrix_view< indices, indices, values, IndexType, ValueType, cusp::device_memory > matrix'),
            Statement('thrust::device_ptr< IndexType > rowptr((IndexType *)_rowptr)'),
            Statement('thrust::device_ptr< IndexType > colidx((IndexType *)_colidx)'),
            Statement('thrust::device_ptr< ValueType > csrdata((ValueType *)_csrdata)'),
            Statement('thrust::device_ptr< ValueType > d_b((ValueType *)_b)'),
            Statement('thrust::device_ptr< ValueType > d_x((ValueType *)_x)'),
            Statement('indices row_offsets(rowptr, rowptr + nrows + 1)'),
            Statement('indices column_indices(colidx, colidx + nnz)'),
            Statement('values matrix_values(csrdata, csrdata + nnz)'),
            Statement('values b(d_b, d_b + nrows)'),
            Statement('values x(d_x, d_x + ncols)'),
            Statement('thrust::fill(x.begin(), x.end(), (ValueType)0)'),
            Statement(
                'matrix A(nrows, ncols, nnz, row_offsets, column_indices, matrix_values)'),
            Statement('cusp::%s_monitor< ValueType > %s' %
                      ('verbose' if parameters['ksp_monitor'] else 'default',
                       monitor)),
            precond_call,
            solve_call
        ]))

    host_mod.add_to_preamble([Include('boost/python/extract.hpp'), Include('string')])
    host_mod.add_to_preamble([Statement('using namespace boost::python')])
    host_mod.add_to_preamble([Statement('using namespace std')])

    nvcc_mod.add_function(nvcc_function)

    host_mod.add_function(
        FunctionBody(
            FunctionDeclaration(Value('void', 'solve'),
                                [Value('object', '_rowptr'),
                                 Value('object', '_colidx'),
                                 Value('object', '_csrdata'),
                                 Value('object', '_b'),
                                 Value('object', '_x'),
                                 Value('object', '_nrows'),
                                 Value('object', '_ncols'),
                                 Value('object', '_nnz')]),
            Block([
                Statement(
                    'CUdeviceptr rowptr = extract<CUdeviceptr>(_rowptr.attr("gpudata"))'),
                Statement(
                    'CUdeviceptr colidx = extract<CUdeviceptr>(_colidx.attr("gpudata"))'),
                Statement(
                    'CUdeviceptr csrdata = extract<CUdeviceptr>(_csrdata.attr("gpudata"))'),
                Statement('CUdeviceptr b = extract<CUdeviceptr>(_b.attr("gpudata"))'),
                Statement('CUdeviceptr x = extract<CUdeviceptr>(_x.attr("gpudata"))'),
                Statement('int nrows = extract<int>(_nrows)'),
                Statement('int ncols = extract<int>(_ncols)'),
                Statement('int nnz = extract<int>(_nnz)'),
                Statement('__cusp_solve(rowptr, colidx, csrdata, b, x, nrows, ncols, nnz)')
            ])))

    nvcc_toolchain.cflags.append('-arch')
    nvcc_toolchain.cflags.append('sm_20')
    nvcc_toolchain.cflags.append('-O3')
    module = nvcc_mod.compile(gcc_toolchain, nvcc_toolchain, debug=configuration["debug"])

    _cusp_cache[cache_key(M.ctype, parameters)] = module
    return module

# FIXME: inherit from base while device gives us the PETSc solver
dataset/ETHPy150Open OP2/PyOP2/pyop2/cuda.py/_cusp_solver