Inverse of a matrix

Inverse of a matrix

[Viadrina]

Hi, I would like to compute a inverse of a square matrix A. Is there any direct method implemented in gsl for 
computing inverses of suqare matrices. I want to compute a quadratic form
		-1
	Q = X'A  X

Could anybody help me with this topic. I'm not very familiar with gsl.

Przem

Re: Inverse of a matrix

[Alan Aspuru-Guzik]

Hi,
What I do in my code is:

typedef struct Determinant {
    gsl_matrix * matrix; //!original matrix
    gsl_matrix * inverse; //!inverse matrix
    gsl_matrix * ludecomp; //!LU decomposition
    int size; //!number of electrons
    SPIN_TYPE spin; //!spin type of the determinant
    double det; //! determinant
    gsl_permutation * perm; //! permutation matrix
} Determinant;

Determinant new_determinant(int size, SPIN_TYPE spin) {
    Determinant d;

    d.matrix = gsl_matrix_calloc(size,size);
    d.inverse = gsl_matrix_calloc(size,size);
    d.ludecomp = gsl_matrix_calloc(size,size);
    d.size = size;
    d.spin = spin;
    d.det = 0;
    d.perm = gsl_permutation_alloc(size);
    return d;
}


gboolean set_determinant_via_LU_decomp(Determinant * D,
                                       gboolean CALC_INV) {
    int s;

    gsl_matrix_memcpy(D->ludecomp,D->matrix);
    gsl_linalg_LU_decomp(D->ludecomp,D->perm,&s);
        D->det =gsl_linalg_LU_det(D->ludecomp,s);
    if(CALC_INV) {
        //g_message("Inverting matrix!");
        gsl_linalg_LU_invert(D->ludecomp,D->perm,D->inverse);

    }

    return TRUE;
}

Greetings,
Alan

On Tue, 28 May 2002, Viadrina wrote:

> Hi, I would like to compute a inverse of a square matrix A. Is there
> any direct method implemented in gsl for computing inverses of suqare
> matrices. I want to compute a quadratic form
> 		-1
> 	Q = X'A  X
> 
> Could anybody help me with this topic. I'm not very familiar with gsl.
> 
> Przem
> 
> 
> 

-- 

Alan Aspuru-Guzik                    Dios mueve al jugador, y éste, la pieza.
(510)642-5911 UC Berkeley           ¿Qué Dios detrás de Dios la trama empieza
(925)422-8739 LLNL                de polvo y tiempo y sueño y agonías? -Borges

Re: Inverse of a matrix

[Alan Aspuru-Guzik]

Hi,
What I do in my code is:

typedef struct Determinant {
    gsl_matrix * matrix; //!original matrix
    gsl_matrix * inverse; //!inverse matrix
    gsl_matrix * ludecomp; //!LU decomposition
    int size; //!number of electrons
    SPIN_TYPE spin; //!spin type of the determinant
    double det; //! determinant
    gsl_permutation * perm; //! permutation matrix
} Determinant;

Determinant new_determinant(int size, SPIN_TYPE spin) {
    Determinant d;

    d.matrix = gsl_matrix_calloc(size,size);
    d.inverse = gsl_matrix_calloc(size,size);
    d.ludecomp = gsl_matrix_calloc(size,size);
    d.size = size;
    d.spin = spin;
    d.det = 0;
    d.perm = gsl_permutation_alloc(size);
    return d;
}


gboolean set_determinant_via_LU_decomp(Determinant * D,
                                       gboolean CALC_INV) {
    int s;

    gsl_matrix_memcpy(D->ludecomp,D->matrix);
    gsl_linalg_LU_decomp(D->ludecomp,D->perm,&s);
        D->det =gsl_linalg_LU_det(D->ludecomp,s);
    if(CALC_INV) {
        //g_message("Inverting matrix!");
        gsl_linalg_LU_invert(D->ludecomp,D->perm,D->inverse);

    }

    return TRUE;
}

Greetings,
Alan

On Tue, 28 May 2002, Viadrina wrote:

> Hi, I would like to compute a inverse of a square matrix A. Is there
> any direct method implemented in gsl for computing inverses of suqare
> matrices. I want to compute a quadratic form
> 		-1
> 	Q = X'A  X
> 
> Could anybody help me with this topic. I'm not very familiar with gsl.
> 
> Przem
> 
> 
> 

-- 

Alan Aspuru-Guzik                    Dios mueve al jugador, y éste, la pieza.
(510)642-5911 UC Berkeley           ¿Qué Dios detrás de Dios la trama empieza
(925)422-8739 LLNL                de polvo y tiempo y sueño y agonías? -Borges

Inverse of a matrix

[Viadrina]

Hi, I would like to compute a inverse of a square matrix A. Is there any direct method implemented in gsl for 
computing inverses of suqare matrices. I want to compute a quadratic form
		-1
	Q = X'A  X

Could anybody help me with this topic. I'm not very familiar with gsl.

Przem

Inverse of a matrix

[Viadukt]

[-- Attachment #1: Type: text/plain, Size: 163 bytes --]

Hi,
					-1				     -1	
I computed an inverse of a matrix A  but after the multiplication AA  I didn't get the identity.

Could somebody help me, please.

Thans



[-- Attachment #2: inverse.c --]
[-- Type: application/octet-stream, Size: 1379 bytes --]

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

#include<gsl/gsl_matrix.h>
#include<gsl/gsl_blas.h>
#include<gsl/gsl_linalg.h>
#include<gsl/gsl_permutation.h>

#define FIRST_DIMENSION 4
#define SECOND_DIMENSION 4

int main()

	{
	int signum;
	size_t i, j;
	gsl_matrix *A, *B, *C, *D;
	clock_t current_time, last_time;
	
	gsl_permutation *p;

	last_time = clock();

	A  = gsl_matrix_calloc(FIRST_DIMENSION,  SECOND_DIMENSION);
    B  = gsl_matrix_calloc(SECOND_DIMENSION,  FIRST_DIMENSION);
    C  = gsl_matrix_calloc(FIRST_DIMENSION,  FIRST_DIMENSION);

	p = gsl_permutation_alloc(FIRST_DIMENSION);

	for(i=0; i<FIRST_DIMENSION; i++)
	{
		for(j=0; j<SECOND_DIMENSION; j++)
		{
		if(i == j)
			gsl_matrix_set(A, i, j, 1);
		else gsl_matrix_set(A, i, j, 0.3);
        gsl_matrix_set(B, i, j, 2);
		}
	}

	gsl_linalg_LU_decomp(A, p, &signum);
	gsl_linalg_LU_invert(A, p, B);

	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, A, B, 0, C);

	gsl_linalg_cholesky_decomp(A);

	for(i=0; i<FIRST_DIMENSION; i++)
		{
		for(j=0; j<FIRST_DIMENSION; j++)
			{
			printf("%.15lf  ", gsl_matrix_get(C, i, j));
			}
		printf("\n");
		}

	current_time = clock();
	printf("Time elapsed:\t%d\n", (current_time - last_time));

	gsl_matrix_free(A);
	gsl_matrix_free(B);
    gsl_matrix_free(C);

	return 0;

	}

Re: Inverse of a matrix

[Fleur Kelpin]

Hi,

From your code:

	gsl_linalg_LU_decomp(A, p, &signum);
	gsl_linalg_LU_invert(A, p, B);

	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, A, B, 0, C);

From the manual for gsl_linalg_LU_decomp:

On output the diagonal and upper triangular part of the input matrix A
contain the matrix U. The lower triangular part of the input matrix
(excluding the diagonal) contains L. The diagonal elements of L are unity,
and are not stored.

So the LU decomposition of A changes the contents of A.
If you wish to check the results, you will have to make a copy of A before
the decomposition.

Hope this helps,
Fleur

> Hi,
> 					-1				     -1
> I computed an inverse of a matrix A  but after the multiplication AA  I didn't get the identity.
>
> Could somebody help me, please.
>
> Thans
>
>
>

Re: Inverse of a matrix

[Fleur Kelpin]

Hi,

From your code:

	gsl_linalg_LU_decomp(A, p, &signum);
	gsl_linalg_LU_invert(A, p, B);

	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, A, B, 0, C);

From the manual for gsl_linalg_LU_decomp:

On output the diagonal and upper triangular part of the input matrix A
contain the matrix U. The lower triangular part of the input matrix
(excluding the diagonal) contains L. The diagonal elements of L are unity,
and are not stored.

So the LU decomposition of A changes the contents of A.
If you wish to check the results, you will have to make a copy of A before
the decomposition.

Hope this helps,
Fleur

> Hi,
> 					-1				     -1
> I computed an inverse of a matrix A  but after the multiplication AA  I didn't get the identity.
>
> Could somebody help me, please.
>
> Thans
>
>
>

Inverse of a matrix

[Viadukt]

[-- Attachment #1: Type: text/plain, Size: 163 bytes --]

Hi,
					-1				     -1	
I computed an inverse of a matrix A  but after the multiplication AA  I didn't get the identity.

Could somebody help me, please.

Thans



[-- Attachment #2: inverse.c --]
[-- Type: application/octet-stream, Size: 1379 bytes --]

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

#include<gsl/gsl_matrix.h>
#include<gsl/gsl_blas.h>
#include<gsl/gsl_linalg.h>
#include<gsl/gsl_permutation.h>

#define FIRST_DIMENSION 4
#define SECOND_DIMENSION 4

int main()

	{
	int signum;
	size_t i, j;
	gsl_matrix *A, *B, *C, *D;
	clock_t current_time, last_time;
	
	gsl_permutation *p;

	last_time = clock();

	A  = gsl_matrix_calloc(FIRST_DIMENSION,  SECOND_DIMENSION);
    B  = gsl_matrix_calloc(SECOND_DIMENSION,  FIRST_DIMENSION);
    C  = gsl_matrix_calloc(FIRST_DIMENSION,  FIRST_DIMENSION);

	p = gsl_permutation_alloc(FIRST_DIMENSION);

	for(i=0; i<FIRST_DIMENSION; i++)
	{
		for(j=0; j<SECOND_DIMENSION; j++)
		{
		if(i == j)
			gsl_matrix_set(A, i, j, 1);
		else gsl_matrix_set(A, i, j, 0.3);
        gsl_matrix_set(B, i, j, 2);
		}
	}

	gsl_linalg_LU_decomp(A, p, &signum);
	gsl_linalg_LU_invert(A, p, B);

	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, A, B, 0, C);

	gsl_linalg_cholesky_decomp(A);

	for(i=0; i<FIRST_DIMENSION; i++)
		{
		for(j=0; j<FIRST_DIMENSION; j++)
			{
			printf("%.15lf  ", gsl_matrix_get(C, i, j));
			}
		printf("\n");
		}

	current_time = clock();
	printf("Time elapsed:\t%d\n", (current_time - last_time));

	gsl_matrix_free(A);
	gsl_matrix_free(B);
    gsl_matrix_free(C);

	return 0;

	}