(a_gpu, k=None, modes='exact', return_amplitudes=False, return_vandermonde=False, handle=None)
Dynamic Mode Decomposition.

    Dynamic Mode Decomposition (DMD) is a data processing algorithm which
    allows to decompose a matrix `a` in space and time.
    The matrix `a` is decomposed as `a = FBV`, where the columns of `F`
    contain the dynamic modes. The modes are ordered corresponding
    to the amplitudes stored in the diagonal matrix `B`. `V` is a Vandermonde
    matrix describing the temporal evolution.

    Parameters
    ----------
    a_gpu : pycuda.gpuarray.GPUArray
        Real/complex input matrix  `a` with dimensions `(m, n)`.
    k : int, optional
        If `k < (n-1)` low-rank Dynamic Mode Decomposition is computed.
    modes : `{'standard', 'exact'}`
        'standard' : uses the standard definition to compute the dynamic modes,
                    `F = U * W`.
        'exact' : computes the exact dynamic modes, `F = Y * V * (S**-1) * W`.
    return_amplitudes : bool `{True, False}`
        True: return amplitudes in addition to dynamic modes.
    return_vandermonde : bool `{True, False}`
        True: return Vandermonde matrix in addition to dynamic modes and amplitudes.
    handle : int
        CUBLAS context. If no context is specified, the default handle from
        `skcuda.misc._global_cublas_handle` is used.

    Returns
    -------
    f_gpu : pycuda.gpuarray.GPUArray
        Matrix containing the dynamic modes of shape `(m, n-1)`  or `(m, k)`.
    b_gpu : pycuda.gpuarray.GPUArray
        1-D array containing the amplitudes of length `min(n-1, k)`.
    v_gpu : pycuda.gpuarray.GPUArray
        Vandermonde matrix of shape `(n-1, n-1)`  or `(k, n-1)`.

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

    References
    ----------
    M. R. Jovanovic, P. J. Schmid, and J. W. Nichols.
    "Low-rank and sparse dynamic mode decomposition."
    Center for Turbulence Research Annual Research Briefs (2012): 139-152.

    J. H. Tu, et al.
    "On dynamic mode decomposition: theory and applications."
    arXiv preprint arXiv:1312.0041 (2013).


     Examples
     --------
     >>> #Numpy
     >>> import numpy as np
     >>> #Plot libs
     >>> import matplotlib.pyplot as plt
     >>> from mpl_toolkits.mplot3d import Axes3D
     >>> from matplotlib import cm
     >>> #GPU DMD libs
     >>> import pycuda.gpuarray as gpuarray
     >>> import pycuda.autoinit
     >>> from skcuda import linalg, rlinalg
     >>> linalg.init()

     >>> # Define time and space discretizations
     >>> x=np.linspace( -15, 15, 200)
     >>> t=np.linspace(0, 8*np.pi , 80)
     >>> dt=t[2]-t[1]
     >>> X, T = np.meshgrid(x,t)
     >>> # Create two patio-temporal patterns
     >>> F1 = 0.5* np.cos(X)*(1.+0.* T)
     >>> F2 = ( (1./np.cosh(X)) * np.tanh(X)) *(2.*np.exp(1j*2.8*T))
     >>> # Add both signals
     >>> F = (F1+F2)

     >>> #Plot dataset
     >>> fig = plt.figure()
     >>> ax = fig.add_subplot(231, projection='3d')
     >>> ax = fig.gca(projection='3d')
     >>> surf = ax.plot_surface(X, T, F, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
     >>> ax.set_zlim(-1, 1)
     >>> plt.title('F')
     >>> ax = fig.add_subplot(232, projection='3d')
     >>> ax = fig.gca(projection='3d')
     >>> surf = ax.plot_surface(X, T, F1, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
     >>> ax.set_zlim(-1, 1)
     >>> plt.title('F1')
     >>> ax = fig.add_subplot(233, projection='3d')
     >>> ax = fig.gca(projection='3d')
     >>> surf = ax.plot_surface(X, T, F2, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
     >>> ax.set_zlim(-1, 1)
     >>> plt.title('F2')

     >>> #Dynamic Mode Decomposition
     >>> F_gpu = np.array(F.T, np.complex64, order='F')
     >>> F_gpu = gpuarray.to_gpu(F_gpu)
     >>> Fmodes_gpu, b_gpu, V_gpu, omega_gpu = linalg.dmd(F_gpu, k=2, modes='exact', return_amplitudes=True, return_vandermonde=True)
     >>> omega = omega_gpu.get()
     >>> plt.scatter(omega.real, omega.imag, marker='o', c='r')


     >>> #Recover original signal
     >>> F1tilde = np.dot(Fmodes_gpu[:,0:1].get() , np.dot(b_gpu[0].get(), V_gpu[0:1,:].get() ) )
     >>> F2tilde = np.dot(Fmodes_gpu[:,1:2].get() , np.dot(b_gpu[1].get(), V_gpu[1:2,:].get() ) )

     >>> #Plot DMD modes
     >>> #Mode 0
     >>> ax = fig.add_subplot(235, projection='3d')
     >>> ax = fig.gca(projection='3d')
     >>> surf = ax.plot_surface(X[0:F1tilde.shape[1],:], T[0:F1tilde.shape[1],:], F1tilde.T, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
     >>> ax.set_zlim(-1, 1)
     >>> plt.title('F1_tilde')
     >>> #Mode 1
     >>> ax = fig.add_subplot(236, projection='3d')
     >>> ax = fig.gca(projection='3d')
     >>> surf = ax.plot_surface(X[0:F2tilde.shape[1],:], T[0:F2tilde.shape[1],:], F2tilde.T, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
     >>> ax.set_zlim(-1, 1)
     >>> plt.title('F2_tilde')
     >>> plt.show()
Dynamic Mode Decomposition.
Dynamic Mode Decomposition .
def dmd(a_gpu, k=None, modes='exact', return_amplitudes=False, return_vandermonde=False, handle=None):
    """
    Dynamic Mode Decomposition.

    Dynamic Mode Decomposition (DMD) is a data processing algorithm which
    allows to decompose a matrix `a` in space and time.
    The matrix `a` is decomposed as `a = FBV`, where the columns of `F`
    contain the dynamic modes. The modes are ordered corresponding
    to the amplitudes stored in the diagonal matrix `B`. `V` is a Vandermonde
    matrix describing the temporal evolution.

    Parameters
    ----------
    a_gpu : pycuda.gpuarray.GPUArray
        Real/complex input matrix  `a` with dimensions `(m, n)`.
    k : int, optional
        If `k < (n-1)` low-rank Dynamic Mode Decomposition is computed.
    modes : `{'standard', 'exact'}`
        'standard' : uses the standard definition to compute the dynamic modes,
                    `F = U * W`.
        'exact' : computes the exact dynamic modes, `F = Y * V * (S**-1) * W`.
    return_amplitudes : bool `{True, False}`
        True: return amplitudes in addition to dynamic modes.
    return_vandermonde : bool `{True, False}`
        True: return Vandermonde matrix in addition to dynamic modes and amplitudes.
    handle : int
        CUBLAS context. If no context is specified, the default handle from
        `skcuda.misc._global_cublas_handle` is used.

    Returns
    -------
    f_gpu : pycuda.gpuarray.GPUArray
        Matrix containing the dynamic modes of shape `(m, n-1)`  or `(m, k)`.
    b_gpu : pycuda.gpuarray.GPUArray
        1-D array containing the amplitudes of length `min(n-1, k)`.
    v_gpu : pycuda.gpuarray.GPUArray
        Vandermonde matrix of shape `(n-1, n-1)`  or `(k, n-1)`.

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

    References
    ----------
    M. R. Jovanovic, P. J. Schmid, and J. W. Nichols.
    "Low-rank and sparse dynamic mode decomposition."
    Center for Turbulence Research Annual Research Briefs (2012): 139-152.

    J. H. Tu, et al.
    "On dynamic mode decomposition: theory and applications."
    arXiv preprint arXiv:1312.0041 (2013).


     Examples
     --------
     >>> #Numpy
     >>> import numpy as np
     >>> #Plot libs
     >>> import matplotlib.pyplot as plt
     >>> from mpl_toolkits.mplot3d import Axes3D
     >>> from matplotlib import cm
     >>> #GPU DMD libs
     >>> import pycuda.gpuarray as gpuarray
     >>> import pycuda.autoinit
     >>> from skcuda import linalg, rlinalg
     >>> linalg.init()

     >>> # Define time and space discretizations
     >>> x=np.linspace( -15, 15, 200)
     >>> t=np.linspace(0, 8*np.pi , 80)
     >>> dt=t[2]-t[1]
     >>> X, T = np.meshgrid(x,t)
     >>> # Create two patio-temporal patterns
     >>> F1 = 0.5* np.cos(X)*(1.+0.* T)
     >>> F2 = ( (1./np.cosh(X)) * np.tanh(X)) *(2.*np.exp(1j*2.8*T))
     >>> # Add both signals
     >>> F = (F1+F2)

     >>> #Plot dataset
     >>> fig = plt.figure()
     >>> ax = fig.add_subplot(231, projection='3d')
     >>> ax = fig.gca(projection='3d')
     >>> surf = ax.plot_surface(X, T, F, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
     >>> ax.set_zlim(-1, 1)
     >>> plt.title('F')
     >>> ax = fig.add_subplot(232, projection='3d')
     >>> ax = fig.gca(projection='3d')
     >>> surf = ax.plot_surface(X, T, F1, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
     >>> ax.set_zlim(-1, 1)
     >>> plt.title('F1')
     >>> ax = fig.add_subplot(233, projection='3d')
     >>> ax = fig.gca(projection='3d')
     >>> surf = ax.plot_surface(X, T, F2, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
     >>> ax.set_zlim(-1, 1)
     >>> plt.title('F2')

     >>> #Dynamic Mode Decomposition
     >>> F_gpu = np.array(F.T, np.complex64, order='F')
     >>> F_gpu = gpuarray.to_gpu(F_gpu)
     >>> Fmodes_gpu, b_gpu, V_gpu, omega_gpu = linalg.dmd(F_gpu, k=2, modes='exact', return_amplitudes=True, return_vandermonde=True)
     >>> omega = omega_gpu.get()
     >>> plt.scatter(omega.real, omega.imag, marker='o', c='r')


     >>> #Recover original signal
     >>> F1tilde = np.dot(Fmodes_gpu[:,0:1].get() , np.dot(b_gpu[0].get(), V_gpu[0:1,:].get() ) )
     >>> F2tilde = np.dot(Fmodes_gpu[:,1:2].get() , np.dot(b_gpu[1].get(), V_gpu[1:2,:].get() ) )

     >>> #Plot DMD modes
     >>> #Mode 0
     >>> ax = fig.add_subplot(235, projection='3d')
     >>> ax = fig.gca(projection='3d')
     >>> surf = ax.plot_surface(X[0:F1tilde.shape[1],:], T[0:F1tilde.shape[1],:], F1tilde.T, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
     >>> ax.set_zlim(-1, 1)
     >>> plt.title('F1_tilde')
     >>> #Mode 1
     >>> ax = fig.add_subplot(236, projection='3d')
     >>> ax = fig.gca(projection='3d')
     >>> surf = ax.plot_surface(X[0:F2tilde.shape[1],:], T[0:F2tilde.shape[1],:], F2tilde.T, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
     >>> ax.set_zlim(-1, 1)
     >>> plt.title('F2_tilde')
     >>> plt.show()
    """

    #*************************************************************************
    #***        Author: N. Benjamin Erichson <nbe@st-andrews.ac.uk>        ***
    #***                              <2015>                               ***
    #***                       License: BSD 3 clause                       ***
    #*************************************************************************

    if not _has_cula:
        raise NotImplementedError('CULA not installed')

    if handle is None:
        handle = misc._global_cublas_handle

    alloc = misc._global_cublas_allocator

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

    if data_type == np.complex64:
        cula_func_gesvd = cula.culaDeviceCgesvd
        cublas_func_gemm = cublas.cublasCgemm
        cublas_func_dgmm = cublas.cublasCdgmm
        cula_func_gels = cula.culaDeviceCgels
        copy_func = cublas.cublasCcopy
        transpose_func = cublas.cublasCgeam
        alpha = np.complex64(1.0)
        beta = np.complex64(0.0)
        TRANS_type = 'C'
        isreal = False
    elif data_type == np.float32:
        cula_func_gesvd = cula.culaDeviceSgesvd
        cublas_func_gemm = cublas.cublasSgemm
        cublas_func_dgmm = cublas.cublasSdgmm
        cula_func_gels = cula.culaDeviceSgels
        copy_func = cublas.cublasScopy
        transpose_func = cublas.cublasSgeam
        alpha = np.float32(1.0)
        beta = np.float32(0.0)
        TRANS_type = 'T'
        isreal = True
    else:
        if cula._libcula_toolkit == 'standard':
            if data_type == np.complex128:
                cula_func_gesvd = cula.culaDeviceZgesvd
                cublas_func_gemm = cublas.cublasZgemm
                cublas_func_dgmm = cublas.cublasZdgmm
                cula_func_gels = cula.culaDeviceZgels
                copy_func = cublas.cublasZcopy
                transpose_func = cublas.cublasZgeam
                alpha = np.complex128(1.0)
                beta = np.complex128(0.0)
                TRANS_type = 'C'
                isreal = False
            elif data_type == np.float64:
                cula_func_gesvd = cula.culaDeviceDgesvd
                cublas_func_gemm = cublas.cublasDgemm
                cublas_func_dgmm = cublas.cublasDdgmm
                cula_func_gels = cula.culaDeviceDgels
                copy_func = cublas.cublasDcopy
                transpose_func = cublas.cublasDgeam
                alpha = np.float64(1.0)
                beta = np.float64(0.0)
                TRANS_type = 'T'
                isreal = True
            else:
                raise ValueError('unsupported type')
            real_type = np.float64
        else:
            raise ValueError('double precision not supported')

    #CUDA assumes that arrays are stored in column-major order
    m, n = a_gpu.shape
    nx = n-1
    #Set k
    if k == None : k = nx
    if k > nx or k < 1: raise ValueError('k is not valid')

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Split data into lef and right snapshot sequence
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Note: we need a copy of X_gpu, because SVD destroys X_gpu
    #While Y_gpu is just a pointer
    X_gpu = gpuarray.empty((m, n), data_type, order="F", allocator=alloc)
    copy_func(handle, X_gpu.size, int(a_gpu.gpudata), 1, int(X_gpu.gpudata), 1)
    X_gpu = X_gpu[:, :nx]
    Y_gpu = a_gpu[:, 1:]

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Singular Value Decomposition
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #gesvd(jobu, jobvt, m, n, int(a), lda, int(s), int(u), ldu, int(vt), ldvt)
    #Parameters
    #----------
    #a : pycuda.gpuarray.GPUArray of shape (m, n)
    #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, nx)`.
    #    If 'O', return the `u` matrix with shape `(m, nx) without
    #    allocating a new matrix.
    #jobvt : {'A', 'S', 'O', 'N'}
    #    If 'A', return the full `vh` matrix with shape `(nx, nx)`.
    #    If 'S', return the `vh` matrix with shape `(nx, nx)`.
    #    If 'O', return the `vh` matrix with shape `(nx, nx) without
    #    allocating a new matrix.
    #
    #Returns
    #-------
    #u : pycuda.gpuarray.GPUArray
    #    Unitary matrix of shape `(m, m)` or `(m, nx)`
    #s : pycuda.gpuarray.GPUArray
    #    Array containing the singular values, sorted such that `s[i] >= s[i+1]`.
    #    `s` is of length `min(m, nx)`.
    #v : pycuda.gpuarray.GPUArray
    #    Unitary matrix of shape `(nx, nx)` or `(nx, nx)`

    #Allocate s, U, Vt for economic SVD
    #Note: singular values are always real
    #Allocate s, U, Vt for economic SVD
    #Note: singular values are always real
    s_gpu = gpuarray.empty(nx, real_type, order="F", allocator=alloc)
    U_gpu = gpuarray.empty((m,nx), data_type, order="F", allocator=alloc)
    Vh_gpu = gpuarray.empty((nx,nx), data_type, order="F", allocator=alloc)

    #Economic SVD
    cula_func_gesvd('S', 'S', m, nx, int(X_gpu.gpudata), m, int(s_gpu.gpudata),
                    int(U_gpu.gpudata), m, int(Vh_gpu.gpudata), nx)

    #Low-rank DMD: trancate SVD if k < nx

    if k != nx:
        s_gpu = s_gpu[:k]
        U_gpu = U_gpu[: , :k]
        #Vt_gpu = Vt_gpu[:k , : ]
        Vh_gpu = Vh_gpu[:k , : ]


    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Solve the LS problem to find estimate for M using the pseudo-inverse
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #real: M = U.T * Y * Vt.T * S**-1
    #complex: M = U.H * Y * Vt.H * S**-1
    #Let G = Y * Vt.H * S**-1, hence M = M * G

    #Allocate G and M
    G_gpu = gpuarray.empty((m,k), data_type, order="F", allocator=alloc)
    M_gpu = gpuarray.empty((k,k), data_type, order="F", allocator=alloc)

    #i) s = s **-1 (inverse)
    if data_type == np.complex64 or data_type == np.complex128:
        s_gpu = 1/s_gpu
        s_gpu = s_gpu + 1j * gpuarray.zeros_like(s_gpu)
    else:
        s_gpu = 1.0/s_gpu


    #ii) real/complex: scale Vs =  Vt* x diag(s**-1)
    Vs_gpu = gpuarray.empty((nx,k), data_type, order="F", allocator=alloc)
    lda = max(1, Vh_gpu.strides[1] // Vh_gpu.dtype.itemsize)
    ldb = max(1, Vs_gpu.strides[1] // Vs_gpu.dtype.itemsize)
    transpose_func(handle, TRANS_type, TRANS_type, nx, k,
                   alpha, int(Vh_gpu.gpudata), lda, beta, int(Vh_gpu.gpudata), lda,
                   int(Vs_gpu.gpudata), ldb)


    cublas_func_dgmm(handle, 'r', nx, k, int(Vs_gpu.gpudata), nx,
                     int(s_gpu.gpudata), 1 , int(Vs_gpu.gpudata), nx)


    #iii) real: G = Y * Vs , complex: G = Y x Vs
    cublas_func_gemm(handle, 'n', 'n', m, k, nx, alpha,
                     int(Y_gpu.gpudata), m, int(Vs_gpu.gpudata), nx,
                        beta, int(G_gpu.gpudata), m )


    #iv) real/complex: M = U* x G
    cublas_func_gemm(handle, TRANS_type, 'n', k, k, m, alpha,
                     int(U_gpu.gpudata), m, int(G_gpu.gpudata), m,
                    beta, int(M_gpu.gpudata), k )

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Eigen Decomposition
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Note: If a_gpu is real the imag part is omitted
    Vr_gpu, w_gpu = eig(M_gpu, 'N', 'V', 'F')
    omega = cumath.log(w_gpu)

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Compute DMD Modes
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    F_gpu = gpuarray.empty((m,k), data_type, order="F", allocator=alloc)
    modes = modes.lower()
    if modes == 'exact': #Compute (exact) DMD modes: F = Y * V * S**-1 * W = G * W
        cublas_func_gemm(handle, 'n', 'n', m, k, k, alpha,
                         G_gpu.gpudata, m, Vr_gpu.gpudata, k,
                         beta, G_gpu.gpudata, m  )
        F_gpu_temp = G_gpu

    elif modes == 'standard': #Compute (standard) DMD modes: F = U * W
        cublas_func_gemm(handle, 'n', 'n', m, k, k,
                         alpha, U_gpu.gpudata, m, Vr_gpu.gpudata, k,
                         beta, U_gpu.gpudata, m  )
        F_gpu_temp = U_gpu
    else:
        raise ValueError('Type of modes is not supported, choose "exact" or "standard".')

    #Copy is required, because gels destroys input
    copy_func(handle, F_gpu_temp.size, int(F_gpu_temp.gpudata),
              1, int(F_gpu.gpudata), 1)

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Compute amplitueds b using least-squares: Fb=x1
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if return_amplitudes==True:
        #x1_gpu = a_gpu[:,0].copy()
        x1_gpu = gpuarray.empty(m, data_type, order="F", allocator=alloc)
        copy_func(handle, x1_gpu.size, int(a_gpu[:,0].gpudata), 1, int(x1_gpu.gpudata), 1)
        cula_func_gels( 'N', m, k, int(1) , F_gpu_temp.gpudata, m, x1_gpu.gpudata, m)
        b_gpu = x1_gpu

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Compute Vandermonde matrix (CPU)
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if return_vandermonde==True:
        V_gpu = vander(w_gpu, n=nx)

    # Free internal CULA memory:
    cula.culaFreeBuffers()

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Return
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if return_amplitudes==True and return_vandermonde==True:
        return F_gpu, b_gpu[:k], V_gpu, omega
    elif return_amplitudes==True and return_vandermonde==False:
        return F_gpu, b_gpu[:k], omega
    elif return_amplitudes==False and return_vandermonde==True:
        return F_gpu, V_gpu, omega
    else:
        return F_gpu, omega
def dmd ( a_gpu , k = None , modes = 'exact' , return_amplitudes = False , return_vandermonde = False , handle = None ) : #************************************************************************* #***        Author: N. Benjamin Erichson <nbe@st-andrews.ac.uk>        *** #***                              <2015>                               *** #***                       License: BSD 3 clause                       *** #************************************************************************* if not _has_cula : raise NotImplementedError ( 'CULA not installed' ) if handle is None : handle = misc . _global_cublas_handle alloc = misc . _global_cublas_allocator # The free version of CULA only supports single precision floating data_type = a_gpu . dtype . type real_type = np . float32 if data_type == np . complex64 : cula_func_gesvd = cula . culaDeviceCgesvd cublas_func_gemm = cublas . cublasCgemm cublas_func_dgmm = cublas . cublasCdgmm cula_func_gels = cula . culaDeviceCgels copy_func = cublas . cublasCcopy transpose_func = cublas . cublasCgeam alpha = np . complex64 ( 1.0 ) beta = np . complex64 ( 0.0 ) TRANS_type = 'C' isreal = False elif data_type == np . float32 : cula_func_gesvd = cula . culaDeviceSgesvd cublas_func_gemm = cublas . cublasSgemm cublas_func_dgmm = cublas . cublasSdgmm cula_func_gels = cula . culaDeviceSgels copy_func = cublas . cublasScopy transpose_func = cublas . cublasSgeam alpha = np . float32 ( 1.0 ) beta = np . float32 ( 0.0 ) TRANS_type = 'T' isreal = True else : if cula . _libcula_toolkit == 'standard' : if data_type == np . complex128 : cula_func_gesvd = cula . culaDeviceZgesvd cublas_func_gemm = cublas . cublasZgemm cublas_func_dgmm = cublas . cublasZdgmm cula_func_gels = cula . culaDeviceZgels copy_func = cublas . cublasZcopy transpose_func = cublas . cublasZgeam alpha = np . complex128 ( 1.0 ) beta = np . complex128 ( 0.0 ) TRANS_type = 'C' isreal = False elif data_type == np . float64 : cula_func_gesvd = cula . culaDeviceDgesvd cublas_func_gemm = cublas . cublasDgemm cublas_func_dgmm = cublas . cublasDdgmm cula_func_gels = cula . culaDeviceDgels copy_func = cublas . cublasDcopy transpose_func = cublas . cublasDgeam alpha = np . float64 ( 1.0 ) beta = np . float64 ( 0.0 ) TRANS_type = 'T' isreal = True else : raise ValueError ( 'unsupported type' ) real_type = np . float64 else : raise ValueError ( 'double precision not supported' ) #CUDA assumes that arrays are stored in column-major order m , n = a_gpu . shape nx = n - 1 #Set k if k == None : k = nx if k > nx or k < 1 : raise ValueError ( 'k is not valid' ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Split data into lef and right snapshot sequence #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Note: we need a copy of X_gpu, because SVD destroys X_gpu #While Y_gpu is just a pointer X_gpu = gpuarray . empty ( ( m , n ) , data_type , order = "F" , allocator = alloc ) copy_func ( handle , X_gpu . size , int ( a_gpu . gpudata ) , 1 , int ( X_gpu . gpudata ) , 1 ) X_gpu = X_gpu [ : , : nx ] Y_gpu = a_gpu [ : , 1 : ] #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Singular Value Decomposition #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #gesvd(jobu, jobvt, m, n, int(a), lda, int(s), int(u), ldu, int(vt), ldvt) #Parameters #---------- #a : pycuda.gpuarray.GPUArray of shape (m, n) #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, nx)`. #    If 'O', return the `u` matrix with shape `(m, nx) without #    allocating a new matrix. #jobvt : {'A', 'S', 'O', 'N'} #    If 'A', return the full `vh` matrix with shape `(nx, nx)`. #    If 'S', return the `vh` matrix with shape `(nx, nx)`. #    If 'O', return the `vh` matrix with shape `(nx, nx) without #    allocating a new matrix. # #Returns #------- #u : pycuda.gpuarray.GPUArray #    Unitary matrix of shape `(m, m)` or `(m, nx)` #s : pycuda.gpuarray.GPUArray #    Array containing the singular values, sorted such that `s[i] >= s[i+1]`. #    `s` is of length `min(m, nx)`. #v : pycuda.gpuarray.GPUArray #    Unitary matrix of shape `(nx, nx)` or `(nx, nx)` #Allocate s, U, Vt for economic SVD #Note: singular values are always real #Allocate s, U, Vt for economic SVD #Note: singular values are always real s_gpu = gpuarray . empty ( nx , real_type , order = "F" , allocator = alloc ) U_gpu = gpuarray . empty ( ( m , nx ) , data_type , order = "F" , allocator = alloc ) Vh_gpu = gpuarray . empty ( ( nx , nx ) , data_type , order = "F" , allocator = alloc ) #Economic SVD cula_func_gesvd ( 'S' , 'S' , m , nx , int ( X_gpu . gpudata ) , m , int ( s_gpu . gpudata ) , int ( U_gpu . gpudata ) , m , int ( Vh_gpu . gpudata ) , nx ) #Low-rank DMD: trancate SVD if k < nx if k != nx : s_gpu = s_gpu [ : k ] U_gpu = U_gpu [ : , : k ] #Vt_gpu = Vt_gpu[:k , : ] Vh_gpu = Vh_gpu [ : k , : ] #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Solve the LS problem to find estimate for M using the pseudo-inverse #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #real: M = U.T * Y * Vt.T * S**-1 #complex: M = U.H * Y * Vt.H * S**-1 #Let G = Y * Vt.H * S**-1, hence M = M * G #Allocate G and M G_gpu = gpuarray . empty ( ( m , k ) , data_type , order = "F" , allocator = alloc ) M_gpu = gpuarray . empty ( ( k , k ) , data_type , order = "F" , allocator = alloc ) #i) s = s **-1 (inverse) if data_type == np . complex64 or data_type == np . complex128 : s_gpu = 1 / s_gpu s_gpu = s_gpu + 1j * gpuarray . zeros_like ( s_gpu ) else : s_gpu = 1.0 / s_gpu #ii) real/complex: scale Vs =  Vt* x diag(s**-1) Vs_gpu = gpuarray . empty ( ( nx , k ) , data_type , order = "F" , allocator = alloc ) lda = max ( 1 , Vh_gpu . strides [ 1 ] // Vh_gpu . dtype . itemsize ) ldb = max ( 1 , Vs_gpu . strides [ 1 ] // Vs_gpu . dtype . itemsize ) transpose_func ( handle , TRANS_type , TRANS_type , nx , k , alpha , int ( Vh_gpu . gpudata ) , lda , beta , int ( Vh_gpu . gpudata ) , lda , int ( Vs_gpu . gpudata ) , ldb ) cublas_func_dgmm ( handle , 'r' , nx , k , int ( Vs_gpu . gpudata ) , nx , int ( s_gpu . gpudata ) , 1 , int ( Vs_gpu . gpudata ) , nx ) #iii) real: G = Y * Vs , complex: G = Y x Vs cublas_func_gemm ( handle , 'n' , 'n' , m , k , nx , alpha , int ( Y_gpu . gpudata ) , m , int ( Vs_gpu . gpudata ) , nx , beta , int ( G_gpu . gpudata ) , m ) #iv) real/complex: M = U* x G cublas_func_gemm ( handle , TRANS_type , 'n' , k , k , m , alpha , int ( U_gpu . gpudata ) , m , int ( G_gpu . gpudata ) , m , beta , int ( M_gpu . gpudata ) , k ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Eigen Decomposition #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Note: If a_gpu is real the imag part is omitted Vr_gpu , w_gpu = eig ( M_gpu , 'N' , 'V' , 'F' ) omega = cumath . log ( w_gpu ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Compute DMD Modes #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ F_gpu = gpuarray . empty ( ( m , k ) , data_type , order = "F" , allocator = alloc ) modes = modes . lower ( ) if modes == 'exact' : #Compute (exact) DMD modes: F = Y * V * S**-1 * W = G * W cublas_func_gemm ( handle , 'n' , 'n' , m , k , k , alpha , G_gpu . gpudata , m , Vr_gpu . gpudata , k , beta , G_gpu . gpudata , m ) F_gpu_temp = G_gpu elif modes == 'standard' : #Compute (standard) DMD modes: F = U * W cublas_func_gemm ( handle , 'n' , 'n' , m , k , k , alpha , U_gpu . gpudata , m , Vr_gpu . gpudata , k , beta , U_gpu . gpudata , m ) F_gpu_temp = U_gpu else : raise ValueError ( 'Type of modes is not supported, choose "exact" or "standard".' ) #Copy is required, because gels destroys input copy_func ( handle , F_gpu_temp . size , int ( F_gpu_temp . gpudata ) , 1 , int ( F_gpu . gpudata ) , 1 ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Compute amplitueds b using least-squares: Fb=x1 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if return_amplitudes == True : #x1_gpu = a_gpu[:,0].copy() x1_gpu = gpuarray . empty ( m , data_type , order = "F" , allocator = alloc ) copy_func ( handle , x1_gpu . size , int ( a_gpu [ : , 0 ] . gpudata ) , 1 , int ( x1_gpu . gpudata ) , 1 ) cula_func_gels ( 'N' , m , k , int ( 1 ) , F_gpu_temp . gpudata , m , x1_gpu . gpudata , m ) b_gpu = x1_gpu #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Compute Vandermonde matrix (CPU) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if return_vandermonde == True : V_gpu = vander ( w_gpu , n = nx ) # Free internal CULA memory: cula . culaFreeBuffers ( ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Return #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if return_amplitudes == True and return_vandermonde == True : return F_gpu , b_gpu [ : k ] , V_gpu , omega elif return_amplitudes == True and return_vandermonde == False : return F_gpu , b_gpu [ : k ] , omega elif return_amplitudes == False and return_vandermonde == True : return F_gpu , V_gpu , omega else : return F_gpu , omega