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

    Performs an in-place Cholesky factorization on the matrix `a`
    such that `a = x*x.T` or `x.T*x`, if the lower='L' or upper='U'
    triangle of `a` is used, respectively.

    Parameters
    ----------
    a_gpu : pycuda.gpuarray.GPUArray
        Input matrix of shape `(m, m)` to decompose.
    uplo : {'U', 'L'}
        Use upper or lower (default) triangle of 'a_gpu'
    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.0],[0.0,7.0]])
    >>> a = np.asarray(a, np.float64)
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> cho_factor(a_gpu)
    >>> np.allclose(a_gpu.get(), scipy.linalg.cho_factor(a)[0])
    True
Cholesky factorization.
Cholesky factorization .
def cho_factor(a_gpu, uplo='L', lib='cusolver'):
    """
    Cholesky factorization.

    Performs an in-place Cholesky factorization on the matrix `a`
    such that `a = x*x.T` or `x.T*x`, if the lower='L' or upper='U'
    triangle of `a` is used, respectively.

    Parameters
    ----------
    a_gpu : pycuda.gpuarray.GPUArray
        Input matrix of shape `(m, m)` to decompose.
    uplo : {'U', 'L'}
        Use upper or lower (default) triangle of 'a_gpu'
    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.0],[0.0,7.0]])
    >>> a = np.asarray(a, np.float64)
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> cho_factor(a_gpu)
    >>> np.allclose(a_gpu.get(), scipy.linalg.cho_factor(a)[0])
    True
    """

    alloc = misc._global_cublas_allocator

    data_type = a_gpu.dtype.type

    if lib == 'cula':
        if not _has_cula:
            raise NotImplementedError('CULA not installed')

        real_type = np.float32
        if cula._libcula_toolkit == 'standard':
            if data_type == np.complex64:
                func = cula.culaDeviceCpotrf
            elif data_type == np.float32:
                func = cula.culaDeviceSpotrf
            elif data_type == np.complex128:
                func = cula.culaDeviceZpotrf
            elif data_type == np.float64:
                func = cula.culaDeviceDpotrf
            else:
                raise ValueError('unsupported type')
            real_type = np.float64
        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.cusolverDnCpotrf
            bufsize = cusolver.cusolverDnCpotrf_bufferSize
        elif data_type == np.float32:
            func = cusolver.cusolverDnSpotrf
            bufsize = cusolver.cusolverDnSpotrf_bufferSize
        elif data_type == np.complex128:
            real_type = np.float64
            func = cusolver.cusolverDnZpotrf
            bufsize = cusolver.cusolverDnZpotrf_bufferSize
        elif data_type == np.float64:
            real_type = np.float64
            func = cusolver.cusolverDnDpotrf
            bufsize = cusolver.cusolverDnDpotrf_bufferSize
        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:
    n, m = a_gpu.shape
    if (n!=m):
        raise ValueError('Matrix must be symmetric positive-definite')

    # Set the leading dimension of the input matrix:
    lda = max(1, m)

    # Factorize and check error status:
    if lib == 'cula':
        func(uplo, n, int(a_gpu.gpudata), lda)

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

        # Allocate working space:
        Lwork = bufsize(misc._global_cusolver_handle, uplo, n, int(a_gpu.gpudata), lda)
        Work = gpuarray.empty(Lwork, data_type, allocator=alloc)
        devInfo = gpuarray.empty(1, np.int32, allocator=alloc)

        func(misc._global_cusolver_handle, uplo, n, int(a_gpu.gpudata), lda,
             int(Work.gpudata), Lwork, int(devInfo.gpudata))

        # Free working space:
        del Work, devInfo
def cho_factor ( a_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' ) real_type = np . float32 if cula . _libcula_toolkit == 'standard' : if data_type == np . complex64 : func = cula . culaDeviceCpotrf elif data_type == np . float32 : func = cula . culaDeviceSpotrf elif data_type == np . complex128 : func = cula . culaDeviceZpotrf elif data_type == np . float64 : func = cula . culaDeviceDpotrf else : raise ValueError ( 'unsupported type' ) real_type = np . float64 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 . cusolverDnCpotrf bufsize = cusolver . cusolverDnCpotrf_bufferSize elif data_type == np . float32 : func = cusolver . cusolverDnSpotrf bufsize = cusolver . cusolverDnSpotrf_bufferSize elif data_type == np . complex128 : real_type = np . float64 func = cusolver . cusolverDnZpotrf bufsize = cusolver . cusolverDnZpotrf_bufferSize elif data_type == np . float64 : real_type = np . float64 func = cusolver . cusolverDnDpotrf bufsize = cusolver . cusolverDnDpotrf_bufferSize 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: n , m = a_gpu . shape if ( n != m ) : raise ValueError ( 'Matrix must be symmetric positive-definite' ) # Set the leading dimension of the input matrix: lda = max ( 1 , m ) # Factorize and check error status: if lib == 'cula' : func ( uplo , n , int ( a_gpu . gpudata ) , lda ) # Free internal CULA memory: cula . culaFreeBuffers ( ) else : # CUSOLVER expects uplo to be an int rather than a char: uplo = cublas . _CUBLAS_FILL_MODE [ uplo ] # Allocate working space: Lwork = bufsize ( misc . _global_cusolver_handle , uplo , n , int ( a_gpu . gpudata ) , lda ) Work = gpuarray . empty ( Lwork , data_type , allocator = alloc ) devInfo = gpuarray . empty ( 1 , np . int32 , allocator = alloc ) func ( misc . _global_cusolver_handle , uplo , n , int ( a_gpu . gpudata ) , lda , int ( Work . gpudata ) , Lwork , int ( devInfo . gpudata ) ) # Free working space: del Work , devInfo