LAPACKE C ++ Real Matrix Inversion - PullRequest
Новая
       75

LAPACKE C ++ Real Matrix Inversion

1 голос
/ 24 января 2020

Я пытаюсь инвертировать реальную матрицу в 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 обратный (python прямо здесь, потому что я умножил матрицу на обратное, и результат очень близок к Indentity). В выводе LAPACKE я умножаю матрицу на ее обратную и получаю диагонали, равные 1, но не диагонали go до значений до 0,17, что огромно по сравнению с 0. Есть ли способ заставить программу LAPACKE дать лучший результат? Спасибо!

1 Ответ

0 голосов
/ 25 января 2020

Для матрицы с большим детерминантом вы могли бы изменить масштаб ввода, вычислить обратное значение, а затем изменить масштаб вывода обратно. Вот очень простой 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];
    }
}
...