Calcular la inversa de una matriz usando lapack en C

Tenemos la solución a este enigma, al menos eso esperamos. Si sigues con interrogantes dínoslo, que con gusto te ayudaremos

Solución:

Aquí está el código de trabajo para calcular la inversa de una matriz usando lapack en C / C ++:

#include 

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 %fn", A[0], A[1]);
    printf("%f %fn", A[2], A[3]);

    return 0;

Primero, M tiene que ser bidimensional array, igual que double M[3][3]. Tu array es, matemáticamente hablando, un vector 1×9, que no es invertible.

• N es un puntero a un int para el orden de la matriz; en este caso, N = 3.

• A es un puntero a la factorización LU de la matriz, que puede obtener ejecutando la rutina LAPACK dgetrf.

• LDA es un número entero para el “elemento principal” de la matriz, lo que le permite seleccionar un subconjunto de una matriz más grande si solo desea invertir una pequeña parte. Si desea invertir toda la matriz, LDA debería ser igual a N.

• IPIV son los índices pivote de la matriz, en otras palabras, es una lista de instrucciones de qué filas intercambiar para invertir la matriz. IPIV debe ser generado por la rutina LAPACK dgetrf.

• LWORK y WORK son los “espacios de trabajo” utilizados por LAPACK. Si está invirtiendo toda la matriz, LWORK debe ser un int igual a N ^ 2, y WORK debe ser un doble array con elementos LWORK.

• INFO es solo una variable de estado para indicarle si la operación se completó correctamente. Dado que no todas las matrices son invertibles, le recomiendo que envíe esto a algún tipo de sistema de verificación de errores. INFO = 0 para una operación exitosa, INFO = -i si el i-ésimo argumento tenía un valor de entrada incorrecto e INFO> 0 si la matriz no es invertible.

Entonces, para su código, haría algo como esto:

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

    

Aquí hay una versión funcional de lo anterior usando la interfaz OpenBlas para LAPACKE. Enlace con biblioteca openblas (LAPACKE ya está incluido)

#include 
#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');

Ejemplo:

% 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 

Al final de la página puedes encontrar las explicaciones de otros creadores, tú igualmente puedes dejar el tuyo si lo deseas.

---

---

Utiliza Nuestro Buscador