Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
56 lines (43 sloc) 1.49 KB
/* Easel's interfaces to the GNU Scientific Library
*/
#include "esl_config.h"
#ifdef HAVE_LIBGSL
#include <stdlib.h>
#include "easel/easel.h"
#include "easel/dmatrix.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_eigen.h>
int
esl_GSL_MatrixInversion(ESL_DMATRIX *A, ESL_DMATRIX **ret_Ai)
{
ESL_DMATRIX *Ai; /* RETURN: A^{-1} */
gsl_matrix_view Av; /* gsl view of matrix A */
gsl_matrix *LU; /* LU decomposition of A */
gsl_matrix *Aiv; /* gsl version of A^{-1} */
gsl_permutation *permute;
int signum;
int i,j;
Ai = esl_dmx_Alloc(A->n, A->m);
/* Invert U to get Ui, using LU decomposition.
*/
Av = gsl_matrix_view_array(A->mx[0], A->n, A->n);
LU = gsl_matrix_alloc(A->n, A->n);
Aiv = gsl_matrix_alloc(A->n, A->n); /* U^{-1}: inverse of U */
permute = gsl_permutation_alloc(A->n);
gsl_matrix_memcpy(LU, &Av.matrix);
if (gsl_linalg_LU_decomp(LU, permute, &signum) != 0) ESL_EXCEPTION(eslEUNKNOWN, "gsl failed");
if (gsl_linalg_LU_invert(LU, permute, Aiv) != 0) ESL_EXCEPTION(eslEUNKNOWN, "gsl failed");
gsl_matrix_free(LU);
gsl_permutation_free(permute);
/* recover the matrix from gsl.
*/
for (i = 0; i < A->n; i++)
for (j = 0; j < A->n; j++)
Ai->mx[i][j] = gsl_matrix_get(Aiv, i, j);
gsl_matrix_free(Aiv);
ret->Ai = Ai;
return eslOK;
}
#endif /*HAVE_LIBGSL*/
You can’t perform that action at this time.