(a_gpu, b_gpu, uplo='L', lib='cusolver')
Cholesky solver.

    Solve a system of equations via Cholesky factorization,
    i.e. `a*x = b`.
    Overwrites `b` to give `inv(a)*b`, and overwrites the chosen triangle
    of `a` with factorized triangle.

    Parameters
    ----------
    a : pycuda.gpuarray.GPUArray
        Input matrix of shape `(n, n)` to decompose.
    b : pycuda.gpuarray.GPUArray
        Input matrix of shape `(n, m)` to decompose. If m > 1,
        `a` and `b` must be stored in Fortran order.
    uplo: chr
        Use the upper='U' or lower='L' (default) triangle of `a`.
    lib : str
        Library to use. May be either 'cula' or 'cusolver' (default).

    Notes
    -----
    If using CULA, double precision is only supported if the standard version of the
    CULA Dense toolkit is installed.

    Examples
    --------
    >>> import pycuda.gpuarray as gpuarray
    >>> import pycuda.autoinit
    >>> import numpy as np
    >>> import scipy.linalg
    >>> import skcuda.linalg as linalg
    >>> linalg.init()
    >>> a = np.array([[3, 0], [0, 7]]).asarray(np.float64)
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> b = np.array([11, 19]).astype(np.float64)
    >>> b_gpu  = gpuarray.to_gpu(b)
    >>> cho_solve(a_gpu, b_gpu)
    >>> np.allclose(b_gpu.get(), scipy.linalg.cho_solve(scipy.linalg.cho_factor(a), b))
    True
Cholesky solver.
Cholesky solver .
def cho_solve(a_gpu, b_gpu, uplo='L', lib='cusolver'):
    """
    Cholesky solver.

    Solve a system of equations via Cholesky factorization,
    i.e. `a*x = b`.
    Overwrites `b` to give `inv(a)*b`, and overwrites the chosen triangle
    of `a` with factorized triangle.

    Parameters
    ----------
    a : pycuda.gpuarray.GPUArray
        Input matrix of shape `(n, n)` to decompose.
    b : pycuda.gpuarray.GPUArray
        Input matrix of shape `(n, m)` to decompose. If m > 1,
        `a` and `b` must be stored in Fortran order.
    uplo: chr
        Use the upper='U' or lower='L' (default) triangle of `a`.
    lib : str
        Library to use. May be either 'cula' or 'cusolver' (default).

    Notes
    -----
    If using CULA, double precision is only supported if the standard version of the
    CULA Dense toolkit is installed.

    Examples
    --------
    >>> import pycuda.gpuarray as gpuarray
    >>> import pycuda.autoinit
    >>> import numpy as np
    >>> import scipy.linalg
    >>> import skcuda.linalg as linalg
    >>> linalg.init()
    >>> a = np.array([[3, 0], [0, 7]]).asarray(np.float64)
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> b = np.array([11, 19]).astype(np.float64)
    >>> b_gpu  = gpuarray.to_gpu(b)
    >>> cho_solve(a_gpu, b_gpu)
    >>> np.allclose(b_gpu.get(), scipy.linalg.cho_solve(scipy.linalg.cho_factor(a), b))
    True
    """

    alloc = misc._global_cublas_allocator

    data_type = a_gpu.dtype.type
    if lib == 'cula':
        if not _has_cula:
            raise NotImplementedError('CULA not installed')

        if cula._libcula_toolkit == 'standard':
            if data_type == np.complex64:
                func = cula.culaDeviceCposv
            elif data_type == np.float32:
                func = cula.culaDeviceSposv
            elif data_type == np.complex128:
                func = cula.culaDeviceZposv
            elif data_type == np.float64:
                func = cula.culaDeviceDposv
            else:
                raise ValueError('unsupported type')
        else:
            raise ValueError('Cholesky factorization not included in CULA Dense Free version')

    elif lib == 'cusolver':
        if not _has_cusolver:
            raise NotImplementedError('CUSOLVER not installed')
        cusolverHandle = misc._global_cusolver_handle

        if data_type == np.complex64:
            func = cusolver.cusolverDnCpotrs
        elif data_type == np.float32:
            func = cusolver.cusolverDnSpotrs
        elif data_type == np.complex128:
            func = cusolver.cusolverDnZpotrs
        elif data_type == np.float64:
            func = cusolver.cusolverDnDpotrs
        else:
            raise ValueError('unsupported type')

    else:
        raise ValueError('invalid library specified')

    # Since CUDA assumes that arrays are stored in column-major
    # format, the input matrix is assumed to be transposed:
    na, ma = a_gpu.shape

    if (na!=ma):
        raise ValueError('Matrix must be symmetric positive-definite')

    if a_gpu.flags.c_contiguous != b_gpu.flags.c_contiguous:
        raise ValueError('unsupported combination of input order')

    b_shape = b_gpu.shape
    if len(b_shape) == 1:
        b_shape = (b_shape[0], 1)

    if a_gpu.flags.f_contiguous:
        lda = max(1, na)
        ldb = max(1, b_shape[0])
    else:
        lda = max(1, ma)
        ldb = lda
        if b_shape[1] > 1:
            raise ValueError('only vectors allowed in c-order RHS')

    if lib == 'cula':
        # Assuming we are only solving for a vector. Hence, nrhs = 1
        func(uplo, na, b_shape[1], int(a_gpu.gpudata), lda,
             int(b_gpu.gpudata), ldb)

        # Free internal CULA memory:
        cula.culaFreeBuffers()
    else:
        # CUSOLVER expects uplo to be an int rather than a char:
        uplo = cublas._CUBLAS_FILL_MODE[uplo]

        # Since CUSOLVER doesn't implement POSV as of 8.0, we need to factor the
        # given matrix before calling POTRS:
        cho_factor(a_gpu, uplo, lib)

        # Assuming we are only solving for a vector. Hence, nrhs = 1
        devInfo = gpuarray.empty(1, np.int32, allocator=alloc)
        func(cusolverHandle, uplo, na, b_shape[1], int(a_gpu.gpudata), lda,
             int(b_gpu.gpudata), ldb, int(devInfo.gpudata))
def cho_solve ( a_gpu , b_gpu , uplo = 'L' , lib = 'cusolver' ) : alloc = misc . _global_cublas_allocator data_type = a_gpu . dtype . type if lib == 'cula' : if not _has_cula : raise NotImplementedError ( 'CULA not installed' ) if cula . _libcula_toolkit == 'standard' : if data_type == np . complex64 : func = cula . culaDeviceCposv elif data_type == np . float32 : func = cula . culaDeviceSposv elif data_type == np . complex128 : func = cula . culaDeviceZposv elif data_type == np . float64 : func = cula . culaDeviceDposv else : raise ValueError ( 'unsupported type' ) else : raise ValueError ( 'Cholesky factorization not included in CULA Dense Free version' ) elif lib == 'cusolver' : if not _has_cusolver : raise NotImplementedError ( 'CUSOLVER not installed' ) cusolverHandle = misc . _global_cusolver_handle if data_type == np . complex64 : func = cusolver . cusolverDnCpotrs elif data_type == np . float32 : func = cusolver . cusolverDnSpotrs elif data_type == np . complex128 : func = cusolver . cusolverDnZpotrs elif data_type == np . float64 : func = cusolver . cusolverDnDpotrs else : raise ValueError ( 'unsupported type' ) else : raise ValueError ( 'invalid library specified' ) # Since CUDA assumes that arrays are stored in column-major # format, the input matrix is assumed to be transposed: na , ma = a_gpu . shape if ( na != ma ) : raise ValueError ( 'Matrix must be symmetric positive-definite' ) if a_gpu . flags . c_contiguous != b_gpu . flags . c_contiguous : raise ValueError ( 'unsupported combination of input order' ) b_shape = b_gpu . shape if len ( b_shape ) == 1 : b_shape = ( b_shape [ 0 ] , 1 ) if a_gpu . flags . f_contiguous : lda = max ( 1 , na ) ldb = max ( 1 , b_shape [ 0 ] ) else : lda = max ( 1 , ma ) ldb = lda if b_shape [ 1 ] > 1 : raise ValueError ( 'only vectors allowed in c-order RHS' ) if lib == 'cula' : # Assuming we are only solving for a vector. Hence, nrhs = 1 func ( uplo , na , b_shape [ 1 ] , int ( a_gpu . gpudata ) , lda , int ( b_gpu . gpudata ) , ldb ) # Free internal CULA memory: cula . culaFreeBuffers ( ) else : # CUSOLVER expects uplo to be an int rather than a char: uplo = cublas . _CUBLAS_FILL_MODE [ uplo ] # Since CUSOLVER doesn't implement POSV as of 8.0, we need to factor the # given matrix before calling POTRS: cho_factor ( a_gpu , uplo , lib ) # Assuming we are only solving for a vector. Hence, nrhs = 1 devInfo = gpuarray . empty ( 1 , np . int32 , allocator = alloc ) func ( cusolverHandle , uplo , na , b_shape [ 1 ] , int ( a_gpu . gpudata ) , lda , int ( b_gpu . gpudata ) , ldb , int ( devInfo . gpudata ) )