(a_gpu, jobvl='N', jobvr='V', imag='F', lib='cusolver')
Eigendecomposition of a matrix.

    Compute the eigenvalues `w`  for a real/complex square matrix `a`
    and (optionally) the real left and right eigenvectors `vl`, `vr`.

    Parameters
    ----------
    a_gpu : pycuda.gpuarray.GPUArray
        Real/complex input matrix  `a` with dimensions `(m, n)`.
    jobvl :  {'V', 'N'}
        'V' : returns `vl`, the left eigenvectors of `a` with dimensions `(m, m)`.
        'N' : left eigenvectors are not computed.
    jobvr :  {'V', 'N'}
        'V' : returns `vr`, the right eigenvectors of `a` with dimensions
        `(m, m)`, (default).
        'N' : right eigenvectors are not computed.
    imag :  {'F', 'T'}
         'F' : imaginary parts of a real matrix are not returned (default).
         'T' : returns the imaginary parts of a real matrix
         (only relevant in the case of single/double precision ).
    lib : str
        Library to use. May be either 'cula' or 'cusolver' (default). If using
        'cusolver', only symmetric/Hermitian matrices are supported.

    Returns
    -------
    vr_gpu : pycuda.gpuarray.GPUArray
         The normalized (Euclidean norm equal to 1) right eigenvectors,
         such that the column `vr[:,i]` is the eigenvector corresponding
         to the eigenvalue `w[i]`.
    w_gpu : pycuda.gpuarray.GPUArray
        Array containing the real/complex eigenvalues, not necessarily ordered.
        `w` is of length `m`.
    vl_gpu : pycuda.gpuarray.GPUArray
         The normalized (Euclidean norm equal to 1) left eigenvectors,
         such that the column `vl[:,i]` is the eigenvector corresponding
         to the eigenvalue `w[i]`.

    Notes
    -----
    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.

    Arrays are expected to be stored in column-major order, i.e., order='F'.

    Examples
    --------
    >>> import pycuda.gpuarray as gpuarray
    >>> import pycuda.autoinit
    >>> import numpy as np
    >>> from skcuda import linalg
    >>> linalg.init()
    >>> # Compute right eigenvectors of a symmetric matrix A and verify A*vr = vr*w
    >>> a = np.array(([1,3],[3,5]), np.float32, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> vr_gpu, w_gpu = linalg.eig(a_gpu, 'N', 'V')
    >>> np.allclose(np.dot(a, vr_gpu.get()), np.dot(vr_gpu.get(), np.diag(w_gpu.get())), 1e-4)
    True
    >>> # Compute left eigenvectors of a symmetric matrix A and verify vl.T*A = w*vl.T
    >>> a = np.array(([1,3],[3,5]), np.float32, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> w_gpu, vl_gpu = linalg.eig(a_gpu, 'V', 'N')
    >>> np.allclose(np.dot(vl_gpu.get().T, a), np.dot(np.diag(w_gpu.get()), vl_gpu.get().T), 1e-4)
    True
    >>> # Compute left/right eigenvectors of a symmetric matrix A and verify A = vr*w*vl.T
    >>> a = np.array(([1,3],[3,5]), np.float32, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> vr_gpu, w_gpu, vl_gpu = linalg.eig(a_gpu, 'V', 'V')
    >>> np.allclose(a, np.dot(vr_gpu.get(), np.dot(np.diag(w_gpu.get()), vl_gpu.get().T)), 1e-4)
    True
    >>> # Compute eigenvalues of a square matrix A and verify that trace(A)=sum(w)
    >>> a = np.array(np.random.rand(9,9), np.float32, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> w_gpu = linalg.eig(a_gpu, 'N', 'N')
    >>> np.allclose(np.trace(a), sum(w_gpu.get()), 1e-4)
    True
    >>> # Compute eigenvalues of a real valued matrix A possessing complex e-valuesand
    >>> a = np.array(np.array(([1, -2], [1, 3])), np.float32, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> w_gpu = linalg.eig(a_gpu, 'N', 'N', imag='T')
    True
    >>> # Compute eigenvalues of a complex valued matrix A and verify that trace(A)=sum(w)
    >>> a = np.array(np.random.rand(2,2) + 1j*np.random.rand(2,2), np.complex64, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> w_gpu = linalg.eig(a_gpu, 'N', 'N')
    >>> np.allclose(np.trace(a), sum(w_gpu.get()), 1e-4)
    True
Eigendecomposition of a matrix.
Eigendecomposition of a matrix .
def eig(a_gpu, jobvl='N', jobvr='V', imag='F', lib='cusolver'):
    """
    Eigendecomposition of a matrix.

    Compute the eigenvalues `w`  for a real/complex square matrix `a`
    and (optionally) the real left and right eigenvectors `vl`, `vr`.

    Parameters
    ----------
    a_gpu : pycuda.gpuarray.GPUArray
        Real/complex input matrix  `a` with dimensions `(m, n)`.
    jobvl :  {'V', 'N'}
        'V' : returns `vl`, the left eigenvectors of `a` with dimensions `(m, m)`.
        'N' : left eigenvectors are not computed.
    jobvr :  {'V', 'N'}
        'V' : returns `vr`, the right eigenvectors of `a` with dimensions
        `(m, m)`, (default).
        'N' : right eigenvectors are not computed.
    imag :  {'F', 'T'}
         'F' : imaginary parts of a real matrix are not returned (default).
         'T' : returns the imaginary parts of a real matrix
         (only relevant in the case of single/double precision ).
    lib : str
        Library to use. May be either 'cula' or 'cusolver' (default). If using
        'cusolver', only symmetric/Hermitian matrices are supported.

    Returns
    -------
    vr_gpu : pycuda.gpuarray.GPUArray
         The normalized (Euclidean norm equal to 1) right eigenvectors,
         such that the column `vr[:,i]` is the eigenvector corresponding
         to the eigenvalue `w[i]`.
    w_gpu : pycuda.gpuarray.GPUArray
        Array containing the real/complex eigenvalues, not necessarily ordered.
        `w` is of length `m`.
    vl_gpu : pycuda.gpuarray.GPUArray
         The normalized (Euclidean norm equal to 1) left eigenvectors,
         such that the column `vl[:,i]` is the eigenvector corresponding
         to the eigenvalue `w[i]`.

    Notes
    -----
    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.

    Arrays are expected to be stored in column-major order, i.e., order='F'.

    Examples
    --------
    >>> import pycuda.gpuarray as gpuarray
    >>> import pycuda.autoinit
    >>> import numpy as np
    >>> from skcuda import linalg
    >>> linalg.init()
    >>> # Compute right eigenvectors of a symmetric matrix A and verify A*vr = vr*w
    >>> a = np.array(([1,3],[3,5]), np.float32, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> vr_gpu, w_gpu = linalg.eig(a_gpu, 'N', 'V')
    >>> np.allclose(np.dot(a, vr_gpu.get()), np.dot(vr_gpu.get(), np.diag(w_gpu.get())), 1e-4)
    True
    >>> # Compute left eigenvectors of a symmetric matrix A and verify vl.T*A = w*vl.T
    >>> a = np.array(([1,3],[3,5]), np.float32, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> w_gpu, vl_gpu = linalg.eig(a_gpu, 'V', 'N')
    >>> np.allclose(np.dot(vl_gpu.get().T, a), np.dot(np.diag(w_gpu.get()), vl_gpu.get().T), 1e-4)
    True
    >>> # Compute left/right eigenvectors of a symmetric matrix A and verify A = vr*w*vl.T
    >>> a = np.array(([1,3],[3,5]), np.float32, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> vr_gpu, w_gpu, vl_gpu = linalg.eig(a_gpu, 'V', 'V')
    >>> np.allclose(a, np.dot(vr_gpu.get(), np.dot(np.diag(w_gpu.get()), vl_gpu.get().T)), 1e-4)
    True
    >>> # Compute eigenvalues of a square matrix A and verify that trace(A)=sum(w)
    >>> a = np.array(np.random.rand(9,9), np.float32, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> w_gpu = linalg.eig(a_gpu, 'N', 'N')
    >>> np.allclose(np.trace(a), sum(w_gpu.get()), 1e-4)
    True
    >>> # Compute eigenvalues of a real valued matrix A possessing complex e-valuesand
    >>> a = np.array(np.array(([1, -2], [1, 3])), np.float32, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> w_gpu = linalg.eig(a_gpu, 'N', 'N', imag='T')
    True
    >>> # Compute eigenvalues of a complex valued matrix A and verify that trace(A)=sum(w)
    >>> a = np.array(np.random.rand(2,2) + 1j*np.random.rand(2,2), np.complex64, order='F')
    >>> a_gpu = gpuarray.to_gpu(a)
    >>> w_gpu = linalg.eig(a_gpu, 'N', 'N')
    >>> np.allclose(np.trace(a), sum(w_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.culaDeviceCgeev
            imag='F'
        elif data_type == np.float32:
            func = cula.culaDeviceSgeev
        else:
            if cula._libcula_toolkit == 'standard':
                if data_type == np.complex128:
                    func = cula.culaDeviceZgeev
                    imag='F'
                elif data_type == np.float64:
                    func = cula.culaDeviceDgeev
                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

        # FIXME: Seems like CUSOLVER only handles symmetric or Hermitian matrices,
        # look into cusolverDn<t>sygvd
        if data_type == np.complex64:
            func = cusolver.cusolverDnCheevd
            bufsize = cusolver.cusolverDnCheevd_bufferSize
        elif data_type == np.float32:
            func = cusolver.cusolverDnSsyevd
            bufsize = cusolver.cusolverDnSsyevd_bufferSize
        elif data_type == np.complex128:
            func = cusolver.cusolverDnZheevd
            bufsize = cusolver.cusolverDnZheevd_bufferSize
        elif data_type == np.float64:
            real_type = np.float64
            func = cusolver.cusolverDnDsyevd
            bufsize = cusolver.cusolverDnDsyevd_bufferSize
        else:
            raise ValueError('unsupported type')
    else:
        raise ValueError('invalid library specified')

    # CUDA assumes that arrays are stored in column-major order
    n, m = a_gpu.shape

    #Check input
    if(m!=n): raise ValueError('matrix is not square!')
    jobvl = jobvl.upper()
    jobvr = jobvr.upper()

    if jobvl not in ['N', 'V'] :
        raise ValueError('jobvl has to  be "N" or "V" ')
    if jobvr not in ['N', 'V'] :
        raise ValueError('jobvr has to  be "N" or "V" ')
    if imag not in ['T', 'F'] :
        raise ValueError('imag has to  be "T" or "F" ')
    if lib == 'cula':
        w_gpu = gpuarray.empty(m, data_type, order="F", allocator=alloc)
        # Allocate vl, vr, and w:
        vl_gpu = gpuarray.empty((m,m), data_type, order="F", allocator=alloc)
        vr_gpu = gpuarray.empty((m,m), data_type, order="F", allocator=alloc)

        if data_type in (np.complex64, np.complex128):
            #culaDeviceCgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr)
            func(jobvl, jobvr, m, a_gpu.gpudata, m, w_gpu.gpudata, vl_gpu.gpudata , m , vr_gpu.gpudata, m )

        elif data_type in (np.float32, np.float64):
            wi_gpu = gpuarray.zeros(m, data_type, order="F", allocator=alloc)
            func(jobvl, jobvr, m, a_gpu.gpudata, m, w_gpu.gpudata, wi_gpu.gpudata, vl_gpu.gpudata , m , vr_gpu.gpudata, m )

        if imag == 'T':
            w_gpu = w_gpu + (1j)*wi_gpu

        # Free internal CULA memory:
        cula.culaFreeBuffers()

        if jobvl  == 'N' and jobvr == 'N':
            return w_gpu
        elif jobvl == 'V' and jobvr == 'V':
            return vr_gpu, w_gpu, vl_gpu
        elif jobvl == 'V' and jobvr == 'N':
            return w_gpu, vl_gpu,
        elif jobvl == 'N' and jobvr == 'V':
            return vr_gpu, w_gpu
    elif lib == 'cusolver':
        if data_type in (np.float32,np.complex64):
            eigv_data_type = np.float32
        elif data_type in ( np.float64, np.complex128):
            eigv_data_type = np.float64
        w_gpu = gpuarray.empty(m, eigv_data_type, order="F", allocator=alloc)
        if jobvl == 'V':
            raise NotImplementedError('CUSOLVER supports only right eigenvectors')

        if jobvr == 'V':
            jobz = cusolver._CUSOLVER_EIG_MODE['CUSOLVER_EIG_MODE_VECTOR']
            # Copy a_gpu, so we don't destroy it
            a_copy_gpu = a_gpu.copy()
        else:
            jobz = cusolver._CUSOLVER_EIG_MODE['CUSOLVER_EIG_MODE_NOVECTOR']
            a_copy_gpu = a_gpu

        # Since we have the full matrix and assuming symmetry, fill mode
        # hopefully doesn't matter
        uplo = cublas._CUBLAS_FILL_MODE[0]

        Lwork = bufsize(
            cusolverHandle,
            jobz,
            uplo,
            n,
            a_copy_gpu.gpudata,
            m,
            w_gpu.gpudata,
        )

        Work = gpuarray.empty(Lwork, data_type, allocator=alloc)
        devInfo = gpuarray.empty(1, np.int32, allocator=alloc)
        func(cusolverHandle, jobz, uplo,
             n, a_copy_gpu.gpudata, m, w_gpu.gpudata,
             Work.gpudata, Lwork, devInfo.gpudata)

        if jobz == cusolver._CUSOLVER_EIG_MODE['CUSOLVER_EIG_MODE_VECTOR']:
            return a_copy_gpu, w_gpu
        else:
            return w_gpu
    else:
        raise ValueError('invalid library specified')
def eig ( a_gpu , jobvl = 'N' , jobvr = 'V' , imag = 'F' , 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 . culaDeviceCgeev imag = 'F' elif data_type == np . float32 : func = cula . culaDeviceSgeev else : if cula . _libcula_toolkit == 'standard' : if data_type == np . complex128 : func = cula . culaDeviceZgeev imag = 'F' elif data_type == np . float64 : func = cula . culaDeviceDgeev 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 # FIXME: Seems like CUSOLVER only handles symmetric or Hermitian matrices, # look into cusolverDn<t>sygvd if data_type == np . complex64 : func = cusolver . cusolverDnCheevd bufsize = cusolver . cusolverDnCheevd_bufferSize elif data_type == np . float32 : func = cusolver . cusolverDnSsyevd bufsize = cusolver . cusolverDnSsyevd_bufferSize elif data_type == np . complex128 : func = cusolver . cusolverDnZheevd bufsize = cusolver . cusolverDnZheevd_bufferSize elif data_type == np . float64 : real_type = np . float64 func = cusolver . cusolverDnDsyevd bufsize = cusolver . cusolverDnDsyevd_bufferSize else : raise ValueError ( 'unsupported type' ) else : raise ValueError ( 'invalid library specified' ) # CUDA assumes that arrays are stored in column-major order n , m = a_gpu . shape #Check input if ( m != n ) : raise ValueError ( 'matrix is not square!' ) jobvl = jobvl . upper ( ) jobvr = jobvr . upper ( ) if jobvl not in [ 'N' , 'V' ] : raise ValueError ( 'jobvl has to  be "N" or "V" ' ) if jobvr not in [ 'N' , 'V' ] : raise ValueError ( 'jobvr has to  be "N" or "V" ' ) if imag not in [ 'T' , 'F' ] : raise ValueError ( 'imag has to  be "T" or "F" ' ) if lib == 'cula' : w_gpu = gpuarray . empty ( m , data_type , order = "F" , allocator = alloc ) # Allocate vl, vr, and w: vl_gpu = gpuarray . empty ( ( m , m ) , data_type , order = "F" , allocator = alloc ) vr_gpu = gpuarray . empty ( ( m , m ) , data_type , order = "F" , allocator = alloc ) if data_type in ( np . complex64 , np . complex128 ) : #culaDeviceCgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr) func ( jobvl , jobvr , m , a_gpu . gpudata , m , w_gpu . gpudata , vl_gpu . gpudata , m , vr_gpu . gpudata , m ) elif data_type in ( np . float32 , np . float64 ) : wi_gpu = gpuarray . zeros ( m , data_type , order = "F" , allocator = alloc ) func ( jobvl , jobvr , m , a_gpu . gpudata , m , w_gpu . gpudata , wi_gpu . gpudata , vl_gpu . gpudata , m , vr_gpu . gpudata , m ) if imag == 'T' : w_gpu = w_gpu + ( 1j ) * wi_gpu # Free internal CULA memory: cula . culaFreeBuffers ( ) if jobvl == 'N' and jobvr == 'N' : return w_gpu elif jobvl == 'V' and jobvr == 'V' : return vr_gpu , w_gpu , vl_gpu elif jobvl == 'V' and jobvr == 'N' : return w_gpu , vl_gpu , elif jobvl == 'N' and jobvr == 'V' : return vr_gpu , w_gpu elif lib == 'cusolver' : if data_type in ( np . float32 , np . complex64 ) : eigv_data_type = np . float32 elif data_type in ( np . float64 , np . complex128 ) : eigv_data_type = np . float64 w_gpu = gpuarray . empty ( m , eigv_data_type , order = "F" , allocator = alloc ) if jobvl == 'V' : raise NotImplementedError ( 'CUSOLVER supports only right eigenvectors' ) if jobvr == 'V' : jobz = cusolver . _CUSOLVER_EIG_MODE [ 'CUSOLVER_EIG_MODE_VECTOR' ] # Copy a_gpu, so we don't destroy it a_copy_gpu = a_gpu . copy ( ) else : jobz = cusolver . _CUSOLVER_EIG_MODE [ 'CUSOLVER_EIG_MODE_NOVECTOR' ] a_copy_gpu = a_gpu # Since we have the full matrix and assuming symmetry, fill mode # hopefully doesn't matter uplo = cublas . _CUBLAS_FILL_MODE [ 0 ] Lwork = bufsize ( cusolverHandle , jobz , uplo , n , a_copy_gpu . gpudata , m , w_gpu . gpudata , ) Work = gpuarray . empty ( Lwork , data_type , allocator = alloc ) devInfo = gpuarray . empty ( 1 , np . int32 , allocator = alloc ) func ( cusolverHandle , jobz , uplo , n , a_copy_gpu . gpudata , m , w_gpu . gpudata , Work . gpudata , Lwork , devInfo . gpudata ) if jobz == cusolver . _CUSOLVER_EIG_MODE [ 'CUSOLVER_EIG_MODE_VECTOR' ] : return a_copy_gpu , w_gpu else : return w_gpu else : raise ValueError ( 'invalid library specified' )