4

我试图反转的矩阵是：

    [ 1  0  1]
A = [ 2  0  1]
    [-1  1  1]

真正的逆是：

       [-1  1  0]
A^-1 = [-3  2  1]
       [ 2 -1  0]

使用 Python 的 numpy.linalg.inv，我得到了正确的答案。我的矩阵逆例程之一使用 dgetri_，它是：

void compute_matrix_inverse_dbl( double* matrix,
                                 int order,
                                 double * inverse )
{

    int N, lwork;
    int success;
    int *pivot;
    double* workspace;

    //===Allocate Space===//
    pivot = malloc(order * order * order * sizeof(*pivot));
    workspace = malloc(order * order * sizeof(*workspace));

    //===Run Setup===//
    N = order;
    copy_array_dbl(matrix, order*order, inverse);
    lwork = order*order;

    //===Factor Matrix===//
    dgetrf_(&N,&N,inverse,&N,pivot,&success);

    //===Compute Inverse===//
    dgetri_(&N, inverse, &N, pivot, workspace, &lwork, &success);

    //===Clean Up===//
    free(workspace);
    free(pivot);

    return;
  }

使用这个例程，我得到：

       [-1   1 +-e1 ]
A^-1 = [-3   2   1  ]
       [ 2  -1 +-e2 ]

其中 e1 和 e2 以及机器精度为 1e-16 的小数。现在也许 dgetri_ 不是最好用的。但是，当我通过 zgeqrf_ 和 zungqr_ 使用 QR 分解进行反转时，我得到了类似的答案。当我使用 dgesvd_ 使用 SVD 进行逆时，我也得到了类似的答案。

Python 似乎使用了一个名为 _umath_linalg.inv 的例程。所以我有几个问题：

  • 那个例程是做什么的？
  • 我可以使用什么 CBLAS/LAPACK 例程来反转此矩阵并获得类似 CBLAS/LAPACK 的结果（这样 e1 和 e2 被适当的零替换）？
4

1 回答 1

1

根据描述，这似乎是scipy.linalg.invnumpy.linalg.inv的精简版：

该模块是 SciPy 中 linalg.py 模块的精简版，其中包含 LAPACK 库的高级 Python 接口。

查看scipy.linalg.inv，它调用getrf, then getri。

于 2018-10-30T16:54:28.180 回答