Re: [Help-gsl] gsl_linalg_LU_invert

[Matthias Kramm]

On Tue, Apr 10, 2007 at 10:20:43AM -0500, Jordi Gutierrez Hermoso wrote:
> On 10/04/07, Matthias Kramm <address@hidden> wrote:
> >the online documentation gives the impression that
> >gsl_linalg_LU_invert() will invert a matrix previously
> >postprocessed with gsl_linalg_LU_decomp().
> 
> Yes, this is what it does. Do you really need to invert a matrix? This
> is hardly ever desirable or needed, you understand.

I believe I do. I need to solve the generalized eigenvector
problem A * x = lambda * B * x with B positive definite, and
am currently using Algorithm 8.7.1 from Golub & Van Loans 
Matrix Computations, which needs to compute an inverse.

> No, this is not true. gsl_linalg_LU_decomp() overwrites the original
> matrix with its LU factorisation, and it saves both the L and the U
> part. For the GSL, L has ones in its diagonal and U has the pivots of
> the original matrix A in its diagonal. The ones of L are not stored,
> of course. The diagonal of the gsl_matrix* A will contain the diagonal
> of U after the call to gsl_linalg_LU_decomp().

Ah! I completely missed that the return matrix of LU_decomp() is
not a triangular matrix. 
Smart way to save space, actually. *bows*

> No, both of the functions you mentioned are working correctly. You can
> output the inverse you computed to satisfy yourself that it does. Your
> bug lies elsewhere, in using gsl_matrix_view_array twice to create
> what you think are two different copies of the same matrix, but in
> fact both refer to the same already modified matrix!

Ok, that would explain the problems I was having.

> Further, you should be using gsl_blas_dgemm instead of
> gsl_linalg_matmult. BLAS functions seem overly complicated, but a lot
> of thought has gone into them. 

Shouldn't gsl_linalg_matmult wrap the BLAS function, if it's so much
better? Might make things easier for new gsl users. 
I mean, it's not exactly obvious from the name that gsl_blas_dgemm()
performs a matrix multiplication.

> I have commented and fixed the bug in your code, which I attach. I

Thanks! Your help is much appreciated.

Greetings

Matthias