I am trying to write Cython functions to compute inverses and obtain solutions to linear systems. I know this is what Numpy does. But I need use these functions in another loop I'm writing in Cython. And so, I don't want to use Python objects in the loop.
I've written an inverse and solve function, and they work some of the time. But other times they fail. I think I may be misunderstanding something fundamental. Anyway, here is a MWE. All of the files are in the same directory. First, the file mopp.pxd:
#cython: boundscheck=False, wraparound=False, nonecheck=False
from scipy.linalg.cython_lapack cimport dgetri, dgetrf, dgesv
cimport cython
cdef inline void solve_c(double[:, ::1] A, double[:, ::1] B, double[:, ::1] out,
int[::1] ipiv, double[:, ::1] Awork):
'''solve system A*X=B [WARNING: Not working for B with 2+ columns]
Parameters
----------
A : memoryview
n x n
B : memoryview
n x r
out : memoryview
n x r, for output storage
ipiv : memoryview
length n integer vector
Awork : memoryview
n x n vector to use for work
'''
cdef int LDA = A.shape[0], N = A.shape[1]
cdef int LDB = B.shape[0], NRHS = B.shape[1]
cdef int info
Awork[...] = A
out[...] = B
dgesv(&N, &NRHS, &Awork[0,0], &LDA, &ipiv[0], &out[0,0], &LDB, &info)
cdef inline void inv_c(double[:, ::1] A, double[:, ::1] B,
double[:, ::1] work, int[::1] ipiv):
'''invert float type square matrix A [WARNING: Not working randomly?]
Parameters
----------
A : memoryview (numpy array)
n x n array to invert
B : memoryview (numpy array)
n x n array to use within the function, function
will modify this matrix in place to become the inverse of A
work : memoryview (numpy array)
n x n array to use within the function
ipiv : memoryview (numpy array)
length n array to use within function
'''
cdef int n = A.shape[0], info, lwork
B[...] = A
dgetrf(&n, &n, &B[0, 0], &n, &ipiv[0], &info)
dgetri(&n, &B[0,0], &n, &ipiv[0], &work[0,0], &lwork, &info)
Then a file containing wrappers to make these functions available to Python for testing, called mopp_wrappers.pyx:
#cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True
from mopp cimport inv_c, solve_c
import numpy as np
def inv(A, B, work, ipiv):
inv_c(A, B, work, ipiv)
def solve(A, B, out, ipiv, Awork):
solve_c(A, B, out, ipiv, Awork)
And I have the associated setup file setup.py:
from distutils.core import setup
from Cython.Build import cythonize
setup(name="mopp_wrappers", ext_modules=cythonize('mopp_wrappers.pyx', annotate=True))
I then compile by running the following in the terminal:
python setup.py build_ext --inplace
Finally, I test the functions in Python with test.py. I just don't understand why the functions work in some scenarios but not others.
import numpy as np
from mopp_wrappers import inv, solve
def mysolve(A, b):
out = np.zeros((A.shape[0], b.shape[1]))
ipiv = np.zeros(3, dtype=np.int32)
Awork = np.zeros(A.shape)
solve(A, b, out, ipiv, Awork)
return out
def myinv(A):
B = np.zeros(A.shape)
work = np.zeros(A.shape)
ipiv = np.zeros(A.shape[0], dtype=np.int32)
inv(A, B, work, ipiv)
return B
np.random.seed(100)
M = np.random.multivariate_normal(mean=np.zeros(3), cov=np.eye(3), size=20)
A = M.T.dot(M)
# inverse is the same
print('Successful computation of inverse:')
print(np.linalg.inv(A))
print(myinv(A))
print('')
# inverse is different
B = np.zeros(A.shape)
work = np.zeros(A.shape)
ipiv = np.zeros(A.shape[0], dtype=np.int32)
inv(A, B, work, ipiv)
print('Failure to compute inverse:')
print(B)
print(np.linalg.inv(A))
print('')
# solution is the same
b = np.array([1,2,3]).reshape((-1, 1)).astype(float)
print('Successful computation of solution, 1-dim b matrix')
print(mysolve(A, b))
print(np.linalg.solve(A, b))
print('')
# solution is different
b = np.hstack([b, b])
print('Failure to compute solution, 2-dim b matrix')
print(mysolve(A, b))
print(np.linalg.solve(A, b))
where I get the following output:
Successful computation of inverse:
[[ 0.0587434 0.01003223 0.01867754]
[ 0.01003223 0.06387784 -0.00035565]
[ 0.01867754 -0.00035565 0.06977158]]
[[ 0.0587434 0.01003223 0.01867754]
[ 0.01003223 0.06387784 -0.00035565]
[ 0.01867754 -0.00035565 0.06977158]]
Failure to compute inverse:
[[ 1.91799368e+01 -1.58548362e-01 -2.68503736e-01]
[-3.04094756e+00 1.56553248e+01 5.09740667e-03]
[-5.14988469e+00 7.98015569e-02 1.43324831e+01]]
[[ 0.0587434 0.01003223 0.01867754]
[ 0.01003223 0.06387784 -0.00035565]
[ 0.01867754 -0.00035565 0.06977158]]
Successful computation of solution, 1-dim b matrix
[[0.13484049]
[0.13672096]
[0.22728098]]
[[0.13484049]
[0.13672096]
[0.22728098]]
Failure to compute solution, 2-dim b matrix
[[0.10613072 0.07319877]
[0.15786505 0.20361612]
[0.21063103 0.24560286]]
[[0.13484049 0.13484049]
[0.13672096 0.13672096]
[0.22728098 0.22728098]]
