Ich möchte gerne die Inverse eines Generals berechnen können NxN Matrix in C/C++ mit Lapack.
Mein Verständnis ist, dass der Weg, eine Inversion in Lapack durchzuführen, die Verwendung von ist dgetri Funktion, aber ich kann nicht herausfinden, was all ihre Argumente sein sollen.
Hier ist der Code, den ich habe:
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;
}
Wie würden Sie es vervollständigen, um die Umkehrung von zu erhalten? 3x3 Matrix M mit dgetri_?
Hier ist der Arbeitscode zum Berechnen der Umkehrung einer Matrix mit Lapack in 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];
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;
}
Erstens muss M ein zweidimensionales Array sein, wie z double M[3][3]. Ihr Array ist mathematisch gesehen ein 1×9-Vektor, der nicht invertierbar ist.
-
N ist ein Zeiger auf ein int für die Ordnung der Matrix – in diesem Fall N=3.
-
A ist ein Zeiger auf die LU-Faktorisierung der Matrix, die Sie durch Ausführen der LAPACK-Routine erhalten können dgetrf.
-
LDA ist eine Ganzzahl für das “führende Element” der Matrix, mit der Sie eine Teilmenge einer größeren Matrix auswählen können, wenn Sie nur ein kleines Stück umkehren möchten. Wenn Sie die gesamte Matrix invertieren möchten, sollte LDA nur gleich N sein.
-
IPIV sind die Pivot-Indizes der Matrix, mit anderen Worten, es ist eine Liste von Anweisungen, welche Zeilen ausgetauscht werden müssen, um die Matrix umzukehren. IPIV sollte von der LAPACK-Routine generiert werden dgetrf.
-
LWORK und WORK sind die von LAPACK verwendeten “Arbeitsbereiche”. Wenn Sie die gesamte Matrix invertieren, sollte LWORK ein int gleich N^2 sein, und WORK sollte ein doppeltes Array mit LWORK-Elementen sein.
-
INFO ist nur eine Statusvariable, die Ihnen mitteilt, ob die Operation erfolgreich abgeschlossen wurde. Da nicht alle Matrizen invertierbar sind, würde ich empfehlen, dass Sie dies an eine Art Fehlerprüfsystem senden. INFO=0 für erfolgreiche Operation, INFO=-i wenn das i-te Argument einen falschen Eingabewert hatte und INFO > 0 wenn die Matrix nicht invertierbar ist.
Also, für Ihren Code würde ich so etwas tun:
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
}
Hier ist eine funktionierende Version des oben Genannten, die die OpenBlas-Schnittstelle zu LAPACKE verwendet. Verknüpfung mit openblas-Bibliothek (LAPACKE ist bereits enthalten)
#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');
}
Beispiel:
% 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
Hier ist eine funktionierende Version des obigen Beispiels von Spencer Nelson. Ein Rätsel dabei ist, dass die Eingabematrix in Zeilenhauptordnung ist, obwohl sie anscheinend die zugrunde liegende Fortran-Routine aufruft dgetri. Ich werde zu der Annahme verleitet, dass alle zugrunde liegenden Fortran-Routinen eine Spalten-Major-Ordnung erfordern, aber ich bin kein Experte für LAPACK, tatsächlich verwende ich dieses Beispiel, um mir dabei zu helfen, es zu lernen. Aber abgesehen von diesem einen Rätsel:
Die Eingabematrix im Beispiel ist singulär. LAPACK versucht Ihnen dies mitzuteilen, indem es a zurückgibt 3 in dem errorHandler. Ich habe die geändert 9 in dieser Matrix zu a 19immer ein errorHandler von 0 Erfolg signalisiert und das Ergebnis mit dem von verglichen Mathematica. Der Vergleich war ebenfalls erfolgreich und bestätigte, dass die Matrix im Beispiel wie dargestellt in Zeilenhauptordnung sein sollte.
Hier ist der Arbeitscode:
#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; }
Ich habe es auf einem Mac wie folgt erstellt und ausgeführt:
gcc main.c -llapacke -llapack
./a.out
Ich habe ein nm in der LAPACKE-Bibliothek und habe Folgendes gefunden:
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
und es sieht so aus, als gäbe es einen LAPACKE [sic] Wrapper, der uns vermutlich davon befreien würde, Adressen überall hin mitnehmen zu müssen, um es Fortran bequemer zu machen, aber ich werde wahrscheinlich nicht dazu kommen, es zu versuchen, weil ich einen Weg nach vorne habe.
BEARBEITEN
Hier ist eine funktionierende Version, die LAPACKE umgeht [sic], unter direkter Verwendung von LAPACK Fortran-Routinen. Ich verstehe nicht, warum eine Row-Major-Eingabe korrekte Ergebnisse liefert, aber ich habe es erneut in Mathematica bestätigt.
#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; }
gebaut und läuft so:
$ 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
BEARBEITEN
Das Mysterium scheint kein Mysterium mehr zu sein. Ich denke, die Berechnungen werden in der Hauptspaltenreihenfolge durchgeführt, wie es sein muss, aber ich gebe die Matrizen sowohl ein als auch drucke ich sie, als ob sie in der Zeilenhauptreihenfolge wären. Ich habe zwei Fehler, die sich gegenseitig aufheben, sodass die Dinge zeilenhaft aussehen, obwohl sie spaltenhaft sind.
