このサイトのコードをテストしていたら以下のようなエラーが吐き出された。

import pycuda.gpuarray as gpuarray
import pycuda.driver as cuda
import pycuda.autoinit
import numpy as np
import scipy.sparse as sp
import ctypes

## Wrap the cuSOLVER cusolverSpDcsrlsvqr() using ctypes
## http://docs.nvidia.com/cuda/cusolver/#cusolver-lt-t-gt-csrlsvqr

# cuSparse
_libcusparse = ctypes.cdll.LoadLibrary('libcusparse.so')
_libcusparse.cusparseCreate.restype = int
_libcusparse.cusparseCreate.argtypes = [ctypes.c_void_p]

_libcusparse.cusparseDestroy.restype = int
_libcusparse.cusparseDestroy.argtypes = [ctypes.c_void_p]

_libcusparse.cusparseCreateMatDescr.restype = int
_libcusparse.cusparseCreateMatDescr.argtypes = [ctypes.c_void_p]


# cuSOLVER
_libcusolver = ctypes.cdll.LoadLibrary('libcusolver.so')

_libcusolver.cusolverSpCreate.restype = int
_libcusolver.cusolverSpCreate.argtypes = [ctypes.c_void_p]

_libcusolver.cusolverSpDestroy.restype = int
_libcusolver.cusolverSpDestroy.argtypes = [ctypes.c_void_p]

_libcusolver.cusolverSpDcsrlsvqr.restype = int
_libcusolver.cusolverSpDcsrlsvqr.argtypes= [ctypes.c_void_p,
                                            ctypes.c_int,
                                            ctypes.c_int,
                                            ctypes.c_void_p,
                                            ctypes.c_void_p,
                                            ctypes.c_void_p,
                                            ctypes.c_void_p,
                                            ctypes.c_void_p,
                                            ctypes.c_double,
                                            ctypes.c_int,
                                            ctypes.c_void_p,
                                            ctypes.c_void_p]

def cuspsolve(A, b):
    Acsr = sp.csr_matrix(A, dtype=float)
    b = np.asarray(b, dtype=float)
    x = np.empty_like(b)

    # Copy arrays to GPU
    dcsrVal = gpuarray.to_gpu(Acsr.data)
    dcsrColInd = gpuarray.to_gpu(Acsr.indices)
    dcsrIndPtr = gpuarray.to_gpu(Acsr.indptr)
    dx = gpuarray.to_gpu(x)
    db = gpuarray.to_gpu(b)

    # Create solver parameters
    m = ctypes.c_int(Acsr.shape[0])  # Need check if A is square
    nnz = ctypes.c_int(Acsr.nnz)
    descrA = ctypes.c_void_p()
    reorder = ctypes.c_int(0)
    tol = ctypes.c_double(1e-10)
    singularity = ctypes.c_int(0)  # -1 if A not singular

    # create cusparse handle
    _cusp_handle = ctypes.c_void_p()
    status = _libcusparse.cusparseCreate(ctypes.byref(_cusp_handle))
    assert(status == 0)
    cusp_handle = _cusp_handle.value

    # create MatDescriptor
    status = _libcusparse.cusparseCreateMatDescr(ctypes.byref(descrA))
    assert(status == 0)

    #create cusolver handle
    _cuso_handle = ctypes.c_void_p()
    status = _libcusolver.cusolverSpCreate(ctypes.byref(_cuso_handle))
    assert(status == 0)
    cuso_handle = _cuso_handle.value

    # Solve
    res=_libcusolver.cusolverSpDcsrlsvqr(cuso_handle,
                                         m,
                                         nnz,
                                         descrA,
                                         int(dcsrVal.gpudata),
                                         int(dcsrIndPtr.gpudata),
                                         int(dcsrColInd.gpudata),
                                         int(db.gpudata),
                                         tol,
                                         reorder,
                                         int(dx.gpudata),
                                         ctypes.byref(singularity))
    assert(res == 0)
    if singularity.value != -1:
        raise ValueError('Singular matrix!')
    x = dx.get()  # Get result as numpy array

    # Destroy handles
    status = _libcusolver.cusolverSpDestroy(cuso_handle)
    assert(status == 0)
    status = _libcusparse.cusparseDestroy(cusp_handle)
    assert(status == 0)

    # Return result
    return x

# Test
if __name__ == '__main__':
    A = np.diag(np.arange(1, 5))
    b = np.ones(4)
    x = cuspsolve(A, b)
    np.testing.assert_almost_equal(x, np.array([1. , 0.5, 0.33333333, 0.25]))

    # Timing comparison
    from scipy.sparse import rand
    from scipy.sparse.linalg import spsolve
    from scipy.sparse import coo_matrix
    import time
    n = 10000
    i = j = np.arange(n)
    diag = np.ones(n)
    A = rand(n, n, density=0.001)
    A = A.tocsr()
    A[i, j] = diag
    b = np.ones(n)

    t0 = time.time()
    x = spsolve(A, b)
    dt1 = time.time() - t0
    print ("scipy.sparse.linalg.spsolve time: %s" %dt1)

    t0 = time.time()
    x = cuspsolve(A, b)
    dt2 = time.time() - t0
    print ("cuspsolve time: %s" %dt2)

    ratio = dt1/dt2
    if ratio > 1:
        print ("CUDA is %s times faster than CPU." %ratio)
    else:
        print ("CUDA is %s times slower than CPU." %(1./ratio))

---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
<ipython-input-7-a34dfa94d72d> in <module>()
     22 
     23 # cuSOLVER
---> 24 _libcusolver = ctypes.cdll.LoadLibrary('libcusolver.so')
     25 
     26 _libcusolver.cusolverSpCreate.restype = int

~/.pyenv/versions/3.6.5/lib/python3.6/ctypes/__init__.py in LoadLibrary(self, name)
    424 
    425     def LoadLibrary(self, name):
--> 426         return self._dlltype(name)
    427 
    428 cdll = LibraryLoader(CDLL)

~/.pyenv/versions/3.6.5/lib/python3.6/ctypes/__init__.py in __init__(self, name, mode, handle, use_errno, use_last_error)
    346 
    347         if handle is None:
--> 348             self._handle = _dlopen(self._name, mode)
    349         else:
    350             self._handle = handle

OSError: /usr/local/cuda-9.0/targets/x86_64-linux/lib/libcusolver.so: undefined symbol: GOMP_critical_end

このエラーを解決する方法がこのサイトに載っていた。

ctypes.CDLL('libgomp.so.1', mode=ctypes.RTLD_GLOBAL)
_libcusolver = ctypes.cdll.LoadLibrary('libcusolver.so')

上記のような形に書き換えてやるとエラーは消える。

import pycuda.gpuarray as gpuarray
import pycuda.driver as cuda
import pycuda.autoinit
import numpy as np
import scipy.sparse as sp
import ctypes

## Wrap the cuSOLVER cusolverSpDcsrlsvqr() using ctypes
## http://docs.nvidia.com/cuda/cusolver/#cusolver-lt-t-gt-csrlsvqr

# cuSparse
_libcusparse = ctypes.cdll.LoadLibrary('libcusparse.so')
_libcusparse.cusparseCreate.restype = int
_libcusparse.cusparseCreate.argtypes = [ctypes.c_void_p]

_libcusparse.cusparseDestroy.restype = int
_libcusparse.cusparseDestroy.argtypes = [ctypes.c_void_p]

_libcusparse.cusparseCreateMatDescr.restype = int
_libcusparse.cusparseCreateMatDescr.argtypes = [ctypes.c_void_p]

# cuSOLVER
ctypes.CDLL('libgomp.so.1', mode=ctypes.RTLD_GLOBAL)
_libcusolver = ctypes.cdll.LoadLibrary('libcusolver.so')

_libcusolver.cusolverSpCreate.restype = int
_libcusolver.cusolverSpCreate.argtypes = [ctypes.c_void_p]

_libcusolver.cusolverSpDestroy.restype = int
_libcusolver.cusolverSpDestroy.argtypes = [ctypes.c_void_p]

_libcusolver.cusolverSpDcsrlsvqr.restype = int
_libcusolver.cusolverSpDcsrlsvqr.argtypes= [ctypes.c_void_p,
                                            ctypes.c_int,
                                            ctypes.c_int,
                                            ctypes.c_void_p,
                                            ctypes.c_void_p,
                                            ctypes.c_void_p,
                                            ctypes.c_void_p,
                                            ctypes.c_void_p,
                                            ctypes.c_double,
                                            ctypes.c_int,
                                            ctypes.c_void_p,
                                            ctypes.c_void_p]

def cuspsolve(A, b):
    Acsr = sp.csr_matrix(A, dtype=float)
    b = np.asarray(b, dtype=float)
    x = np.empty_like(b)

    # Copy arrays to GPU
    dcsrVal = gpuarray.to_gpu(Acsr.data)
    dcsrColInd = gpuarray.to_gpu(Acsr.indices)
    dcsrIndPtr = gpuarray.to_gpu(Acsr.indptr)
    dx = gpuarray.to_gpu(x)
    db = gpuarray.to_gpu(b)

    # Create solver parameters
    m = ctypes.c_int(Acsr.shape[0])  # Need check if A is square
    nnz = ctypes.c_int(Acsr.nnz)
    descrA = ctypes.c_void_p()
    reorder = ctypes.c_int(0)
    tol = ctypes.c_double(1e-10)
    singularity = ctypes.c_int(0)  # -1 if A not singular

    # create cusparse handle
    _cusp_handle = ctypes.c_void_p()
    status = _libcusparse.cusparseCreate(ctypes.byref(_cusp_handle))
    assert(status == 0)
    cusp_handle = _cusp_handle.value

    # create MatDescriptor
    status = _libcusparse.cusparseCreateMatDescr(ctypes.byref(descrA))
    assert(status == 0)

    #create cusolver handle
    _cuso_handle = ctypes.c_void_p()
    status = _libcusolver.cusolverSpCreate(ctypes.byref(_cuso_handle))
    assert(status == 0)
    cuso_handle = _cuso_handle.value

    # Solve
    res=_libcusolver.cusolverSpDcsrlsvqr(cuso_handle,
                                         m,
                                         nnz,
                                         descrA,
                                         int(dcsrVal.gpudata),
                                         int(dcsrIndPtr.gpudata),
                                         int(dcsrColInd.gpudata),
                                         int(db.gpudata),
                                         tol,
                                         reorder,
                                         int(dx.gpudata),
                                         ctypes.byref(singularity))
    assert(res == 0)
    if singularity.value != -1:
        raise ValueError('Singular matrix!')
    x = dx.get()  # Get result as numpy array

    # Destroy handles
    status = _libcusolver.cusolverSpDestroy(cuso_handle)
    assert(status == 0)
    status = _libcusparse.cusparseDestroy(cusp_handle)
    assert(status == 0)

    # Return result
    return x

# Test
if __name__ == '__main__':
    A = np.diag(np.arange(1, 5))
    b = np.ones(4)
    x = cuspsolve(A, b)
    np.testing.assert_almost_equal(x, np.array([1. , 0.5, 0.33333333, 0.25]))

    # Timing comparison
    from scipy.sparse import rand
    from scipy.sparse.linalg import spsolve
    from scipy.sparse import coo_matrix
    import time
    n = 10000
    i = j = np.arange(n)
    diag = np.ones(n)
    A = rand(n, n, density=0.001)
    A = A.tocsr()
    A[i, j] = diag
    b = np.ones(n)

    t0 = time.time()
    x = spsolve(A, b)
    dt1 = time.time() - t0
    print ("scipy.sparse.linalg.spsolve time: %s" %dt1)

    t0 = time.time()
    x = cuspsolve(A, b)
    dt2 = time.time() - t0
    print ("cuspsolve time: %s" %dt2)

    ratio = dt1/dt2
    if ratio > 1:
        print ("CUDA is %s times faster than CPU." %ratio)
    else:
        print ("CUDA is %s times slower than CPU." %(1./ratio))

/root/.pyenv/versions/py365/lib/python3.6/site-packages/scipy/sparse/compressed.py:746: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
  SparseEfficiencyWarning)

scipy.sparse.linalg.spsolve time: 217.97456216812134
cuspsolve time: 117.69748067855835
CUDA is 1.851990041855085 times faster than CPU.