Инверсия матрицы в CBLAS/LAPACK против Python

Question

Матрица, которую я пытаюсь инвертировать:

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

Истинное обратное:

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

Используя Python's 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_ не самый лучший в использовании. Однако, когда я инвертирую, используя QR-разложение через zgeqrf_ и zungqr_, я получаю аналогичный ответ. Когда я использую dgesvd_ для инверсии с использованием SVD, я получаю аналогичный ответ.

Кажется, что Python использует подпрограмму _umath_linalg.inv. Итак, у меня есть несколько вопросов:

• Что делает эта рутина?
• Какую процедуру CBLAS/LAPACK я могу использовать, чтобы инвертировать эту матрицу и получить результат, подобный CBLAS/LAPACK (такой, что e1 и e2 заменяются правильными нулями)?

Answer 1

Кажется, что numpy.linalg.inv облегченная версия scipy.linalg.inv согласно описанию:

Этот модуль является облегченной версией модуля linalg.py в SciPy, который содержит высокоуровневый интерфейс Python для библиотеки LAPACK.

Глядя на scipy.linalg.inv, он звонит getrf, затем getri,