#!/usr/bin/env python
"""
PyCUDA-based randomized linear algebra functions.
"""
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 curandom
from pycuda import cumath
import pycuda.gpuarray as gpuarray
import pycuda.driver as drv
import pycuda.elementwise as el
import pycuda.tools as tools
import numpy as np
from . import cublas
from . import misc
from . import linalg
rand = curandom.MRG32k3aRandomNumberGenerator()
import sys
if sys.version_info < (3,):
range = xrange
class LinAlgError(Exception):
"""Randomized Linear Algebra Error."""
pass
try:
from . import cula
_has_cula = True
except (ImportError, OSError):
_has_cula = False
from .misc import init, add_matvec, div_matvec, mult_matvec
from .linalg import hermitian, transpose
# Get installation location of C headers:
from . import install_headers
[docs]def rsvd(a_gpu, k=None, p=0, q=0, method="standard", handle=None):
"""
Randomized Singular Value Decomposition.
Randomized algorithm for computing the approximate low-rank singular value
decomposition of a rectangular (m, n) matrix `a` with target rank `k << n`.
The input matrix a is factored as `a = U * diag(s) * Vt`. The right singluar
vectors are the columns of the real or complex unitary matrix `U`. The left
singular vectors are the columns of the real or complex unitary matrix `V`.
The singular values `s` are non-negative and real numbers.
The paramter `p` is a oversampling parameter to improve the approximation.
A value between 2 and 10 is recommended.
The paramter `q` specifies the number of normlized power iterations
(subspace iterations) to reduce the approximation error. This is recommended
if the the singular values decay slowly and in practice 1 or 2 iterations
achive good results. However, computing power iterations is increasing the
computational time.
If k > (n/1.5), partial SVD or trancated SVD might be faster.
Parameters
----------
a_gpu : pycuda.gpuarray.GPUArray
Real/complex input matrix `a` with dimensions `(m, n)`.
k : int
`k` is the target rank of the low-rank decomposition, k << min(m,n).
p : int
`p` sets the oversampling parameter (default k=0).
q : int
`q` sets the number of power iterations (default=0).
method : `{'standard', 'fast'}`
'standard' : Standard algorithm as described in [1, 2]
'fast' : Version II algorithm as described in [2]
handle : int
CUBLAS context. If no context is specified, the default handle from
`skcuda.misc._global_cublas_handle` is used.
Returns
-------
u_gpu : pycuda.gpuarray
Right singular values, array of shape `(m, k)`.
s_gpu : pycuda.gpuarray
Singular values, 1-d array of length `k`.
vt_gpu : pycuda.gpuarray
Left singular values, array of shape `(k, n)`.
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'.
Input matrix of shape `(m, n)`, where `n>m` is not supported yet.
References
----------
N. Halko, P. Martinsson, and J. Tropp.
"Finding structure with randomness: probabilistic
algorithms for constructing approximate matrix
decompositions" (2009).
(available at `arXiv <http://arxiv.org/abs/0909.4061>`_).
S. Voronin and P.Martinsson.
"RSVDPACK: Subroutines for computing partial singular value
decompositions via randomized sampling on single core, multi core,
and GPU architectures" (2015).
(available at `arXiv <http://arxiv.org/abs/1502.05366>`_).
Examples
--------
>>> import pycuda.gpuarray as gpuarray
>>> import pycuda.autoinit
>>> import numpy as np
>>> from skcuda import linalg, rlinalg
>>> linalg.init()
>>> rlinalg.init()
>>> #Randomized SVD decomposition of the square matrix `a` with single precision.
>>> #Note: There is no gain to use rsvd if k > int(n/1.5)
>>> a = np.array(np.random.randn(5, 5), np.float32, order='F')
>>> a_gpu = gpuarray.to_gpu(a)
>>> U, s, Vt = rlinalg.rsvd(a_gpu, k=5, method='standard')
>>> np.allclose(a, np.dot(U.get(), np.dot(np.diag(s.get()), Vt.get())), 1e-4)
True
>>> #Low-rank SVD decomposition with target rank k=2
>>> a = np.array(np.random.randn(5, 5), np.float32, order='F')
>>> a_gpu = gpuarray.to_gpu(a)
>>> U, s, Vt = rlinalg.rsvd(a_gpu, k=2, method='standard')
"""
#*************************************************************************
#*** Author: N. Benjamin Erichson <nbe@st-andrews.ac.uk> ***
#*** <September, 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
copy_func = cublas.cublasCcopy
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
copy_func = cublas.cublasScopy
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
copy_func = cublas.cublasZcopy
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
copy_func = cublas.cublasDcopy
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 = np.array(a_gpu.shape, int)
if n>m : raise ValueError('input matrix of shape (m,n), where n>m is not supported')
#Set k
if k == None : raise ValueError('k must be provided')
if k > n or k < 1: raise ValueError('k must be 0 < k <= n')
kt = k
k = k + p
if k > n: k=n
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Generate a random sampling matrix O
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if isreal==False:
Oimag_gpu = gpuarray.empty((n,k), real_type, order="F", allocator=alloc)
Oreal_gpu = gpuarray.empty((n,k), real_type, order="F", allocator=alloc)
O_gpu = gpuarray.empty((n,k), data_type, order="F", allocator=alloc)
rand.fill_uniform(Oimag_gpu)
rand.fill_uniform(Oreal_gpu)
O_gpu = Oreal_gpu + 1j * Oimag_gpu
O_gpu = O_gpu.T * 2 - 1 #Scale to [-1,1]
else:
O_gpu = gpuarray.empty((n,k), real_type, order="F", allocator=alloc)
rand.fill_uniform(O_gpu) #Draw random samples from a ~ Uniform(-1,1) distribution
O_gpu = O_gpu * 2 - 1 #Scale to [-1,1]
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Build sample matrix Y : Y = A * O
#Note: Y should approximate the range of A
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Allocate Y
Y_gpu = gpuarray.zeros((m,k), data_type, order="F", allocator=alloc)
#Dot product Y = A * O
cublas_func_gemm(handle, 'n', 'n', m, k, n, alpha,
a_gpu.gpudata, m, O_gpu.gpudata, n,
beta, Y_gpu.gpudata, m )
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Orthogonalize Y using economic QR decomposition: Y=QR
#If q > 0 perfrom q subspace iterations
#Note: economic QR just returns Q, and destroys Y_gpu
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if q > 0:
Z_gpu = gpuarray.empty((n,k), data_type, order="F", allocator=alloc)
for i in np.arange(1, q+1 ):
if( (2*i-2)%q == 0 ):
Y_gpu = linalg.qr(Y_gpu, 'economic', lib='cula')
cublas_func_gemm(handle, TRANS_type, 'n', n, k, m, alpha,
a_gpu.gpudata, m, Y_gpu.gpudata, m,
beta, Z_gpu.gpudata, n )
if( (2*i-1)%q == 0 ):
Z_gpu = linalg.qr(Z_gpu, 'economic', lib='cula')
cublas_func_gemm(handle, 'n', 'n', m, k, n, alpha,
a_gpu.gpudata, m, Z_gpu.gpudata, n,
beta, Y_gpu.gpudata, m )
#End for
#End if
Q_gpu = linalg.qr(Y_gpu, 'economic', lib='cula')
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Project the data matrix a into a lower dimensional subspace
#B = Q.T * A
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Allocate B
B_gpu = gpuarray.empty((k,n), data_type, order="F", allocator=alloc)
cublas_func_gemm(handle, TRANS_type, 'n', k, n, m, alpha,
Q_gpu.gpudata, m, a_gpu.gpudata, m,
beta, B_gpu.gpudata, k )
if method == 'standard':
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Singular Value Decomposition
#Note: B = U" * S * Vt
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#gesvd(jobu, jobvt, m, n, int(a), lda, int(s), int(u), ldu, int(vt), ldvt)
#Allocate s, U, Vt for economic SVD
#Note: singular values are always real
s_gpu = gpuarray.empty(k, real_type, order="F", allocator=alloc)
U_gpu = gpuarray.empty((k,k), data_type, order="F", allocator=alloc)
Vt_gpu = gpuarray.empty((k,n), data_type, order="F", allocator=alloc)
#Economic SVD
cula_func_gesvd('S', 'S', k, n, int(B_gpu.gpudata), k, int(s_gpu.gpudata),
int(U_gpu.gpudata), k, int(Vt_gpu.gpudata), k)
#Compute right singular vectors as U = Q * U"
cublas_func_gemm(handle, 'n', 'n', m, k, k, alpha,
Q_gpu.gpudata, m, U_gpu.gpudata, k,
beta, Q_gpu.gpudata, m )
U_gpu = Q_gpu #Set pointer
# Free internal CULA memory:
cula.culaFreeBuffers()
#Return
return U_gpu[ : , 0:kt ], s_gpu[ 0:kt ], Vt_gpu[ 0:kt , : ]
elif method == 'fast':
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Orthogonalize B.T using reduced QR decomposition: B.T = Q" * R"
#Note: reduced QR returns Q and R, and destroys B_gpu
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if isreal==True:
B_gpu = transpose(B_gpu) #transpose B
else:
B_gpu = hermitian(B_gpu) #transpose B
Qstar_gpu, Rstar_gpu = linalg.qr(B_gpu, 'reduced', lib='cula')
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Singular Value Decomposition of R"
#Note: R" = U" * S" * Vt"
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#gesvd(jobu, jobvt, m, n, int(a), lda, int(s), int(u), ldu, int(vt), ldvt)
#Allocate s, U, Vt for economic SVD
#Note: singular values are always real
s_gpu = gpuarray.empty(k, real_type, order="F", allocator=alloc)
Ustar_gpu = gpuarray.empty((k,k), data_type, order="F", allocator=alloc)
Vtstar_gpu = gpuarray.empty((k,k), data_type, order="F", allocator=alloc)
#Economic SVD
cula_func_gesvd('A', 'A', k, k, int(Rstar_gpu.gpudata), k, int(s_gpu.gpudata),
int(Ustar_gpu.gpudata), k, int(Vtstar_gpu.gpudata), k)
#Compute right singular vectors as U = Q * Vt.T"
cublas_func_gemm(handle, 'n', TRANS_type, m, k, k, alpha,
Q_gpu.gpudata, m, Vtstar_gpu.gpudata, k,
beta, Q_gpu.gpudata, m )
U_gpu = Q_gpu #Set pointer
#Compute left singular vectors as Vt = U".T * Q".T
Vt_gpu = gpuarray.empty((k,n), data_type, order="F", allocator=alloc)
cublas_func_gemm(handle, TRANS_type, TRANS_type, k, n, k, alpha,
Ustar_gpu.gpudata, k, Qstar_gpu.gpudata, n,
beta, Vt_gpu.gpudata, k )
# Free internal CULA memory:
cula.culaFreeBuffers()
#Return
return U_gpu[ : , 0:kt ], s_gpu[ 0:kt ], Vt_gpu[ 0:kt , : ]
#End if
[docs]def rdmd(a_gpu, k=None, p=5, q=1, modes='exact', method_rsvd='standard', return_amplitudes=False, return_vandermonde=False, handle=None):
"""
Randomized 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.
p : int
`p` sets the oversampling parameter for rSVD (default k=5).
q : int
`q` sets the number of power iterations for rSVD (default=1).
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`.
method_rsvd : `{'standard', 'fast'}`
'standard' : (default) Standard algorithm as described in [1, 2]
'fast' : Version II algorithm as described in [2]
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
----------
N. B. Erichson and C. Donovan.
"Randomized Low-Rank Dynamic Mode Decomposition for Motion Detection"
Under Review.
N. Halko, P. Martinsson, and J. Tropp.
"Finding structure with randomness: probabilistic
algorithms for constructing approximate matrix
decompositions" (2009).
(available at `arXiv <http://arxiv.org/abs/0909.4061>`_).
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()
>>> rlinalg.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 = rlinalg.rdmd(F_gpu, k=2, p=0, q=1, 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 = np.array(a_gpu.shape, int)
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:]
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Randomized Singular Value Decomposition
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
U_gpu, s_gpu, Vh_gpu = rsvd(X_gpu, k=k, p=p, q=q,
method=method_rsvd, handle=handle)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#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 = linalg.eig(M_gpu, 'N', 'V', 'F', lib='cula')
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 = linalg.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
[docs]def cdmd(a_gpu, k=None, c=None, modes='exact', return_amplitudes=False, return_vandermonde=False, handle=None):
"""
Compressed 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.
c : int
`p` sets the number of measurements sensors.
modes : `{'exact'}`
'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
----------
S. L. Brunton, et al.
"Compressed sampling and dynamic mode decomposition."
arXiv preprint arXiv:1312.5186 (2013).
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()
>>> rlinalg.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 = rlinalg.cdmd(F_gpu, k=2, c=20, 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 = np.array(a_gpu.shape, int)
nx = n-1
#Set k
if k == None : k = nx
if k > nx or k < 1: raise ValueError('k is not valid')
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Compress
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if c==None:
Ac_gpu = A
c=m
else:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Generate a random sensing matrix S
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if isreal==False:
Simag_gpu = gpuarray.empty((m,c), real_type, order="F", allocator=alloc)
Sreal_gpu = gpuarray.empty((m,c), real_type, order="F", allocator=alloc)
S_gpu = gpuarray.empty((c,m), data_type, order="F", allocator=alloc)
rand.fill_uniform(Simag_gpu)
rand.fill_uniform(Sreal_gpu)
S_gpu = Sreal_gpu + 1j * Simag_gpu
S_gpu = S_gpu.T * 2 -1 #Scale to [-1,1]
else:
S_gpu = gpuarray.empty((c,m), real_type, order="F", allocator=alloc)
rand.fill_uniform(S_gpu) #Draw random samples from a ~ Uniform(-1,1) distribution
S_gpu = S_gpu * 2 - 1 #Scale to [-1,1]
#Allocate Ac
Ac_gpu = gpuarray.empty((c,n), data_type, order="F", allocator=alloc)
#Compress input matrix
cublas_func_gemm(handle, 'n', 'n', c, n, m, alpha,
int(S_gpu.gpudata), c, int(a_gpu.gpudata), m,
beta, int(Ac_gpu.gpudata), c )
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#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((c, n), data_type, order="F", allocator=alloc)
copy_func(handle, X_gpu.size, int(Ac_gpu.gpudata), 1, int(X_gpu.gpudata), 1)
X_gpu = X_gpu[:, :nx]
Y_gpu = Ac_gpu[:, 1:]
Yorig_gpu = a_gpu[:, 1:]
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Singular Value Decomposition
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Allocate s, U, Vt for economic SVD
#Note: singular values are always real
min_s = min(nx,c)
s_gpu = gpuarray.zeros(min_s, real_type, order="F", allocator=alloc)
U_gpu = gpuarray.zeros((c,min_s), data_type, order="F", allocator=alloc)
Vh_gpu = gpuarray.zeros((min_s,nx), data_type, order="F", allocator=alloc)
#Economic SVD
cula_func_gesvd('S', 'S', c, nx, int(X_gpu.gpudata), c, int(s_gpu.gpudata),
int(U_gpu.gpudata), c, int(Vh_gpu.gpudata), min_s)
#Low-rank DMD: trancate SVD if k < nx
if k != nx:
s_gpu = s_gpu[:k]
U_gpu = U_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.zeros((c,k), data_type, order="F", allocator=alloc)
M_gpu = gpuarray.zeros((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/s_gpu
#ii) real/complex: scale Vs = Vt* x diag(s**-1)
Vs_gpu = gpuarray.zeros((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,
1.0, int(Vh_gpu.gpudata), lda, 0.0, int(Vh_gpu.gpudata), lda,
int(Vs_gpu.gpudata), ldb)
#End Transpose
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', c, k, nx, alpha,
int(Y_gpu.gpudata), c, int(Vs_gpu.gpudata), nx,
beta, int(G_gpu.gpudata), c )
#iv) real/complex: M = U* x G
cublas_func_gemm(handle, TRANS_type, 'n', k, k, c, alpha,
int(U_gpu.gpudata), c, int(G_gpu.gpudata), c,
beta, int(M_gpu.gpudata), k )
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Eigen Decomposition
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Note: If a_gpu is real the imag part is omitted
Vr_gpu, w_gpu = linalg.eig(M_gpu, 'N', 'V', 'F', lib='cula')
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', nx, k, k, alpha,
int(Vs_gpu.gpudata), nx, int(Vr_gpu.gpudata), k,
beta, int(Vs_gpu.gpudata), nx )
cublas_func_gemm(handle, 'n', 'n', m, k, nx, alpha,
Yorig_gpu.gpudata, m, Vs_gpu.gpudata, nx,
beta, F_gpu.gpudata, m )
else:
raise ValueError('Type of modes is not supported, choose "exact" or "standard".')
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Compute amplitueds b using least-squares: Fb=x1
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if return_amplitudes==True:
F_gpu_temp = gpuarray.empty((m,k), data_type, order="F", allocator=alloc)
#Copy is required, because gels destroys input
copy_func(handle, F_gpu.size, int(F_gpu.gpudata),
1, int(F_gpu_temp.gpudata), 1)
#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 = linalg.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
if __name__ == "__main__":
import doctest
doctest.testmod()