浮点矩阵的gsl lu分解和求逆

Question

由于内存限制，我需要使用gsl_matrix_float 而不是gsl_matrix 来存储double 类型的数据。但是，我想使用仅支持gsl_matrix 的gsl_linalg_LU_decomp 和gsl_linalg_LU_invert。而且我没有找到其他支持gsl中浮点版本分解和反转的方法。

有什么办法可以解决这个困境吗？或者我只能从 float 转移到 double 然后再返回？提前致谢！

Answer 1

按照您的建议，您可能做的最好的事情是从 float 转换为 double 并返回。这是执行反转的示例代码（仅给出了基本组件 - 您必须填写空白）：

include <gsl/gsl_blas.h>
include <gsl/gsl_linalg.h>

void matrix_invert(gsl_matrix_float *, gsl_matrix_float *, int);

int main()
{
 gsl_matrix_float *X = gsl_matrix_float_alloc(N, N);
 gsl_matrix_float *invX = gsl_matrix_float_alloc(N, N);

 matrix_invert(X, invX, N); //invM = inv(I)
 return 0;
}

void matrix_invert(gsl_matrix_float *matrix, gsl_matrix_float *inverse, int N)
{
 int i=0,j=0,signum=0;
 gsl_matrix *DM = gsl_matrix_alloc(N, N);
 gsl_matrix *DM_I = gsl_matrix_alloc(N, N);

 for (i=0;i<N;i++)
  for (j=0;j<N;j++)
   gsl_matrix_set(DM, i, j, gsl_matrix_float_get(matrix,i,j));

  gsl_permutation *p = gsl_permutation_alloc(N);
  gsl_linalg_LU_decomp(DM, p, &signum);
  gsl_linalg_LU_invert(DM, p, DM_I);
  gsl_permutation_free(p);
  gsl_matrix_free(DM);

  for (i=0;i<N;i++)
   for (j=0;j<N;j++)
    gsl_matrix_float_set(inverse, i, j, gsl_matrix_get(DM_I,i,j));
}