code: 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
  
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 1 for calling a global function with example usage: svd(a_gpu, jobu, jobvt, lib) and returns: vh_gpu,  s_gpu,  u_gpu, vh_gpu,  s_gpu, s_gpu,  u_gpu, s_gpu