Re: [Help-gsl] gsl_linalg_LU_invert

[Jordi Gutierrez Hermoso]

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 don't quite see how that can be possible, as
> gsl_linalg_LU_decomp() will destroy the original matrix
> (transform it into an upper triagonal), and doesn't
> save enough information to "reconstruct" it.
> (only the permutation of the matrix rows is saved, but not
>  the lower triagonal part (U))

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().

> Is this a bug, or am I overlooking something? In any case,
> I can get neither of these functions to work correctly.

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!

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. See Golub & Van Loan for a good
motivation for the need for BLAS. You will observe that that the GSL
no longer mentions gsl_linalg_matmult in its documentation, and I take
it this means its use is to be considered deprecated.

I have commented and fixed the bug in your code, which I attach. I
also "fixed" other things that aren't strictly speaking bugs, but just
stylistic issues I that IMAO could be done better.

HTH,
- Jordi G. H.

gsl_inverse_test.c
Description: Text Data