(a_gpu, mode='reduced', handle=None, lib='cusolver')
QR Decomposition.

    Factor the real/complex matrix `a` as `QR`, where `Q` is an orthonormal/unitary
    matrix and `R` is an upper triangular matrix.

    Parameters
    ----------
    a_gpu: pycuda.gpuarray.GPUArray
        Real/complex input matrix  `a` with dimensions `(m, n)`.
        `a` is assumed to be `m`>=`n`.
    mode :  {'reduced', 'economic', 'r'}
        'reduced' : returns `Q`, `R` with dimensions `(m, k)` and `(k, n)` (default).
        'economic' : returns `Q` only with dimensions `(m, k)`.
        'r' : returns `R` only with dimensions `(k, n)` with `k`=min`(m,n)`.
    handle : int
        CUBLAS context. If no context is specified, the default handle from
        `skcuda.misc._global_cublas_handle` is used.
    lib : str
        Library to use. May be either 'cula' or 'cusolver' (default).

    Returns
    -------
    q_gpu : pycuda.gpuarray.GPUArray
        Orthonormal/unitary matrix (depending on whether or not `A` is real/complex).
    r_gpu : pycuda.gpuarray.GPUArray
        The upper-triangular matrix.

    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 assumed to be stored in column-major order, i.e., order='F'.

    Examples
    --------
    >>> import pycuda.autoinit
    >>> import pycuda.gpuarray as gpuarray
    >>> import numpy as np
    >>> import skcuda.linalg as linalg
    >>> linalg.init()
    >>> # Rectangular matrix A, np.float32
    >>> A = np.array(np.random.randn(9, 7), np.float32, order='F')
    >>> A_gpu = gpuarray.to_gpu(A)
    >>> Q_gpu, R_gpu = linalg.qr(A_gpu, 'reduced')
    >>> np.allclose(A, np.dot(Q_gpu.get(), R_gpu.get()), 1e-4)
    True
    >>> # Square matrix A, np.complex128
    >>> A = np.random.randn(9, 9) + 1j*np.random.randn(9, 9)
    >>> A = np.asarray(A, np.complex128, order='F')
    >>> A_gpu = gpuarray.to_gpu(A)
    >>> Q_gpu, R_gpu = linalg.qr(A_gpu, 'reduced')
    >>> np.allclose(A, np.dot(Q_gpu.get(), R_gpu.get()), 1e-4)
    True
    >>> np.allclose(np.identity(Q_gpu.shape[0]) + 1j*0, np.dot(Q_gpu.get().conj().T, Q_gpu.get()), 1e-4)
    True
    >>> # Numpy QR and CULA QR
    >>> A = np.array(np.random.randn(9, 7), np.float32, order='F')
    >>> Q, R = np.linalg.qr(A, 'reduced')
    >>> a_gpu = gpuarray.to_gpu(A)
    >>> Q_gpu, R_gpu = linalg.qr(a_gpu, 'reduced')
    >>> np.allclose(Q, Q_gpu.get(), 1e-4)
    True
    >>> np.allclose(R, R_gpu.get(), 1e-4)
    True
QR Decomposition.
QR Decomposition .
def qr(a_gpu, mode='reduced', handle=None, lib='cusolver'):
    """
    QR Decomposition.

    Factor the real/complex matrix `a` as `QR`, where `Q` is an orthonormal/unitary
    matrix and `R` is an upper triangular matrix.

    Parameters
    ----------
    a_gpu: pycuda.gpuarray.GPUArray
        Real/complex input matrix  `a` with dimensions `(m, n)`.
        `a` is assumed to be `m`>=`n`.
    mode :  {'reduced', 'economic', 'r'}
        'reduced' : returns `Q`, `R` with dimensions `(m, k)` and `(k, n)` (default).
        'economic' : returns `Q` only with dimensions `(m, k)`.
        'r' : returns `R` only with dimensions `(k, n)` with `k`=min`(m,n)`.
    handle : int
        CUBLAS context. If no context is specified, the default handle from
        `skcuda.misc._global_cublas_handle` is used.
    lib : str
        Library to use. May be either 'cula' or 'cusolver' (default).

    Returns
    -------
    q_gpu : pycuda.gpuarray.GPUArray
        Orthonormal/unitary matrix (depending on whether or not `A` is real/complex).
    r_gpu : pycuda.gpuarray.GPUArray
        The upper-triangular matrix.

    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 assumed to be stored in column-major order, i.e., order='F'.

    Examples
    --------
    >>> import pycuda.autoinit
    >>> import pycuda.gpuarray as gpuarray
    >>> import numpy as np
    >>> import skcuda.linalg as linalg
    >>> linalg.init()
    >>> # Rectangular matrix A, np.float32
    >>> A = np.array(np.random.randn(9, 7), np.float32, order='F')
    >>> A_gpu = gpuarray.to_gpu(A)
    >>> Q_gpu, R_gpu = linalg.qr(A_gpu, 'reduced')
    >>> np.allclose(A, np.dot(Q_gpu.get(), R_gpu.get()), 1e-4)
    True
    >>> # Square matrix A, np.complex128
    >>> A = np.random.randn(9, 9) + 1j*np.random.randn(9, 9)
    >>> A = np.asarray(A, np.complex128, order='F')
    >>> A_gpu = gpuarray.to_gpu(A)
    >>> Q_gpu, R_gpu = linalg.qr(A_gpu, 'reduced')
    >>> np.allclose(A, np.dot(Q_gpu.get(), R_gpu.get()), 1e-4)
    True
    >>> np.allclose(np.identity(Q_gpu.shape[0]) + 1j*0, np.dot(Q_gpu.get().conj().T, Q_gpu.get()), 1e-4)
    True
    >>> # Numpy QR and CULA QR
    >>> A = np.array(np.random.randn(9, 7), np.float32, order='F')
    >>> Q, R = np.linalg.qr(A, 'reduced')
    >>> a_gpu = gpuarray.to_gpu(A)
    >>> Q_gpu, R_gpu = linalg.qr(a_gpu, 'reduced')
    >>> np.allclose(Q, Q_gpu.get(), 1e-4)
    True
    >>> np.allclose(R, R_gpu.get(), 1e-4)
    True
    """

    alloc = misc._global_cublas_allocator

    if handle is None:
         handle = misc._global_cublas_handle

    data_type = a_gpu.dtype.type

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

        # The free version of CULA only supports single precision floating
        # point numbers:

        real_type = np.float32
        if data_type == np.complex64:
            func_qr = cula.culaDeviceCgeqrf
            func_q = cula.culaDeviceCungqr
            copy_func = cublas.cublasCcopy
            use_double = 0
            use_complex = 1
        elif data_type == np.float32:
            func_qr = cula.culaDeviceSgeqrf
            func_q = cula.culaDeviceSorgqr
            copy_func = cublas.cublasScopy
            use_double = 0
            use_complex = 0
        else:
            if cula._libcula_toolkit == 'standard':
                if data_type == np.complex128:
                    func_qr = cula.culaDeviceZgeqrf
                    func_q = cula.culaDeviceZungqr
                    copy_func = cublas.cublasZcopy
                    use_double = 1
                    use_complex = 1
                elif data_type == np.float64:
                    func_qr = cula.culaDeviceDgeqrf
                    func_q = cula.culaDeviceDorgqr
                    copy_func = cublas.cublasDcopy
                    use_double = 1
                    use_complex = 0
                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_qr = cusolver.cusolverDnCgeqrf
            func_q = cusolver.cusolverDnCungqr
            bufsize_qr = cusolver.cusolverDnCgeqrf_bufferSize
            bufsize_q = cusolver.cusolverDnCungqr_bufferSize
            copy_func = cublas.cublasCcopy
            use_double = 0
            use_complex = 1
        elif data_type == np.float32:
            func_qr = cusolver.cusolverDnSgeqrf
            func_q = cusolver.cusolverDnSorgqr
            bufsize_qr = cusolver.cusolverDnSgeqrf_bufferSize
            bufsize_q = cusolver.cusolverDnSorgqr_bufferSize
            copy_func = cublas.cublasScopy
            use_double = 0
            use_complex = 0
        elif data_type == np.complex128:
            real_type = np.float64
            func_qr = cusolver.cusolverDnZgeqrf
            func_q = cusolver.cusolverDnZungqr
            bufsize_qr = cusolver.cusolverDnZgeqrf_bufferSize
            bufsize_q = cusolver.cusolverDnZungqr_bufferSize
            copy_func = cublas.cublasZcopy
            use_double = 1
            use_complex = 1
        elif data_type == np.float64:
            real_type = np.float64
            func_qr = cusolver.cusolverDnDgeqrf
            func_q = cusolver.cusolverDnDorgqr
            bufsize_qr = cusolver.cusolverDnDgeqrf_bufferSize
            bufsize_q = cusolver.cusolverDnDorgqr_bufferSize
            copy_func = cublas.cublasDcopy
            use_double = 1
            use_complex = 0
        else:
            raise ValueError('unsupported type')

    else:
        raise ValueError('invalid library specified')

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

    if m<n and mode != 'r':
        raise ValueError('if m < n only the mode "r" is supported')

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

    # Set k:
    k = min(m, n)

    # Set the leading dimension and allocate u:
    tau_gpu = gpuarray.empty(k, data_type, allocator=alloc, order='F')

    # Compute QR and check error status:
    if lib == 'cula':
        func_qr(m, n, int(a_gpu.gpudata), lda, int(tau_gpu.gpudata))
    else:
        Lwork = bufsize_qr(cusolverHandle, m, n, int(a_gpu.gpudata), m)      
        workspace_gpu = gpuarray.empty(Lwork, data_type, allocator=alloc)
        devInfo = gpuarray.empty(1, np.int32, allocator=alloc)
        func_qr(cusolverHandle, m, n, int(a_gpu.gpudata), lda, int(tau_gpu.gpudata),
                int(workspace_gpu.gpudata), Lwork, int(devInfo.gpudata))

    if mode != 'economic':
        # Get upper triangular matrix R with dimensions (n,n)
        # Note: _get_tril_kernel returns the upper triangular
        r_gpu = gpuarray.empty((m, n), data_type, allocator=alloc, order='F')
        copy_func(handle, a_gpu.size, int(a_gpu.gpudata), 1, int(r_gpu.gpudata), 1)

        # tril
        dev = misc.get_current_device()
        block_dim, grid_dim = misc.select_block_grid_sizes(dev, r_gpu.shape)
        tril = _get_tril_kernel(use_double, use_complex, cols=m) #cols are here rows
        tril(r_gpu, np.uint32(r_gpu.size), block=block_dim, grid=grid_dim)

        # Mode r
        if mode == 'r':
            return r_gpu[:k, :n]

    # Compute Q and check error status:
    if lib == 'cula':
        func_q(m, n, k, int(a_gpu.gpudata), lda, int(tau_gpu.gpudata))

        # Free internal CULA memory:
        cula.culaFreeBuffers()
    else:
        Lwork = bufsize_q(cusolverHandle, m, n, k, int(a_gpu.gpudata), lda, int(tau_gpu.gpudata))
        workspace_gpu = gpuarray.empty(Lwork, data_type, allocator=alloc)
        # Reuse devInfo allocated earlier:
        func_q(cusolverHandle, m, n, k, int(a_gpu.gpudata), lda,
               int(tau_gpu.gpudata), int(workspace_gpu.gpudata), Lwork,
               int(devInfo.gpudata))
    q_gpu = a_gpu

    # Mode economic
    if mode == 'reduced':
        return q_gpu, r_gpu[:k, :n]
    if mode == 'economic':
        return q_gpu
def qr ( a_gpu , mode = 'reduced' , handle = None , lib = 'cusolver' ) : alloc = misc . _global_cublas_allocator if handle is None : handle = misc . _global_cublas_handle data_type = a_gpu . dtype . type if lib == 'cula' : if not _has_cula : raise NotImplementedError ( 'CULA not installed' ) # The free version of CULA only supports single precision floating # point numbers: real_type = np . float32 if data_type == np . complex64 : func_qr = cula . culaDeviceCgeqrf func_q = cula . culaDeviceCungqr copy_func = cublas . cublasCcopy use_double = 0 use_complex = 1 elif data_type == np . float32 : func_qr = cula . culaDeviceSgeqrf func_q = cula . culaDeviceSorgqr copy_func = cublas . cublasScopy use_double = 0 use_complex = 0 else : if cula . _libcula_toolkit == 'standard' : if data_type == np . complex128 : func_qr = cula . culaDeviceZgeqrf func_q = cula . culaDeviceZungqr copy_func = cublas . cublasZcopy use_double = 1 use_complex = 1 elif data_type == np . float64 : func_qr = cula . culaDeviceDgeqrf func_q = cula . culaDeviceDorgqr copy_func = cublas . cublasDcopy use_double = 1 use_complex = 0 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_qr = cusolver . cusolverDnCgeqrf func_q = cusolver . cusolverDnCungqr bufsize_qr = cusolver . cusolverDnCgeqrf_bufferSize bufsize_q = cusolver . cusolverDnCungqr_bufferSize copy_func = cublas . cublasCcopy use_double = 0 use_complex = 1 elif data_type == np . float32 : func_qr = cusolver . cusolverDnSgeqrf func_q = cusolver . cusolverDnSorgqr bufsize_qr = cusolver . cusolverDnSgeqrf_bufferSize bufsize_q = cusolver . cusolverDnSorgqr_bufferSize copy_func = cublas . cublasScopy use_double = 0 use_complex = 0 elif data_type == np . complex128 : real_type = np . float64 func_qr = cusolver . cusolverDnZgeqrf func_q = cusolver . cusolverDnZungqr bufsize_qr = cusolver . cusolverDnZgeqrf_bufferSize bufsize_q = cusolver . cusolverDnZungqr_bufferSize copy_func = cublas . cublasZcopy use_double = 1 use_complex = 1 elif data_type == np . float64 : real_type = np . float64 func_qr = cusolver . cusolverDnDgeqrf func_q = cusolver . cusolverDnDorgqr bufsize_qr = cusolver . cusolverDnDgeqrf_bufferSize bufsize_q = cusolver . cusolverDnDorgqr_bufferSize copy_func = cublas . cublasDcopy use_double = 1 use_complex = 0 else : raise ValueError ( 'unsupported type' ) else : raise ValueError ( 'invalid library specified' ) # CUDA assumes that arrays are stored in column-major order m , n = a_gpu . shape if m < n and mode != 'r' : raise ValueError ( 'if m < n only the mode "r" is supported' ) # Set the leading dimension of the input matrix: lda = max ( 1 , m ) # Set k: k = min ( m , n ) # Set the leading dimension and allocate u: tau_gpu = gpuarray . empty ( k , data_type , allocator = alloc , order = 'F' ) # Compute QR and check error status: if lib == 'cula' : func_qr ( m , n , int ( a_gpu . gpudata ) , lda , int ( tau_gpu . gpudata ) ) else : Lwork = bufsize_qr ( cusolverHandle , m , n , int ( a_gpu . gpudata ) , m ) workspace_gpu = gpuarray . empty ( Lwork , data_type , allocator = alloc ) devInfo = gpuarray . empty ( 1 , np . int32 , allocator = alloc ) func_qr ( cusolverHandle , m , n , int ( a_gpu . gpudata ) , lda , int ( tau_gpu . gpudata ) , int ( workspace_gpu . gpudata ) , Lwork , int ( devInfo . gpudata ) ) if mode != 'economic' : # Get upper triangular matrix R with dimensions (n,n) # Note: _get_tril_kernel returns the upper triangular r_gpu = gpuarray . empty ( ( m , n ) , data_type , allocator = alloc , order = 'F' ) copy_func ( handle , a_gpu . size , int ( a_gpu . gpudata ) , 1 , int ( r_gpu . gpudata ) , 1 ) # tril dev = misc . get_current_device ( ) block_dim , grid_dim = misc . select_block_grid_sizes ( dev , r_gpu . shape ) tril = _get_tril_kernel ( use_double , use_complex , cols = m ) #cols are here rows tril ( r_gpu , np . uint32 ( r_gpu . size ) , block = block_dim , grid = grid_dim ) # Mode r if mode == 'r' : return r_gpu [ : k , : n ] # Compute Q and check error status: if lib == 'cula' : func_q ( m , n , k , int ( a_gpu . gpudata ) , lda , int ( tau_gpu . gpudata ) ) # Free internal CULA memory: cula . culaFreeBuffers ( ) else : Lwork = bufsize_q ( cusolverHandle , m , n , k , int ( a_gpu . gpudata ) , lda , int ( tau_gpu . gpudata ) ) workspace_gpu = gpuarray . empty ( Lwork , data_type , allocator = alloc ) # Reuse devInfo allocated earlier: func_q ( cusolverHandle , m , n , k , int ( a_gpu . gpudata ) , lda , int ( tau_gpu . gpudata ) , int ( workspace_gpu . gpudata ) , Lwork , int ( devInfo . gpudata ) ) q_gpu = a_gpu # Mode economic if mode == 'reduced' : return q_gpu , r_gpu [ : k , : n ] if mode == 'economic' : return q_gpu