1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
|
/* 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*/
|
