Cálculo de la inversa de una matriz usando lapack en C

9 minutos de lectura

Calculo de la inversa de una matriz usando lapack en
DR

Me gustaría poder calcular el inverso de un general NxN matriz en C/C++ usando lapack.

Tengo entendido que la forma de hacer una inversión en lapack es usando el dgetri Sin embargo, no puedo descifrar cuáles se supone que son todos sus argumentos.

Aquí está el código que tengo:

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;
}

¿Cómo lo completarías para obtener el inverso de la 3x3 matriz M usando dgetri_?

Calculo de la inversa de una matriz usando lapack en
DR

Aquí está el código de trabajo para calcular el inverso de una matriz usando lapack en 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;
}

  • No necesita asignar N+1 (pero solo N sin firmar) int para la variable IPIV. Además, no recomiendo usar este tipo de función para calcular los múltiplos inversos. Asigne los datos que necesita de una vez por todas y gratis solo al final.

    – matovitch

    22 de noviembre de 2013 a las 11:09

  • Puede que necesites delete[] IPIV y delete [] work.

    – Cabaña Reb

    20 de mayo de 2016 a las 2:13

  • No hay lenguaje C/C++. El código que muestras es C++. Pero la pregunta es sobre C.

    – demasiado honesto para este sitio

    23/06/2016 a las 20:40

  • Hacer que este código C sea válido es tan simple como cambiar el new llama a malloc y el delete[] llamadas a frees (y deshacerse de la “C” externa).

    – brote de alfalfa

    12 de enero de 2017 a las 17:48

Calculo de la inversa de una matriz usando lapack en
spencer nelson

Primero, M tiene que ser una matriz bidimensional, como double M[3][3]. Su matriz es, matemáticamente hablando, un vector de 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, que le permite seleccionar un subconjunto de una matriz más grande si 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 número entero igual a N^2, y WORK debe ser una matriz doble con elementos LWORK.

  • INFO es solo una variable de estado que le indica si la operación se completó con éxito. Dado que no todas las matrices son invertibles, recomendaría 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 tu 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

    }

  • En cuanto a 1×9 o 3×3. Realmente no hay ninguna diferencia en términos de diseño de memoria. De hecho, las rutinas BLAS/LAPACK no toman arreglos 2D. Toman arreglos 1d y hacen suposiciones sobre cómo lo dispuso. Sin embargo, tenga en cuenta que BLAS y LAPACK asumirán el orden FORTRAN (las filas cambian más rápido) en lugar del orden C.

    – M. Rocklin

    18/10/2012 a las 18:28

  • Puede que necesites LAPACKE_dgetrf en lugar de la rutina fortran dgetrf_. Ídem dgetri. De lo contrario, debe llamar a todo con direcciones.

    – Cabaña Reb

    20 de mayo de 2016 a las 2:32

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

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

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 

1647544211 775 Calculo de la inversa de una matriz usando lapack en
Reb.Cabina

Aquí hay una versión funcional del ejemplo anterior de Spencer Nelson. Un misterio al respecto es que la matriz de entrada está en orden de fila principal, aunque parece llamar a la rutina fortran subyacente dgetri. Me hacen creer que todas las rutinas fortran subyacentes requieren un orden de columna principal, pero no soy un experto en LAPACK, de hecho, estoy usando este ejemplo para ayudarme a aprenderlo. Pero, ese misterio aparte:

La matriz de entrada en el ejemplo es singular. LAPACK trata de decirle que al devolver un 3 en el errorHandler. cambié el 9 en esa matriz a un 19obteniendo un errorHandler de 0 señalando el éxito, y comparó el resultado con el de Mathematica. La comparación también fue exitosa y confirmó que la matriz del ejemplo debería estar en el orden de las filas principales, tal como se presenta.

Aquí está el código de trabajo:

#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;   }

Lo construí y ejecuté de la siguiente manera en una Mac:

gcc main.c -llapacke -llapack
./a.out

hice un nm en la biblioteca LAPACKE y encontré lo siguiente:

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

y parece que hay un LAPACKE [sic] envoltorio que presumiblemente nos liberaría de tener que llevar direcciones a todas partes para la conveniencia de fortran, pero probablemente no voy a intentarlo porque tengo un camino a seguir.

EDITAR

Aquí hay una versión de trabajo que pasa por alto LAPACKE [sic], utilizando directamente las rutinas LAPACK fortran. No entiendo por qué una entrada de fila principal produce resultados correctos, pero lo confirmé nuevamente en 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;   }

construido y ejecutado así:

$ 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

EDITAR

El misterio ya no parece ser un misterio. Creo que los cálculos se están haciendo en orden de columna principal, como debe ser, pero estoy ingresando e imprimiendo las matrices como si estuvieran en orden de fila principal. Tengo dos errores que se cancelan entre sí, por lo que las cosas se ven como filas aunque sean como columnas.

¿Ha sido útil esta solución?

Esta web utiliza cookies propias y de terceros para su correcto funcionamiento y para fines analíticos y para mostrarte publicidad relacionada con sus preferencias en base a un perfil elaborado a partir de tus hábitos de navegación. Al hacer clic en el botón Aceptar, acepta el uso de estas tecnologías y el procesamiento de tus datos para estos propósitos. Configurar y más información
Privacidad