1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
|
#!/usr/bin/env python
"""
Demo of how to call low-level CUSOLVER wrappers to perform eigen decomposition
for Hermitian matrices.
"""
import numpy as np
import pycuda.autoinit
import pycuda.gpuarray as gpuarray
import skcuda.cusolver as solver
handle = solver.cusolverDnCreate()
x = np.random.randn(1024,1024)+1j*np.random.rand(1024,1024)
x = x+x.conj().T
# Need to reverse dimensions because CUSOLVER expects column-major matrices:
n, m = x.shape
x_gpu = gpuarray.to_gpu(x.T.copy())
# Set up output buffers:
w = gpuarray.empty(n, dtype = np.double)
# Set up parameters
params = solver.cusolverDnCreateSyevjInfo()
solver.cusolverDnXsyevjSetTolerance(params, 1e-7)
solver.cusolverDnXsyevjSetMaxSweeps(params, 15)
# Set up work buffers:
lwork = solver.cusolverDnZheevj_bufferSize(handle, 'CUSOLVER_EIG_MODE_VECTOR',
'u', n, x_gpu.gpudata, m,
w.gpudata, params)
workspace_gpu = gpuarray.zeros(lwork, dtype = x.dtype)
info = gpuarray.zeros(1, dtype = np.int32)
# Compute:
solver.cusolverDnZheevj(handle, 'CUSOLVER_EIG_MODE_VECTOR',
'u', n, x_gpu.gpudata, m,
w.gpudata, workspace_gpu.gpudata,
lwork, info.gpudata, params)
# Print info
print(solver.cusolverDnXsyevjGetSweeps(handle, params))
print(solver.cusolverDnXsyevjGetResidual(handle, params))
# Destroy handle
solver.cusolverDnDestroySyevjInfo(params)
solver.cusolverDnDestroy(handle)
# Check error
Q = x_gpu.get().T
print('maximum error in A * Q - Q * Lambda is: %r' %
np.abs(np.dot(x, Q) - np.dot(Q, np.diag(w.get()))).max())
|
