Вычисление инверсии матрицы с использованием lapack в C
Я хотел бы иметь возможность вычислять обратную матрицу NxN в C / C ++ с использованием lapack.
Насколько я понимаю, способ сделать инверсию в лапаке - использовать функцию dgetri , однако я не могу понять, какими должны быть все ее аргументы.
Вот код, который у меня есть:
void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO);
int main(){
double M [9] = {
1,2,3,
4,5,6,
7,8,9
};
return 0;
}
Как бы вы заполнили его, чтобы получить обратную матрицу M 3x3 , используя dgetri _?
2 ответа
Во-первых, M должен быть двумерным массивом, например double M [3] [3] . Математически ваш массив представляет собой вектор 1x9, который не является обратимым.
N - указатель на int для порядок матрицы - в данном случае N = 3.
A - указатель на LU факторизация матрицы, которая вы можете получить, запустив LAPACK подпрограмма
dgetrf.LDA - целое число для "ведущего элемент »матрицы, позволяющий вы выбираете подмножество большего матрицу, если вы хотите просто инвертировать маленький кусочек. Если вы хотите инвертировать всю матрицу, LDA просто должно быть равным N.
IPIV - это основные индексы матрица, другими словами, это список инструкций, какие строки поменять местами чтобы перевернуть матрицу. IPIV должен быть сгенерирован LAPACK подпрограмма
dgetrf.LWORK и WORK - это "рабочие области" используется LAPACK. Если вы инвертируете всю матрицу, LWORK должен быть int равным N ^ 2, а РАБОТА должна быть двойной массив с элементами LWORK.
INFO - это просто статусная переменная для скажу вам, была ли операция завершено успешно. Поскольку не все матрицы обратимы, я бы рекомендую отправить это некоторым своего рода система проверки ошибок. INFO = 0 для успешной работы, INFO = -i, если i-й аргумент имел неправильное входное значение, и INFO> 0, если матрица необратима.
Итак, для вашего кода я бы сделал что-то вроде этого:
int main(){
double M[3][3] = { {1 , 2 , 3},
{4 , 5 , 6},
{7 , 8 , 9}}
double pivotArray[3]; //since our matrix has three rows
int errorHandler;
double lapackWorkspace[9];
// dgetrf(M,N,A,LDA,IPIV,INFO) means invert LDA columns of an M by N matrix
// called A, sending the pivot indices to IPIV, and spitting error
// information to INFO.
// also don't forget (like I did) that when you pass a two-dimensional array
// to a function you need to specify the number of "rows"
dgetrf_(3,3,M[3][],3,pivotArray[3],&errorHandler);
//some sort of error check
dgetri_(3,M[3][],3,pivotArray[3],9,lapackWorkspace,&errorHandler);
//another error check
}
Вот рабочий код для вычисления обратной матрицы с помощью lapack на C/C++:
#include <cstdio>
extern "C" {
// LU decomoposition of a general matrix
void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);
// generate inverse of a matrix given its LU decomposition
void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO);
}
void inverse(double* A, int N)
{
int *IPIV = new int[N+1];
int LWORK = N*N;
double *WORK = new double[LWORK];
int INFO;
dgetrf_(&N,&N,A,&N,IPIV,&INFO);
dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO);
delete IPIV;
delete WORK;
}
int main(){
double A [2*2] = {
1,2,
3,4
};
inverse(A, 2);
printf("%f %f\n", A[0], A[1]);
printf("%f %f\n", A[2], A[3]);
return 0;
}