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

概述我在这里定义了计算矩阵行列式的函数.但有时我会得到错误的信号.我从 this answer模拟了我的功能. from scipy.linalg.cython_lapack cimport dgetrfcpdef double det_c(double[:, ::1] A, double[:, ::1] work, double[::1] ipiv): '''obtain determi 我在这里定义了计算矩阵行列式的函数.但有时我会得到错误的信号.我从 this answer模拟了我的功能.

from scipy.linalg.cython_lapack cimport dgetrfcpdef 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

解决方法 dgetrf是Fortran子例程,Fortran使用基于1的索引,因此ipiv中的值介于1和n之间(包括1和n).为了解决这个问题,请更改循环中的测试

if j != ipiv[j]:

至

if j != ipiv[j] - 1:

总结

以上是内存溢出为你收集整理的python – 在Cython中使用LAPACK包装器估算行列式的LU分解全部内容，希望文章能够帮你解决python – 在Cython中使用LAPACK包装器估算行列式的LU分解所遇到的程序开发问题。

如果觉得内存溢出网站内容还不错，欢迎将内存溢出网站推荐给程序员好友。