python – 在Cython中使用LAPACK包装器估算行列式的LU分解

from scipy.linalg.cython_lapack cimport dgetrf

cpdef double det_c(double[:,::1] A,double[:,::1] work,double[::1] ipiv):
    '''obtain determinant of float type square matrix A

    Notes
    -----
    As is,this function is not yet computing the sign of the determinant
    correctly,help!

    Parameters
    ----------
    A : memoryview (numpy array)
        n x n array to compute determinant of
    work : memoryview (numpy array)
        n x n array to use within function
    ipiv : memoryview (numpy array)
        length n vector use within function

    Returns
    -------
    detval : float
        determinant of matrix A
    '''

    cdef int n = A.shape[0],info
    work[...] = A

    dgetrf(&n,&n,&work[0,0],&ipiv[0],&info)

    cdef double detval = 1.
    cdef int j

    for j in range(n):
        if j != ipiv[j]:
            detval = -detval*work[j,j]
        else:
            detval = detval*work[j,j]

    return detval

当我测试这个函数并将它与np.linalg.det进行比较时,有时我得到了错误的符号.

>>> a = np.array([[1,2],[3,5.]])
>>> np.linalg.det(a)
>>> -1.0000000000000004
>>> det_c(a,np.zeros((2,2)),np.zeros(2,dtype=np.int32))
>>> 1

其他时候,正确的标志.

>>> b = np.array([[1,2,3],[1,1],[5,6,1.]])
>>> np.linalg.det(b)
>>> -7.999999999999998
>>> det_c(a,np.zeros((3,3)),np.zeros(3,dtype=np.int32))
>>> -8.0