(a_gpu, jobu='A', jobvt='A', lib='cusolver')
Singular Value Decomposition.

    Factors the matrix `a` into two unitary matrices, `u` and `vh`,
    and a 1-dimensional array of real, non-negative singular values,
    `s`, such that `a == dot(u.T, dot(diag(s), vh.T))`.

    Parameters
    ----------
    a : pycuda.gpuarray.GPUArray
        Input matrix of shape `(m, n)` to decompose.
    jobu : {'A', 'S', 'O', 'N'}
        If 'A', return the full `u` matrix with shape `(m, m)`.
        If 'S', return the `u` matrix with shape `(m, k)`.
        If 'O', return the `u` matrix with shape `(m, k) without
        allocating a new matrix.
        If 'N', don't return `u`.
    jobvt : {'A', 'S', 'O', 'N'}
        If 'A', return the full `vh` matrix with shape `(n, n)`.
        If 'S', return the `vh` matrix with shape `(k, n)`.
        If 'O', return the `vh` matrix with shape `(k, n) without
        allocating a new matrix.
        If 'N', don't return `vh`.
    lib : str
        Library to use. May be either 'cula' or 'cusolver' (default).

    Returns
    -------
    u : pycuda.gpuarray.GPUArray
        Unitary matrix of shape `(m, m)` or `(m, k)` depending on
        value of `jobu`.
    s : pycuda.gpuarray.GPUArray
        Array containing the singular values, sorted such that `s[i] >= s[i+1]`.
        `s` is of length `min(m, n)`.
    vh : pycuda.gpuarray.GPUArray
        Unitary matrix of shape `(n, n)` or `(k, n)`, depending
        on `jobvt`.

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

    This function destroys the contents of the input matrix regardless
    of the values of `jobu` and `jobvt`.

    Only one of `jobu` or `jobvt` may be set to `O`, and then only for
    a square matrix.

    The CUSOLVER library in CUDA 7.0 only supports `jobu` == `jobvt` == 'A'.

    Examples
    --------
    >>> import pycuda.autoinit
    >>> import pycuda.gpuarray as gpuarray
    >>> import numpy as np
    >>> import skcuda.linalg as linalg
    >>> linalg.init()
    >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
    >>> a = np.asarray(a, np.complex64)
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> u_gpu, s_gpu, vh_gpu = linalg.svd(a_gpu, 'S', 'S')
    >>> np.allclose(a, np.dot(u_gpu.get(), np.dot(np.diag(s_gpu.get()), vh_gpu.get())), 1e-4)
    True
Singular Value Decomposition.
Singular Value Decomposition .
def svd(a_gpu, jobu='A', jobvt='A', lib='cusolver'):
    """
    Singular Value Decomposition.

    Factors the matrix `a` into two unitary matrices, `u` and `vh`,
    and a 1-dimensional array of real, non-negative singular values,
    `s`, such that `a == dot(u.T, dot(diag(s), vh.T))`.

    Parameters
    ----------
    a : pycuda.gpuarray.GPUArray
        Input matrix of shape `(m, n)` to decompose.
    jobu : {'A', 'S', 'O', 'N'}
        If 'A', return the full `u` matrix with shape `(m, m)`.
        If 'S', return the `u` matrix with shape `(m, k)`.
        If 'O', return the `u` matrix with shape `(m, k) without
        allocating a new matrix.
        If 'N', don't return `u`.
    jobvt : {'A', 'S', 'O', 'N'}
        If 'A', return the full `vh` matrix with shape `(n, n)`.
        If 'S', return the `vh` matrix with shape `(k, n)`.
        If 'O', return the `vh` matrix with shape `(k, n) without
        allocating a new matrix.
        If 'N', don't return `vh`.
    lib : str
        Library to use. May be either 'cula' or 'cusolver' (default).

    Returns
    -------
    u : pycuda.gpuarray.GPUArray
        Unitary matrix of shape `(m, m)` or `(m, k)` depending on
        value of `jobu`.
    s : pycuda.gpuarray.GPUArray
        Array containing the singular values, sorted such that `s[i] >= s[i+1]`.
        `s` is of length `min(m, n)`.
    vh : pycuda.gpuarray.GPUArray
        Unitary matrix of shape `(n, n)` or `(k, n)`, depending
        on `jobvt`.

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

    This function destroys the contents of the input matrix regardless
    of the values of `jobu` and `jobvt`.

    Only one of `jobu` or `jobvt` may be set to `O`, and then only for
    a square matrix.

    The CUSOLVER library in CUDA 7.0 only supports `jobu` == `jobvt` == 'A'.

    Examples
    --------
    >>> import pycuda.autoinit
    >>> import pycuda.gpuarray as gpuarray
    >>> import numpy as np
    >>> import skcuda.linalg as linalg
    >>> linalg.init()
    >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
    >>> a = np.asarray(a, np.complex64)
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> u_gpu, s_gpu, vh_gpu = linalg.svd(a_gpu, 'S', 'S')
    >>> np.allclose(a, np.dot(u_gpu.get(), np.dot(np.diag(s_gpu.get()), vh_gpu.get())), 1e-4)
    True
    """

    alloc = misc._global_cublas_allocator

    # The free version of CULA only supports single precision floating
    # point numbers:
    data_type = a_gpu.dtype.type
    real_type = np.float32

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

        if data_type == np.complex64:
            func = cula.culaDeviceCgesvd
        elif data_type == np.float32:
            func = cula.culaDeviceSgesvd
        else:
            if cula._libcula_toolkit == 'standard':
                if data_type == np.complex128:
                    func = cula.culaDeviceZgesvd
                elif data_type == np.float64:
                    func = cula.culaDeviceDgesvd
                else:
                    raise ValueError('unsupported type')
                real_type = np.float64
            else:
                raise ValueError('double precision not supported')
    elif lib == 'cusolver':
        if not _has_cusolver:
            raise NotImplementedError('CUSOLVER not installed')

        cusolverHandle = misc._global_cusolver_handle

        if data_type == np.complex64:
            func = cusolver.cusolverDnCgesvd
            bufsize = cusolver.cusolverDnCgesvd_bufferSize
        elif data_type == np.float32:
            func = cusolver.cusolverDnSgesvd
            bufsize = cusolver.cusolverDnSgesvd_bufferSize
        elif data_type == np.complex128:
            real_type = np.float64
            func = cusolver.cusolverDnZgesvd
            bufsize = cusolver.cusolverDnZgesvd_bufferSize
        elif data_type == np.float64:
            real_type = np.float64
            func = cusolver.cusolverDnDgesvd
            bufsize = cusolver.cusolverDnDgesvd_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
    square = (n == m)

    # CUSOLVER's gesvd routines only support m >= n as of CUDA 7.5:
    if lib == 'cusolver' and m < n:
        raise ValueError('CUSOLVER only supports a_gpu.shape[1] >= a_gpu.shape[0]')

    # Since the input matrix is transposed, jobu and jobvt must also
    # be switched because the computed matrices will be returned in
    # reversed order:
    jobvt, jobu = jobu, jobvt

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

    # Allocate the array of singular values:
    s_gpu = gpuarray.empty(min(m, n), real_type, allocator=alloc)

    # CUSOLVER in CUDA 7.0 only supports jobu = jobvt = 'A':
    jobu = jobu.upper()
    jobvt = jobvt.upper()
    if lib == 'cusolver' and (jobu != 'A' or jobvt != 'A') and \
      cudart._cudart_version <= 7000:
         raise ValueError("CUSOLVER 7.0 only supports jobu = jobvt = 'A'")

    # Set the leading dimension and allocate u:
    ldu = m
    if jobu == 'A':
        u_gpu = gpuarray.empty((ldu, m), data_type, allocator=alloc)
    elif jobu == 'S':
        u_gpu = gpuarray.empty((min(m, n), ldu), data_type, allocator=alloc)
    elif jobu == 'O':
        if not square:
            raise ValueError('in-place computation of singular vectors '+
                             'of non-square matrix not allowed')
        ldu = a_gpu.shape[1]
        u_gpu = a_gpu
    else:
        ldu = 1
        u_gpu = gpuarray.empty((), data_type, allocator=alloc)

    # Set the leading dimension and allocate vh:
    if jobvt == 'A':
        ldvt = n
        vh_gpu = gpuarray.empty((n, n), data_type, allocator=alloc)
    elif jobvt == 'S':
        ldvt = min(m, n)
        vh_gpu = gpuarray.empty((n, ldvt), data_type, allocator=alloc)
    elif jobvt == 'O':
        if jobu == 'O':
            raise ValueError('jobu and jobvt cannot both be O')
        if not square:
            raise ValueError('in-place computation of singular vectors '+
                             'of non-square matrix not allowed')
        ldvt = a_gpu.shape[1]
        vh_gpu = a_gpu
    else:
        ldvt = 1
        vh_gpu = gpuarray.empty((), data_type, allocator=alloc)

    # Compute SVD and check error status:
    if lib == 'cula':
        func(jobu, jobvt, m, n, int(a_gpu.gpudata),
             lda, int(s_gpu.gpudata), int(u_gpu.gpudata),
             ldu, int(vh_gpu.gpudata), ldvt)

        # Free internal CULA memory:
        cula.culaFreeBuffers()
    else:
        # Allocate working space:
        Lwork = bufsize(misc._global_cusolver_handle, m, n)

        Work = gpuarray.empty(Lwork, data_type, allocator=alloc)
        devInfo = gpuarray.empty(1, np.int32, allocator=alloc)

        # rwork is only needed for complex arrays:
        if data_type != real_type:
            rwork = np.empty(Lwork, real_type).ctypes.data
        else:
            rwork = 0
        func(misc._global_cusolver_handle,
             jobu, jobvt, m, n, int(a_gpu.gpudata),
             lda, int(s_gpu.gpudata), int(u_gpu.gpudata),
             ldu, int(vh_gpu.gpudata), ldvt,
             int(Work.gpudata), Lwork, rwork,
             int(devInfo.gpudata))

        # Free working space:
        del rwork, Work, devInfo

    # Since the input is assumed to be transposed, it is necessary to
    # return the computed matrices in reverse order:
    if jobu in ['A', 'S', 'O'] and jobvt in ['A', 'S', 'O']:
        return vh_gpu, s_gpu, u_gpu
    elif jobu == 'N' and jobvt != 'N':
        return vh_gpu, s_gpu
    elif jobu != 'N' and jobvt == 'N':
        return s_gpu, u_gpu
    else:
        return s_gpu
def svd ( a_gpu , jobu = 'A' , jobvt = 'A' , lib = 'cusolver' ) : alloc = misc . _global_cublas_allocator # The free version of CULA only supports single precision floating # point numbers: data_type = a_gpu . dtype . type real_type = np . float32 if lib == 'cula' : if not _has_cula : raise NotImplementedError ( 'CULA not installed' ) if data_type == np . complex64 : func = cula . culaDeviceCgesvd elif data_type == np . float32 : func = cula . culaDeviceSgesvd else : if cula . _libcula_toolkit == 'standard' : if data_type == np . complex128 : func = cula . culaDeviceZgesvd elif data_type == np . float64 : func = cula . culaDeviceDgesvd else : raise ValueError ( 'unsupported type' ) real_type = np . float64 else : raise ValueError ( 'double precision not supported' ) elif lib == 'cusolver' : if not _has_cusolver : raise NotImplementedError ( 'CUSOLVER not installed' ) cusolverHandle = misc . _global_cusolver_handle if data_type == np . complex64 : func = cusolver . cusolverDnCgesvd bufsize = cusolver . cusolverDnCgesvd_bufferSize elif data_type == np . float32 : func = cusolver . cusolverDnSgesvd bufsize = cusolver . cusolverDnSgesvd_bufferSize elif data_type == np . complex128 : real_type = np . float64 func = cusolver . cusolverDnZgesvd bufsize = cusolver . cusolverDnZgesvd_bufferSize elif data_type == np . float64 : real_type = np . float64 func = cusolver . cusolverDnDgesvd bufsize = cusolver . cusolverDnDgesvd_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 square = ( n == m ) # CUSOLVER's gesvd routines only support m >= n as of CUDA 7.5: if lib == 'cusolver' and m < n : raise ValueError ( 'CUSOLVER only supports a_gpu.shape[1] >= a_gpu.shape[0]' ) # Since the input matrix is transposed, jobu and jobvt must also # be switched because the computed matrices will be returned in # reversed order: jobvt , jobu = jobu , jobvt # Set the leading dimension of the input matrix: lda = max ( 1 , m ) # Allocate the array of singular values: s_gpu = gpuarray . empty ( min ( m , n ) , real_type , allocator = alloc ) # CUSOLVER in CUDA 7.0 only supports jobu = jobvt = 'A': jobu = jobu . upper ( ) jobvt = jobvt . upper ( ) if lib == 'cusolver' and ( jobu != 'A' or jobvt != 'A' ) and cudart . _cudart_version <= 7000 : raise ValueError ( "CUSOLVER 7.0 only supports jobu = jobvt = 'A'" ) # Set the leading dimension and allocate u: ldu = m if jobu == 'A' : u_gpu = gpuarray . empty ( ( ldu , m ) , data_type , allocator = alloc ) elif jobu == 'S' : u_gpu = gpuarray . empty ( ( min ( m , n ) , ldu ) , data_type , allocator = alloc ) elif jobu == 'O' : if not square : raise ValueError ( 'in-place computation of singular vectors ' + 'of non-square matrix not allowed' ) ldu = a_gpu . shape [ 1 ] u_gpu = a_gpu else : ldu = 1 u_gpu = gpuarray . empty ( ( ) , data_type , allocator = alloc ) # Set the leading dimension and allocate vh: if jobvt == 'A' : ldvt = n vh_gpu = gpuarray . empty ( ( n , n ) , data_type , allocator = alloc ) elif jobvt == 'S' : ldvt = min ( m , n ) vh_gpu = gpuarray . empty ( ( n , ldvt ) , data_type , allocator = alloc ) elif jobvt == 'O' : if jobu == 'O' : raise ValueError ( 'jobu and jobvt cannot both be O' ) if not square : raise ValueError ( 'in-place computation of singular vectors ' + 'of non-square matrix not allowed' ) ldvt = a_gpu . shape [ 1 ] vh_gpu = a_gpu else : ldvt = 1 vh_gpu = gpuarray . empty ( ( ) , data_type , allocator = alloc ) # Compute SVD and check error status: if lib == 'cula' : func ( jobu , jobvt , m , n , int ( a_gpu . gpudata ) , lda , int ( s_gpu . gpudata ) , int ( u_gpu . gpudata ) , ldu , int ( vh_gpu . gpudata ) , ldvt ) # Free internal CULA memory: cula . culaFreeBuffers ( ) else : # Allocate working space: Lwork = bufsize ( misc . _global_cusolver_handle , m , n ) Work = gpuarray . empty ( Lwork , data_type , allocator = alloc ) devInfo = gpuarray . empty ( 1 , np . int32 , allocator = alloc ) # rwork is only needed for complex arrays: if data_type != real_type : rwork = np . empty ( Lwork , real_type ) . ctypes . data else : rwork = 0 func ( misc . _global_cusolver_handle , jobu , jobvt , m , n , int ( a_gpu . gpudata ) , lda , int ( s_gpu . gpudata ) , int ( u_gpu . gpudata ) , ldu , int ( vh_gpu . gpudata ) , ldvt , int ( Work . gpudata ) , Lwork , rwork , int ( devInfo . gpudata ) ) # Free working space: del rwork , Work , devInfo # Since the input is assumed to be transposed, it is necessary to # return the computed matrices in reverse order: if jobu in [ 'A' , 'S' , 'O' ] and jobvt in [ 'A' , 'S' , 'O' ] : return vh_gpu , s_gpu , u_gpu elif jobu == 'N' and jobvt != 'N' : return vh_gpu , s_gpu elif jobu != 'N' and jobvt == 'N' : return s_gpu , u_gpu else : return s_gpu