Skip to content

Commit a6037d8

Browse files
wingkitlee0wingkitlee0
committed
added ssyedx and ssyedx_m (and demo), still need more tests and docs
1 parent 15a69af commit a6037d8

File tree

2 files changed

+184
-28
lines changed

2 files changed

+184
-28
lines changed

‎demos/magma/demo_ssyedx_m.py‎

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
"""
2+
Demo of using ssyedx_m
3+
"""
4+
import numpy as np
5+
from skcuda import magma
6+
import time
7+
8+
typedict = {'s': np.float32, 'd': np.float64, 'c': np.complex64, 'z': np.complex128}
9+
typedict_= {v: k for k, v in typedict.items()}
10+
11+
def eigs(a_gpu, k=None, which='LM', imag=False, return_eigenvectors=True):
12+
"""
13+
Driver for eigenvalue solver for symmetric matrix
14+
"""
15+
16+
if len(a_gpu.shape) != 2:
17+
raise ValueError("M needs to be a rank 2 square array for eig.")
18+
19+
magma.magma_init()
20+
ngpu=1
21+
22+
dtype = a_gpu.dtype.type
23+
t = typedict_[dtype]
24+
N = a_gpu.shape[1]
25+
26+
if k is None: k=int(np.ceil(N/2))
27+
28+
if 'L' in which:
29+
a_gpu = -a_gpu
30+
31+
if return_eigenvectors:
32+
jobz = 'V'
33+
else:
34+
jobz = 'N'
35+
36+
rnge = 'I'; uplo = 'U'
37+
vl = 0; vu = 1e10
38+
il = 1; iu = k; m = 0
39+
print(f"k = {k}, iu = {iu}")
40+
41+
w_gpu = np.empty((N,), dtype, order='F') # eigenvalues
42+
43+
if t == 's':
44+
nb = magma.magma_get_ssytrd_nb(N)
45+
elif t == 'd':
46+
nb = magma.magma_get_dsytrd_nb(N)
47+
else:
48+
raise ValueError('unsupported type')
49+
50+
lwork = N*(1 + 2*nb)
51+
if jobz:
52+
lwork = max(lwork, 1+6*N+2*N**2)
53+
work = np.empty(lwork, dtype)
54+
liwork = 3+5*N if jobz else 1
55+
iwork = np.empty(liwork, dtype)
56+
57+
if t == 's':
58+
status = magma.magma_ssyevdx_m(ngpu, jobz, rnge, uplo, N, a_gpu.ctypes.data, N,
59+
vl, vu, il, iu, m,
60+
w_gpu.ctypes.data, work.ctypes.data, lwork, iwork.ctypes.data, liwork)
61+
elif t == 'd':
62+
status = magma.magma_dsyevdx_m(ngpu, jobz, rnge, uplo, N, a_gpu.ctypes.data, N,
63+
vl, vu, il, iu, m,
64+
w_gpu.ctypes.data, work.ctypes.data, lwork, iwork.ctypes.data, liwork)
65+
else:
66+
raise ValueError('unsupported type')
67+
68+
if 'L' in which:
69+
w_gpu = -w_gpu
70+
71+
magma.magma_finalize()
72+
73+
if jobz:
74+
return w_gpu, a_gpu
75+
else:
76+
return w_gpu
77+
78+
79+
if __name__=='__main__':
80+
import sys
81+
N = int(sys.argv[1])
82+
83+
# not symmetric, but only side of the diagonal is used
84+
M_gpu = np.random.random((N, N))
85+
M_gpu = M_gpu.astype(np.float32)
86+
M_cpu = M_gpu.copy()
87+
88+
# GPU
89+
t1 = time.time()
90+
W_gpu, V_gpu = eigs(M_gpu)
91+
W_gpu = np.sort(W_gpu)
92+
t2 = time.time()
93+
94+
# CPU
95+
t3 = time.time()
96+
W_cpu, V_cpu = np.linalg.eigh(M_cpu)
97+
W_cpu = np.sort(W_cpu)
98+
t4 = time.time()
99+
100+
print("First 10 eigenvalues")
101+
print(f"GPU: {W_gpu[:10]}")
102+
print(f"CPU: {W_cpu[:10]}")
103+
print("Time")
104+
print(f"GPU: {t2-t1}")
105+
print(f"CPU: {t4-t3}")

‎skcuda/magma.py‎

Lines changed: 79 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -3114,12 +3114,12 @@ def magma_zposv_gpu(uplo, n, nhrs, a_gpu, lda, b_gpu, ldb):
31143114

31153115
# SGESV, DGESV, CGESV, ZGESV
31163116
_libmagma.magma_sgesv_gpu.restype = int
3117-
_libmagma.magma_sgesv_gpu.argtypes = [c_int_type,
3118-
c_int_type,
3119-
ctypes.c_void_p,
3120-
c_int_type,
3121-
ctypes.c_void_p,
3122-
ctypes.c_void_p,
3117+
_libmagma.magma_sgesv_gpu.argtypes = [c_int_type, # n
3118+
c_int_type, # nrhs
3119+
ctypes.c_void_p, # dA
3120+
c_int_type, # lddA
3121+
ctypes.c_void_p, # ipiv
3122+
ctypes.c_void_p, # dB
31233123
c_int_type,
31243124
ctypes.c_void_p]
31253125
def magma_sgesv_gpu(n, nhrs, A, lda, ipiv, B, ldb):
@@ -3716,40 +3716,91 @@ def magma_zheevd_m(ngpu, jobz, uplo, n, A, lda, w, work, lwork,
37163716
lrwork, int(iwork), liwork, ctypes.byref(info))
37173717
magmaCheckStatus(status)
37183718

3719+
# SSYEVDX
3720+
# TO BE TESTED
3721+
_libmagma.magma_ssyevdx.restype = int
3722+
_libmagma.magma_ssyevdx.argtypes = [ c_int_type, # jobz
3723+
c_int_type, # rnge
3724+
c_int_type, # uplo
3725+
c_int_type, # n
3726+
ctypes.c_void_p, # A
3727+
c_int_type, # lda
3728+
ctypes.c_float, # vl
3729+
ctypes.c_float, # vu
3730+
c_int_type, # il
3731+
c_int_type, # iu
3732+
ctypes.c_void_p, # m
3733+
ctypes.c_void_p, # w
3734+
ctypes.c_void_p, # work
3735+
c_int_type, # lwork
3736+
ctypes.c_void_p, # iwork
3737+
c_int_type, # liwork
3738+
ctypes.c_void_p] # info
3739+
def magma_ssyevdx(jobz, rnge, uplo, n, A, lda,
3740+
vl, vu, il, iu, m,
3741+
w, work, lwork, iwork, liwork):
3742+
"""
3743+
Compute eigenvalues of real symmetric matrix.
3744+
Multi-GPU, data on host
3745+
"""
3746+
3747+
# _XXX_conversion[] returns integer according to magma_types.h
3748+
jobz = _vec_conversion[jobz]
3749+
rnge = _range_conversion[rnge]
3750+
uplo = _uplo_conversion[uplo]
3751+
3752+
mfound = c_int_type()
3753+
info = c_int_type()
3754+
status = _libmagma.magma_ssyevdx(jobz, rnge, uplo, n, int(A), lda,
3755+
vl, vu, il, iu, ctypes.byref(mfound),
3756+
int(w), int(work), lwork, int(iwork), liwork,
3757+
ctypes.byref(info))
3758+
m = mfound
3759+
magmaCheckStatus(status)
37193760

37203761
# SSYEVDX_M, DSYEVDX_M, CHEEVDX_M, ZHEEVDX_M
3721-
_libmagma.magma_ssyevd_m.restype = int
3722-
_libmagma.magma_ssyevd_m.argtypes = [c_int_type,
3723-
c_int_type,
3724-
c_int_type,
3725-
c_int_type,
3726-
c_int_type,
3727-
ctypes.c_void_p,
3728-
c_int_type,
3729-
ctypes.c_float,
3730-
ctypes.c_float,
3731-
c_int_type,
3732-
c_int_type,
3733-
ctypes.c_void_p,
3734-
ctypes.c_void_p,
3735-
ctypes.c_void_p,
3736-
c_int_type,
3737-
ctypes.c_void_p,
3738-
c_int_type,
3739-
ctypes.c_void_p]
3762+
# TO BE TESTED
3763+
_libmagma.magma_ssyevdx_m.restype = int
3764+
_libmagma.magma_ssyevdx_m.argtypes = [c_int_type, # ngpu
3765+
c_int_type, # jobz
3766+
c_int_type, # rnge
3767+
c_int_type, # uplo
3768+
c_int_type, # n
3769+
ctypes.c_void_p, # A
3770+
c_int_type, # lda
3771+
ctypes.c_float, # vl
3772+
ctypes.c_float, # vu
3773+
c_int_type, # il
3774+
c_int_type, # iu
3775+
ctypes.c_void_p, # m
3776+
ctypes.c_void_p, # w
3777+
ctypes.c_void_p, # work
3778+
c_int_type, # lwork
3779+
ctypes.c_void_p, # iwork
3780+
c_int_type, # liwork
3781+
ctypes.c_void_p] # info
37403782
def magma_ssyevdx_m(ngpu, jobz, rnge, uplo, n, A, lda,
37413783
vl, vu, il, iu, m,
37423784
w, work, lwork, iwork, liwork):
37433785
"""
37443786
Compute eigenvalues of real symmetric matrix.
37453787
Multi-GPU, data on host
3788+
3789+
source: ssyedx_m.cpp
37463790
"""
37473791

3792+
# _XXX_conversion[] returns integer according to magma_types.h
3793+
jobz = _vec_conversion[jobz]
3794+
rnge = _range_conversion[rnge]
37483795
uplo = _uplo_conversion[uplo]
3796+
3797+
mfound = c_int_type() # output
37493798
info = c_int_type()
3750-
status = _libmagma.magma_ssyevdx_m(ngpu, jobz, uplo, n, int(A), lda,
3751-
int(w), int(work),
3752-
lwork, int(iwork), liwork, ctypes.byref(info))
3799+
status = _libmagma.magma_ssyevdx_m(ngpu, jobz, rnge, uplo, n, int(A), lda,
3800+
vl, vu, il, iu, ctypes.byref(mfound),
3801+
int(w), int(work), lwork, int(iwork), liwork,
3802+
ctypes.byref(info))
3803+
m = mfound
37533804
magmaCheckStatus(status)
37543805

37553806
# SYMMETRIZE

0 commit comments

Comments
 (0)