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

resolver matrices en c

Question

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

Mi entendimiento es que la manera de hacer una inversión en lapack es mediante el uso de la función dgetri , sin embargo, no puedo entender qué 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 matriz M de 3x3 usando dgetri_?

Answer 1

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

Answer 2

Aquí hay una versión de trabajo de lo anterior usando la interfaz OpenBlas para LAPACKE. Enlace con la librería openblas (LAPACKE ya está contenido)

#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 

Answer 3

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

La matriz de entrada en el ejemplo es singular. LAPACK intenta decirle que devuelve un 3 en el errorHandler . Cambié el 9 en esa matriz a un 19 , obteniendo un errorHandler de 0 éxito de señalización y errorHandler el resultado con Mathematica . La comparación también fue exitosa y confirmó que la matriz en el ejemplo debería estar en el orden de la fila principal, como se presentó.

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 lo 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 envoltorio LAPACKE [sic] que presumiblemente nos libera de tener que tomar direcciones en todas partes para la conveniencia de Fortran, pero probablemente no voy a intentarlo porque tengo una manera de avanzar.

EDITAR

Aquí hay una versión de trabajo que pasa por alto a LAPACKE [sic], usando las rutinas fortan de LAPACK directamente. No entiendo por qué una entrada de fila mayor 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 de esta manera:

$ 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 realizando en orden de columna mayor, como deben ser, pero estoy ingresando e imprimiendo las matrices como si estuvieran en orden de fila mayor. Tengo dos errores que se anulan entre sí, por lo que las cosas se ven como filas, aunque sean de columna.

Answer 4

Primero, M tiene que ser una matriz bidimensional, como double M[3][3] . Su matriz es, matemáticamente hablando, un vector 1x9, 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 el dgetrf rutina LAPACK.

• LDA es un número entero para el "elemento inicial" de la matriz, 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 dinámicos de la matriz, en otras palabras, es una lista de instrucciones sobre qué filas intercambiar para invertir la matriz. IPIV debe ser generado por la rutina dgetrf LAPACK.

• LWORK y WORK son los "espacios de trabajo" utilizados por LAPACK. Si está invirtiendo la matriz completa, LWORK debería ser un int igual a N ^ 2, y WORK debería ser una matriz doble con elementos LWORK.

• INFO es solo una variable de estado para decirle si la operación se completó con éxito. Dado que no todas las matrices son invertibles, recomendaría que lo envíe a algún tipo de sistema de comprobación de errores. INFO = 0 para una operación exitosa, INFO = -i si el argumento i''th 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

    }