gsl lu разложение и инверсия для матрицы с плавающей точкой

Question

Из-за ограничения памяти, мне нужно использовать gsl_matrix_float вместо gsl_matrix который хранит данные типа double. Тем не менее, я хочу использовать gsl_linalg_LU_decomp а также gsl_linalg_LU_invert которые только поддерживают gsl_matrix, И я не нашел какой-то другой метод, который поддерживает декомпозицию и инверсию версии с плавающей точкой в ​​gsl.

Есть ли способ решить эту дилемму? Или я могу только переводить с поплавка на двойной, а потом обратно? Заранее спасибо!

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));
}