code: 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
  
imports: from __future__ import absolute_import, division
  from pprint import pprint
  from string import Template
  from pycuda.tools import context_dependent_memoize
  from pycuda.compiler import SourceModule
  from pycuda.reduction import ReductionKernel
  from pycuda import cumath
  import pycuda.gpuarray
  import pycuda.driver
  import pycuda.elementwise
  import pycuda.tools
  import numpy
  from .None import cublas
  from .None import cudart
  from .None import misc
  from .None import cusolver
  import sys
  from .misc import init, shutdown, add_matvec, div_matvec, mult_matvec
  from .None import install_headers
label_desc: Use image node_id 35 for calling a global function with example usage: dmd(a_gpu, k, modes, return_amplitudes, return_vandermonde, handle) and returns: F_gpu,  b_gpu,  V_gpu,  omega, F_gpu,  b_gpu,  omega, F_gpu,  V_gpu,  omega, F_gpu,  omega