Ecco una versione funzionante dell'esempio di Spencer Nelson sopra. Un mistero è che la matrice di input è in ordine di riga principale, anche se sembra chiamare la routine di fortran sottostante dgetri. Sono portato a credere che tutte le routine di fortran sottostanti richiedano un ordine di colonna principale, ma non sono un esperto di LAPACK, infatti, sto usando questo esempio per aiutarmi ad apprenderlo. Ma, a parte questo mistero:
La matrice di input nell'esempio è singolare. LAPACK prova a dirti che restituendo uno 3 nel errorHandler. Ho modificato il 9 in quella matrice in un 19, ottenendo un errorHandler di 0 segnalazione riuscita e confrontato il risultato con quello da Mathematica. Anche il confronto ha avuto esito positivo e ha confermato che la matrice nell'esempio dovrebbe essere in ordine di riga maggiore, come presentato.
Ecco il codice di lavoro:
#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; }
ho costruito e fatto funzionare come segue su un Mac:
gcc main.c -llapacke -llapack
./a.out
ho fatto un nm sulla libreria LAPACKE e trovato il seguente:
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
e sembra che ci sia un wrapper LAPACKE [sic] che presumibilmente allevierà noi di dover prendere indirizzi ovunque per comodità di Fortran, ma probabilmente non ho intenzione di andare in giro a provarlo perché ho una strada da percorrere.
EDIT
Ecco una versione di lavoro che bypassa LAPACKE [sic], utilizzando direttamente le routine FORTRAN LAPACK. Non capisco perché un input di riga principale produca risultati corretti, ma l'ho confermato di nuovo in 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; }
costruite e gestite in questo modo:
$ 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
Il mistero non è più sembra essere un mistero. Penso che i calcoli vengano eseguiti in ordine di colonna principale, come devono, ma sto sia inserendo e stampando le matrici come se fossero in ordine di riga principale. Ho due bug che si annullano a vicenda, quindi le cose sembrano in fila, anche se sono colonne-ish.
Per quanto riguarda 1x9 o 3x3. Non c'è davvero alcuna differenza in termini di layout di memoria. In effetti le routine BLAS/LAPACK non utilizzano array 2D. Prendono 1d array e fanno ipotesi su come lo hai esposto. Attenzione però che BLAS e LAPACK assumeranno l'ordine FORTRAN (le righe cambiano più velocemente) anziché l'ordinamento C. – MRocklin
Potrebbe essere necessario 'LAPACKE_dgetrf' piuttosto che la routine fortran' dgetrf_'. Idem 'dgetri'. Altrimenti devi chiamare tutto con gli indirizzi. –