Inversa y determinante de una matriz con Lapackpp

La inversa y el determinante de una matriz se pueden obtener con las librerías de álgebra lineal de Lapackpp usando la factorización LU, que combina técnicas de pivoteo que permiten evitar la división por cero en pasos intermedios de cálculo. La matriz de la factorización LU (inferior: Low; superior: Up) puede ser separada en dos matrices triangulares; donde una de ellas (la matriz L) tiene todos sus elementos en la diagonal iguales a 1 y, por tanto, la matriz U puede usarse para la estimación del determinante (en realidad, la función miembro de la Clase en Lapackpp produce una matriz LU integrada en una sóla por lo que el cálculo del determinante se hace directamente en esta matriz).

El cálculo del determinante se obtiene primero que la inversa, a partir de LUFactorizeIP, porque la matriz factorizada original requerida para el cálculo del determinante y que también será usada por LaLUInverseIP (conjuntamente con sus pivotes), será destruida en este último paso (para la estimación correcta del determinante hay que multiplicarlo por -1 elevado al número de filas o de columnas de la matriz más 1 ).

En el ejemplo, vamos a usar dos matrices: una matriz simétrica de Hilbert [aij = 1/(i + j)] y otra donde se truncan los coeficientes de la primera en su tercera cifra decimal. Ambas son pobremente condicionadas. La condición de una matriz se estima a partir del producto de la norma de ésta y de su inversa. Si la condición es cercana a 1 la matriz está bien condicionada mientras que si ésta es >> 1 (> 100 para una matriz pequeña) se dice que está mal condicionada. En Lapackpp existen funciones para evaluar diferentes tipos de norma. Aquí vamos a emplear la norma p=1(Blas_Norm1).

El código es el siguiente:

#include <iostream>
#include <stdio.h>
#include <lapackpp/laslv.h>
#include <lapackpp/blas3pp.h>
#include <cmath>

using namespace std;

int main(){

double a[]={ 1/2.0,  1/3.0,  1/4.0,
             1/3.0,  1/4.0,  1/5.0,
             1/4.0,  1/5.0,  1/6.0};

double b[]={ 0.500,  0.333,  0.250,
             0.333,  0.250,  0.200,
             0.250,  0.200,  0.166};

	LaGenMatDouble A(a,3,3), B(b,3,3); // Declara objetos de la Clase Matriz

	double normaA, normaB, normaAinv, normaBinv;

	normaA = Blas_Norm1(A);

	cout<<"La matriz original A es:\n";
	
	cout << A;

	cout<<"\nLa matriz original B es:\n";
	
	cout << B;

	normaB = Blas_Norm1(B);

	LaVectorLongInt pivots1(A.cols());
	LaVectorLongInt pivots2(B.cols());

	LUFactorizeIP(A, pivots1);
	LUFactorizeIP(B, pivots2);

	double detA = 1, detB = 1;

	for (int i = 0; i < A.cols(); ++i){

	    detA *= A(i, i); detB *= B(i,i);

	}
	
	LaLUInverseIP(A, pivots1);
	LaLUInverseIP(B, pivots2);

	cout << "\nEste es el determinante de A\n";

	cout << pow(-1, A.cols()+1)*detA << "\n";

	cout << "\nEste es el determinante de B\n";

	cout << pow(-1, B.cols()+1)*detB << "\n";
	
	cout << "\nEsta es la matriz inversa de A\n";

	cout << A;

	cout << "\nEsta es la matriz inversa de B\n";

	cout << B;

	cout << "\nnormas de A y de B; respectivamente\n";

	cout << normaA << " " << normaB << "\n";

	normaAinv = Blas_Norm1(A); normaBinv = Blas_Norm1(B);

	cout << "\nnormas de las inversas de A y de B; respectivamente\n";

	cout << normaAinv << " " << normaBinv << "\n";

	cout << "\ncondicion de las matrices A y B; respectivamente\n";

	cout << normaA*normaAinv << " " << normaB*normaBinv << "\n";

	return 0;

}

y ésta es la salida:

La matriz original A es:
0.5  0.333333  0.25
0.333333  0.25  0.2
0.25  0.2  0.166667

La matriz original B es:
0.5  0.333  0.25
0.333  0.25  0.2
0.25  0.2  0.166

Este es el determinante de A
2.31481e-05

Este es el determinante de B
1.7426e-05

Esta es la matriz inversa de A
72  -240  180
-240  900  -720
180  -720  600

Esta es la matriz inversa de B
86.0783  -302.881  235.281
-302.881  1176.4  -961.207
235.281  -961.207  809.767

normas de A y de B; respectivamente
1.08333 1.083

normas de las inversas de A y de B; respectivamente
1860 2440.49

condicion de las matrices A y B; respectivamente
2015 2643.05

Observen la relativa gran diferencia entre ambas inversas y lo pobremente condicionadas que están las dos matrices originales.

Esta entrada fue publicada en Código C++ y etiquetada , . Guarda el enlace permanente.

4 respuestas a Inversa y determinante de una matriz con Lapackpp

  1. Cristobal dijo:

    Muy bien🙂
    Una cosa, no tengo mucho tiempo libre, el otro día intenté hacer lo de los valores propios cuando son complejos y me daba error de librerías, ¿puedes ver si a ti también te da error? Exactamente usando LaGenMatComplex o algo así

    Otra cosa, Low y Up se traducen por Inferior y Superior respectivamente, se debe a lo de matriz triangular inferior y superior😉

  2. OK. Si me puedes mandar el ejemplo que te dió el error te lo agradecería. Así iría al grano.

    “Otra cosa, Low y Up se traducen por Inferior y Superior respectivamente, se debe a lo de matriz triangular inferior y superior ;-)”

    Lo peor es que sé que esa es la traducción y por estar mirando porque no me salía lo del determinante (hasta que lo encontré) lo traduje a la usanza coloquial. Gracias por la aclaratoria. Ya lo corregí.

    Saludos

  3. Cristobal dijo:

    Vale, mira con éste simple código ya me falla http://www.mediafire.com/?vfwdxcklzt3e3ea

    He intentado declarar matrices complejas con complex.h pero no hay forma, no está implementado. Así que he cogido y he hecho lo mismo que con la matriz double de lapackpp pero compleja y me da error.

    Gracias🙂

  4. Tu código así:

    #include <iostream>
    #include <complex>
    #include <lapackpp/laslv.h>
    //#include <lapackpp/gmc.h>
    #include <lapackpp/blas3pp.h>
    
    
    using namespace std;
    
    int main(int argc, char *argv[])
    {
    	double preal=0.,pimag=1.;
    	double preal1=1.,pimag1=-2.;
    	complex <double_t> c1(preal,pimag);
    	complex <double_t> c2(preal1,pimag1);
    	cout<<real(c1)<<" + "<<imag(c1)<<"i"<<endl;
    	cout<<c1+c2<<endl;
    	double A_real[]={ 1.0,  2.0,  3.0,
    				2.0,  1.0,  4.0,
    				3.0,  4.0,  1.0}; 
    	double A_imag[]={ 1.0,  2.0,  3.0,
    				2.0,  1.0,  4.0,
    				3.0,  4.0,  1.0}; 
    	
    	LaGenMatDouble A(A_real,3,3);
    //	LaGenMatComplex A(A_real,A_imag,3,3);
              
    	return 0;
    }
    

    ya compila y se ejecuta. Eso me llevó a ver detalladamente en <lapackpp/gmc.h> y observé lo siguiente:

    00044 #ifndef LA_COMPLEX_SUPPORT
    00045 /* An application must define LA_COMPLEX_SUPPORT if it wants to use
    00046  * complex numbers here. */
    00047 # error "The macro LA_COMPLEX_SUPPORT needs to be defined if you want to use complex-valued matrices."
    00048 #endif
    

    Hay que definir la macro LA_COMPLEX_SUPPORT e incluirlo antes de las librerías de Lapackpp en la aplicación; según pude leer. Ahora hay que saber como.

    Saludos

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s