Coming from the as a reference: How to include LAPACK functions with complex inputs · Issue #5301 · numba/numba · GitHub
I have the following attempt at calling dsyevr from scipy using numba:
import ctypes
import numpy as np
import numba as nb
from numba import vectorize, njit, jit
from numba.extending import get_cython_function_address
addr = get_cython_function_address('scipy.linalg.cython_lapack','dsyevr')
_PTR, _dble, _char, _int = ctypes.POINTER, ctypes.c_double, ctypes.c_char, ctypes.c_int
_ptr_select = ctypes.c_voidp
_ptr_dble = _PTR(_dble)
_ptr_char = _PTR(_char)
_ptr_int = _PTR(_int)
## Retrieve function address
functype = ctypes.CFUNCTYPE(None,
_ptr_int, # JOBVS
_ptr_int, # RANGE
_ptr_int, # UPLO
_ptr_int, # N
_ptr_dble, # A
_ptr_int, # LDA
_ptr_dble, # VL
_ptr_dble, # VU
_ptr_int, # IL
_ptr_int, # IU
_ptr_dble, # ABSTOL
_ptr_int, # M
_ptr_dble, # W
_ptr_dble, # Z
_ptr_int, # LDZ
_ptr_int, # ISUPPZ
_ptr_dble, # WORK
_ptr_int, # LWORK
_ptr_int, # IWORK
_ptr_int, # LIWORK
_ptr_int, # INFO
)
dsyevr_fn = functype(addr)
_ORD_JOBVS = ord('N')
_ORD_RNG = ord('I')
_ORD_UPLO = ord('U')
@nb.jit(nb.types.Tuple((nb.float64[:], nb.float64[:,:], nb.int32))(nb.float64[:,:], nb.int32, nb.int32), nopython=True)
def numba_dsyevr(x, I1, I2):
assert (x.shape[0] == x.shape[1]) and x[0,0] == 0.0
assert I1 <= I2 and I1 >= 1 and I2 <= x.shape[0]
## Setup arguments
i1, i2 = I1, I2
n, m = x.shape[0], abs(i2-i1)+1
JOBVS = np.array([_ORD_JOBVS], np.int32)
RNG = np.array([_ORD_RNG], np.int32)
UPLO = np.array([_ORD_UPLO], np.int32)
N = np.array([n], np.int32)
A = x.copy() # in & out
LDA = np.array([n], np.int32)
VL, VU, IL, IU = np.array(0, np.float64), np.array(0, np.float64), np.array(i1, np.int32), np.array(i2, np.int32)
ABS_TOL = np.array([0.0], np.float64)
M = np.array([m], np.int32)
W, Z, LDZ, ISUPPZ, WORK, LWORK = np.zeros(n, np.float64), np.zeros((n, m), np.float64).T, np.array([n], np.int32), np.zeros(2*m, np.int32), np.zeros(26*n, np.float64), np.array([26*n], np.int32)
IWORK, LIWORK = np.zeros(10*n, np.int32), np.array([10*n], np.int32)
INFO = np.empty(1, dtype=np.int32)
## Call the function
# jobvs = ctypes.POINTER(ctypes.c_int)()
dsyevr_fn(
JOBVS.ctypes,
RNG.ctypes,
UPLO.ctypes,
N.ctypes,
A.view(np.float64).ctypes,
LDA.ctypes,
VL.ctypes,
VU.ctypes,
IL.ctypes,
IU.ctypes,
ABS_TOL.ctypes,
M.ctypes,
W.view(np.float64).ctypes,
Z.view(np.float64).ctypes,
LDZ.ctypes,
ISUPPZ.view(np.int32).ctypes,
WORK.view(np.float64).ctypes,
LWORK.ctypes,
IWORK.view(np.int32).ctypes,
LIWORK.ctypes,
INFO.ctypes
)
return((W, Z.T, INFO[0]))
I’m trying to do this to improve an MDS calculation. The function compiles, returns INFO=0, and I have the tolerance down to 0.0 (which the docs say you can do). I’ve tripled checked the types, the sizes of the inputs, etc… nonetheless, I’m getting negative eigenvalues w/ the following:
import numpy as np
from scipy.spatial.distance import squareform, pdist, cdist
a = np.random.uniform(size=(15,8))
x = squareform(pdist(a))**2
z = numba_dsyevr(x, 1, x.shape[0]-1)
Does anyone have any idea why this is happening or possible things I can do to debug this? With INFO returning 0 I’m not sure how to approach