我希望能夠使用lapack在C/C++中計算一般的NxN矩陣的逆矩陣。使用lapack在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;
}
你將如何完成它使用dgetri_獲得3x3矩陣M的逆？
首先，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的 int，並且WORK應該是 帶有LWORK元素的雙數組。
INFO只是一個狀態變量，以 告訴你操作 是否成功完成。由於不是所有的 矩陣都是可逆的，所以我會推薦你把這個發送給一些 類型的錯誤檢查系統。成功操作時INFO = 0，如果第i個參數的輸入值不正確，則INFO = -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
}
下面是在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;
}
下面是Spencer Nelson上面的示例的工作版本。關於它的一個謎就是輸入矩陣按照行優先順序排列，即使它看起來叫做底層Fortran例程dgetri。我被引導認爲所有Fortran的例程都需要列主要的順序，但我並不是LAPACK方面的專家，事實上，我正在使用這個例子來幫助我學習它。但是，拋開那個謎團：
示例中的輸入矩陣是單數。 LAPACK試圖告訴你，通過在errorHandler中返回一個3。我將該矩陣中的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
我做的LAPACKE庫中nm，發現如下：
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 [原文如此]包裝，可能會緩解我們不得不爲了fortran的方便而到處尋址，但我可能不會試着去嘗試它，因爲我有前進的方向。
EDIT
這裏是繞過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; }
構建並運行這樣的：
$ 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
編輯
謎似乎不再是一個謎。我認爲計算是按照列主要的順序完成的，因爲它們必須這樣做，但我同時輸入和打印這些矩陣，就好像它們是按照主要順序排列的一樣。我有兩個互相抵消的錯誤，所以事情看起來很划算，儘管它們是列式的。
以上是使用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
關於1X9或3×3。內存佈局方面確實沒有任何區別。事實上，BLAS/LAPACK例程不需要2d陣列。他們採用一維數組，並假設你如何佈置它。請注意BLAS和LAPACK將假設FORTRAN排序（行更改最快）而不是C排序。 – MRocklin 2012-10-18 18:28:09
您可能需要'LAPACKE_dgetrf'而不是fortran例程'dgetrf_'。同上'dgetri'。否則你必須使用地址來調用所有內容。 – 2016-05-20 02:32:48