Вычисление обратной матрицы с использованием Lapack в C
Я хотел бы иметь возможность вычислить обратное к общему NxN Матрица в C/C++ с использованием Lapack.
Насколько я понимаю, что способ сделать инверсию в 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;
}
Как бы вы завершили это, чтобы получить обратное 3x3 матрица М с помощью dgetri_?
4 ответа
Во-первых, 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 должен быть целым числом, равным N^2, и WORK должен быть двойным массивом с элементами 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;
}
Вот рабочая версия выше, используя интерфейс OpenBlas для LAPACKE. Связь с библиотекой openblas (LAPACKE уже содержится)
#include <stdio.h>
#include "cblas.h"
#include "lapacke.h"
// inplace inverse n x n matrix A.
// matrix A is Column Major (i.e. firts line, second line ... *not* C[][] order)
// returns:
// ret = 0 on success
// ret < 0 illegal argument value
// ret > 0 singular matrix
lapack_int matInv(double *A, unsigned n)
{
int ipiv[n+1];
lapack_int ret;
ret = LAPACKE_dgetrf(LAPACK_COL_MAJOR,
n,
n,
A,
n,
ipiv);
if (ret !=0)
return ret;
ret = LAPACKE_dgetri(LAPACK_COL_MAJOR,
n,
A,
n,
ipiv);
return ret;
}
int main()
{
double A[] = {
0.378589, 0.971711, 0.016087, 0.037668, 0.312398,
0.756377, 0.345708, 0.922947, 0.846671, 0.856103,
0.732510, 0.108942, 0.476969, 0.398254, 0.507045,
0.162608, 0.227770, 0.533074, 0.807075, 0.180335,
0.517006, 0.315992, 0.914848, 0.460825, 0.731980
};
for (int i=0; i<25; i++) {
if ((i%5) == 0) putchar('\n');
printf("%+12.8f ",A[i]);
}
putchar('\n');
matInv(A,5);
for (int i=0; i<25; i++) {
if ((i%5) == 0) putchar('\n');
printf("%+12.8f ",A[i]);
}
putchar('\n');
}
Пример:
% g++ -I [OpenBlas path]/include/ example.cpp [OpenBlas path]/lib/libopenblas.a
% a.out
+0.37858900 +0.97171100 +0.01608700 +0.03766800 +0.31239800
+0.75637700 +0.34570800 +0.92294700 +0.84667100 +0.85610300
+0.73251000 +0.10894200 +0.47696900 +0.39825400 +0.50704500
+0.16260800 +0.22777000 +0.53307400 +0.80707500 +0.18033500
+0.51700600 +0.31599200 +0.91484800 +0.46082500 +0.73198000
+0.24335255 -2.67946180 +3.57538817 +0.83711880 +0.34704217
+1.02790497 -1.05086895 -0.07468137 +0.71041070 +0.66708313
-0.21087237 -4.47765165 +1.73958308 +1.73999641 +3.69324020
-0.14100897 +2.34977565 -0.93725915 +0.47383541 -2.15554470
-0.26329660 +6.46315378 -4.07721533 -3.37094863 -2.42580445
Вот рабочая версия примера Спенсера Нельсона выше. Одна загадка в этом заключается в том, что матрица ввода находится в главном порядке строк, хотя она, кажется, вызывает основную процедуру на Фортране. dgetri, Я склонен полагать, что все базовые процедуры фортрана требуют порядка главных колонн, но я не эксперт по LAPACK, на самом деле, я использую этот пример, чтобы помочь мне изучить его. Но это одна загадка в стороне:
Входная матрица в примере является единственной. LAPACK пытается сказать вам это, возвращая 3 в errorHandler, Я изменил 9 в этой матрице к 19, получив errorHandler из 0 сигнал успеха и сравнил результат с Mathematica, Сравнение также было успешным и подтвердило, что матрица в примере должна быть в главном порядке строк, как представлено.
Вот рабочий код:
#include <stdio.h>
#include <stddef.h>
#include <lapacke.h>
int main() {
int N = 3;
int NN = 9;
double M[3][3] = { {1 , 2 , 3},
{4 , 5 , 6},
{7 , 8 , 9} };
int 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_(&N, &N, M[0], &N, pivotArray, &errorHandler);
printf ("dgetrf eh, %d, should be zero\n", errorHandler);
dgetri_(&N, M[0], &N, pivotArray, lapackWorkspace, &NN, &errorHandler);
printf ("dgetri eh, %d, should be zero\n", errorHandler);
for (size_t row = 0; row < N; ++row)
{ for (size_t col = 0; col < N; ++col)
{ printf ("%g", M[row][col]);
if (N-1 != col)
{ printf (", "); } }
if (N-1 != row)
{ printf ("\n"); } }
return 0; }
Я построил и запустил его на Mac следующим образом:
gcc main.c -llapacke -llapack
./a.out
Я сделал nm в библиотеке LAPACKE и обнаружил следующее:
liblapacke.a(lapacke_dgetri.o):
U _LAPACKE_dge_nancheck
0000000000000000 T _LAPACKE_dgetri
U _LAPACKE_dgetri_work
U _LAPACKE_xerbla
U _free
U _malloc
liblapacke.a(lapacke_dgetri_work.o):
U _LAPACKE_dge_trans
0000000000000000 T _LAPACKE_dgetri_work
U _LAPACKE_xerbla
U _dgetri_
U _free
U _malloc
и похоже, что существует оболочка LAPACKE [sic], которая, вероятно, избавит нас от необходимости брать адреса везде для удобства Фортрана, но я, вероятно, не собираюсь пытаться это сделать, потому что у меня есть путь вперед.
РЕДАКТИРОВАТЬ
Вот рабочая версия, которая обходит LAPACKE [sic], напрямую используя подпрограммы LAPACK. Я не понимаю, почему ввод основной строки дает правильные результаты, но я подтвердил это снова в Mathematica.
#include <stdio.h>
#include <stddef.h>
int main() {
int N = 3;
int NN = 9;
double M[3][3] = { {1 , 2 , 3},
{4 , 5 , 6},
{7 , 8 , 19} };
int pivotArray[3]; //since our matrix has three rows
int errorHandler;
double lapackWorkspace[9];
/* from http://www.netlib.no/netlib/lapack/double/dgetrf.f
SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
*
* -- LAPACK routine (version 3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
* ..
* .. Array Arguments ..
INTEGER IPIV( * )
DOUBLE PRECISION A( LDA, * )
*/
extern void dgetrf_ (int * m, int * n, double * A, int * LDA, int * IPIV,
int * INFO);
/* from http://www.netlib.no/netlib/lapack/double/dgetri.f
SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
*
* -- LAPACK routine (version 3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LWORK, N
* ..
* .. Array Arguments ..
INTEGER IPIV( * )
DOUBLE PRECISION A( LDA, * ), WORK( * )
*/
extern void dgetri_ (int * n, double * A, int * LDA, int * IPIV,
double * WORK, int * LWORK, int * INFO);
// 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_(&N, &N, M[0], &N, pivotArray, &errorHandler);
printf ("dgetrf eh, %d, should be zero\n", errorHandler);
dgetri_(&N, M[0], &N, pivotArray, lapackWorkspace, &NN, &errorHandler);
printf ("dgetri eh, %d, should be zero\n", errorHandler);
for (size_t row = 0; row < N; ++row)
{ for (size_t col = 0; col < N; ++col)
{ printf ("%g", M[row][col]);
if (N-1 != col)
{ printf (", "); } }
if (N-1 != row)
{ printf ("\n"); } }
return 0; }
построен и работает так:
$ gcc foo.c -llapack
$ ./a.out
dgetrf eh, 0, should be zero
dgetri eh, 0, should be zero
-1.56667, 0.466667, 0.1
1.13333, 0.0666667, -0.2
0.1, -0.2, 0.1
РЕДАКТИРОВАТЬ
Тайна больше не кажется загадкой. Я думаю, что вычисления выполняются в основном порядке столбцов, как они должны, но я и ввод, и печать матриц, как если бы они были в основном порядке строк. У меня есть две ошибки, которые взаимно компенсируют друг друга, поэтому вещи выглядят как строки, даже если они столбцы.