LAPACKE C++ Обращение вещественной матрицы

Я пытаюсь инвертировать реальную матрицу в С++ LAPACKE. У меня есть такая же функция для сложных матриц, и она работает. Но реальный случай дает неверный ответ. Вот моя функция:

void inv(std::vector<std::vector<double>> &ans, std::vector<std::vector<double>> MAT){

    int N = MAT.size();

    int *IPIV = new int[N];

    double * arr = new double[N*N];
    for (int i = 0; i<N; i++){
        for (int j = 0; j<N; j++){
            int idx = i*N + j;
            arr[idx] = MAT[i][j];
        }
    }

    LAPACKE_dgetrf(LAPACK_ROW_MAJOR, N, N, arr, N, IPIV);
    LAPACKE_dgetri(LAPACK_ROW_MAJOR, N, arr, N, IPIV);

     for (int i = 0; i<N; i++){
        for (int j = 0; j<N; j++){
            int idx = i*N + j;
            ans[i][j] = arr[idx];
        }
    }

    delete[] IPIV;
    delete[] arr;
}

Я попытался инвертировать матрицу двойников 24 на 24. В то время как программа, кажется, почти готова, инверсия еще не совсем завершена, и она сильно отличается от того, что дает мне python linalg inverse (python здесь, потому что я умножил матрицу на инверсию, и результат очень близок к идентичности). В выводе LAPACKE я умножаю матрицу на обратную и получаю, что диагонали равны 1, но недиагонали достигают значений до 0,17, что огромно по сравнению с 0. Есть ли способ заставить программу LAPACKE давать лучшее результат? Спасибо!


c++ matrix matrix-inverse lapack lapacke
person Imran Alkhatib    schedule 24.01.2020    source источник
comment
Это четко определенная матрица? Чему равен определитель матрицы?   -  person Severin Pappadeux    schedule 24.01.2020
comment
2,2864066779666567e+37. У него большие значения, поэтому я думаю, почему у программы проблемы?   -  person Imran Alkhatib    schedule 24.01.2020
comment
Да, это может объяснить эффект. Вы можете попытаться умножить входную матрицу (каждый член) на (1/25), что уменьшит определитель в 25 ^ {- 24} = 2,8e-34 (если я не ошибаюсь), таким образом сделав определитель входной матрицы около 1000. Затем вычислите обратное значение и умножьте его на 1/25. Я добавил в ответ простой код Python.   -  person Severin Pappadeux    schedule 25.01.2020
comment
Взгляните на источник linalg.inv() numpy, в частности на строку 1555 - 1693 -1727 umath_linalg.c, он решает A.X=Id, используя dgesv() вместо вызова dgetrf() и dgetri(). Scipy использует dgetrf() / dgetri(). Не могли бы вы попробовать использовать scipy.linalg.inv() scipy в python и посмотреть, приведет ли это к правильному выводу?   -  person francis    schedule 25.01.2020

Ответы (1)


arrow_upward
-1
arrow_downward

Для матрицы с большим определителем вы можете изменить масштаб ввода, вычислить инверсию, а затем снова масштабировать вывод. Вот очень простой пример Python, ваш масштабный коэффициент должен быть 1/25 или около того, чтобы получить общий множитель (1/25)24=2,8e-34, что делает определитель входной матрицы около 1000 .

import numpy as np

scale = 0.5

i = np.array([[1,2],[3,4]])
print(i)
print(np.linalg.det(i))
print("-----------------------------------")
x = np.multiply(i, [scale]) # rescale matrix
print(x)
print(np.linalg.det(x))     # determinant should be less
print("-----------------------------------")
y = np.linalg.inv(x)
print(y)
print(np.linalg.det(y))
print("-----------------------------------")
o = np.multiply(y, [scale]) # rescale matrix
print(o)
print(np.linalg.det(o))
print(np.dot(i, o))

Вы можете поиграть с scale и убедиться, что любое значение кода возвращает правильно инвертированную входную матрицу.

Таким образом, ваш код будет

double scale = 1.0/25.0;

for (int i = 0; i<N; i++){
    for (int j = 0; j<N; j++){
        int idx = i*N + j;
        arr[idx] = scale*MAT[i][j];
    }
}

....

for (int i = 0; i<N; i++){
    for (int j = 0; j<N; j++){
        int idx = i*N + j;
        ans[i][j] = scale * arr[idx];
    }
}
person Severin Pappadeux    schedule 24.01.2020