Re: [Help-gsl] LU decomposition

[John D Lamb]

On Fri, 2007-03-02 at 09:15 +0200, Keymone wrote:
> question is not how to decompose matrix
> gsl_linalg_LU_decomp stores LU decomposition into matrix A
> and i'm asking how to restore previous matrix without copying it before
> decomposition

I don't think there's an explicit function to do this. You could use
cblas and create your own matrix from the matrix LU created by the
decomposition, but a more efficient approach would contain something
like the following:

size_t n, i, j, k; /* dimension of matrix LU (and A) */
double l_ik, u_kj, a_ij;

...

for( i = 0; i < n; ++i )
    for( j = 0; j < n; ++j ){
        a_ij = 0;
        for( k = 0; k < n; ++k ){
            l_ik = i < k : 0 ? (i==k : 1 ? gsl_matrix_get( LU, i, k ));
            u_kj = j < k : 0 ? gsl_matrix_get( LU, k, j );
            a_ij += l_ij * u_kj;
        }
        gsl_matrix_set( A, i, j, a_ij );
    }

I haven't tested this.

-- 
JDL