C语言 在 C 中使用 lapack 计算矩阵的逆
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/3519959/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Computing the inverse of a matrix using lapack in C
提问by dzhelil
I would like to be able to compute the inverse of a general NxNmatrix in C/C++ using lapack.
我希望能够NxN使用 lapack计算C/C++ 中一般矩阵的逆。
My understanding is that the way to do an inversion in lapack is by using the dgetrifunction, however, I can't figure out what all of its arguments are supposed to be.
我的理解是在 lapack 中进行反演的方法是使用该dgetri函数，但是，我无法弄清楚它的所有参数应该是什么。
Here is the code I have:
这是我的代码：
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;
}
How would you complete it to obtain the inverse of the 3x3matrix M using dgetri_?
您将如何3x3使用 dgetri_完成它以获取矩阵 M的逆矩阵？
采纳答案by spencer nelson
First, M has to be a two-dimensional array, like double M[3][3]. Your array is, mathematically speaking, a 1x9 vector, which is not invertible.
首先，M 必须是一个二维数组，比如double M[3][3]。从数学上讲，您的数组是一个不可逆的 1x9 向量。
N is a pointer to an int for the order of the matrix - in this case, N=3.
A is a pointer to the LU factorization of the matrix, which you can get by running the LAPACK routine
dgetrf.LDA is an integer for the "leading element" of the matrix, which lets you pick out a subset of a bigger matrix if you want to just invert a little piece. If you want to invert the whole matrix, LDA should just be equal to N.
IPIV is the pivot indices of the matrix, in other words, it's a list of instructions of what rows to swap in order to invert the matrix. IPIV should be generated by the LAPACK routine
dgetrf.LWORK and WORK are the "workspaces" used by LAPACK. If you are inverting the whole matrix, LWORK should be an int equal to N^2, and WORK should be a double array with LWORK elements.
INFO is just a status variable to tell you whether the operation completed successfully. Since not all matrices are invertible, I would recommend that you send this to some sort of error-checking system. INFO=0 for successful operation, INFO=-i if the i'th argument had an incorrect input value, and INFO > 0 if the matrix is not invertible.
N 是指向矩阵顺序的 int 的指针 - 在这种情况下，N = 3。
A 是指向矩阵的 LU 分解的指针，您可以通过运行 LAPACK 例程来获得它
dgetrf。LDA 是矩阵“前导元素”的整数，如果您只想反转一小部分，它可以让您挑选出更大矩阵的子集。如果要反转整个矩阵，LDA 应该等于 N。
IPIV 是矩阵的主元索引，换句话说，它是交换行以反转矩阵的指令列表。IPIV 应该由 LAPACK 例程生成
dgetrf。LWORK 和 WORK 是 LAPACK 使用的“工作区”。如果要反转整个矩阵，LWORK 应该是一个等于 N^2 的 int，而 WORK 应该是一个包含 LWORK 元素的双数组。
INFO 只是一个状态变量，用于告诉您操作是否成功完成。由于并非所有矩阵都是可逆的，我建议您将其发送到某种错误检查系统。INFO=0 表示操作成功，INFO=-i 如果第 i 个参数的输入值不正确，INFO > 0 如果矩阵不可逆。
So, for your code, I would do something like this:
所以，对于你的代码，我会做这样的事情：
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
}
回答by dzhelil
Here is the working code for computing the inverse of a matrix using lapack in C/C++:
这是在 C/C++ 中使用 lapack 计算矩阵逆的工作代码：
#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;
}
回答by PolarBear2015
Here is a working version of the above using OpenBlas interface to LAPACKE. Link with openblas library (LAPACKE is already contained)
这是使用 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');
}
Example:
例子：
% 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
回答by Reb.Cabin
Here is a working version of Spencer Nelson's example above. One mystery about it is that the input matrix is in row-major order, even though it appears to call the underlying fortran routine dgetri. I am led to believe that all the underlying fortran routines require column-major order, but I am no expert on LAPACK, in fact, I'm using this example to help me learn it. But, that one mystery aside:
这是上面 Spencer Nelson 示例的工作版本。关于它的一个谜是输入矩阵是按行优先顺序排列的，即使它似乎调用了底层的 fortran 例程dgetri。我相信所有底层的 fortran 例程都需要列主序，但我不是 LAPACK 的专家，事实上，我正在使用这个例子来帮助我学习它。但是，撇开一个谜：
The input matrix in the example is singular. LAPACK tries to tell you that by returning a 3in the errorHandler. I changed the 9in that matrix to a 19, getting an errorHandlerof 0signalling success, and compared the result to that from Mathematica. The comparison was also successful and confirmed that the matrix in the example should be in row-major order, as presented.
示例中的输入矩阵是奇异矩阵。LAPACK 试图通过3在errorHandler. 我9将该矩阵中的更改为 a 19，获得了成功errorHandler的0信号，并将结果与来自 的结果进行了比较Mathematica。比较也成功并确认示例中的矩阵应按行优先顺序排列，如所示。
Here is the working code:
这是工作代码：
#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; }
I built and ran it as follows on a Mac:
我在 Mac 上按如下方式构建和运行它：
gcc main.c -llapacke -llapack
./a.out
I did an nmon the LAPACKE library and found the following:
我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
and it looks like there is a LAPACKE [sic] wrapper that would presumably relieve us of having to take addresses everywhere for fortran's convenience, but I am probably not going to get around to trying it because I have a way forward.
看起来有一个 LAPACKE [原文如此] 包装器可以让我们不必为了 fortran 的方便而到处获取地址，但我可能不会尝试尝试它，因为我有前进的道路。
EDIT
编辑
Here is a working version that bypasses LAPACKE [sic], using LAPACK fortran routines directly. I do not understand why a row-major input produces correct results, but I confirmed it again in Mathematica.
这是一个绕过 LAPACKE [原文如此] 的工作版本，直接使用 LAPACK fortran 例程。我不明白为什么行主要输入会产生正确的结果，但我在 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; }
built and run like this:
像这样构建和运行：
$ 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
EDIT
编辑
The mystery no longer appears to be a mystery. I think the computations are being done in column-major order, as they must, but I am both inputting and printing the matrices as if they were in row-major order. I have two bugs that cancel each other out so things look row-ish even though they're column-ish.
这个谜似乎不再是谜。我认为计算是按列优先顺序进行的，因为它们必须这样做，但是我输入和打印矩阵都好像它们是按行优先顺序一样。我有两个相互抵消的错误，所以即使它们是列式的，事情看起来也是行式的。
