Вот набросал код, иллюстрирующий проблему. Перемножаю две матрицы 2x2, значения элементов матрицы произвольны.
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
gsl_matrix *MatrixMultiply(gsl_matrix *mat1, gsl_matrix *mat2)
{
if (!mat1 || !mat2)
{
printf("one empty matrix detected!\n");
return NULL;
}
gsl_matrix *r=gsl_matrix_calloc(mat1->size1, mat2->size2);
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, mat1, mat2, 0.0, r);
return r;
}
void MatrixPrint(gsl_matrix *mat)
{
if (!mat)
{
printf("empty matrix!\n");
return;
}
/* gsl_matrix_fprintf(stdout, mat, "%.10e"); */
size_t i=0;
do
{
size_t j=0;
do
printf("%17.10e ", gsl_matrix_get(mat, i, j));
while(++j<mat->size2);
putc('\n', stdout);
}
while(++i<mat->size1);
}
gsl_matrix *InvertMatrix(gsl_matrix *mat)
{
if(!mat || mat->size1!=mat->size2)
return NULL;
int s;
gsl_matrix *r=gsl_matrix_calloc(mat->size1, mat->size1);
gsl_permutation *perm=gsl_permutation_alloc(mat->size1);
gsl_linalg_LU_decomp(mat, perm, &s);
gsl_linalg_LU_invert(mat, perm, r);
return r;
}
int main(int argc, char **argv)
{
gsl_matrix *T1=gsl_matrix_calloc(2, 2);
gsl_matrix *T2=gsl_matrix_calloc(2, 2);
gsl_matrix_set(T1, 0, 0, 7.0);
gsl_matrix_set(T1, 0, 1, 3.0);
gsl_matrix_set(T1, 1, 0, -2.0);
gsl_matrix_set(T1, 1, 1, 11.0);
gsl_matrix_set(T2, 0, 0, -3.0);
gsl_matrix_set(T2, 0, 1, 4.0);
gsl_matrix_set(T2, 1, 0, 8.0);
gsl_matrix_set(T2, 1, 1, 2.0);
MatrixPrint(T1);
puts("");
MatrixPrint(InvertMatrix(T1));
puts("");
MatrixPrint(MatrixMultiply(T1, InvertMatrix(T1)));
return 0;
}
Вот что на выходе
$ ./test
7.0000000000e+00 3.0000000000e+00
-2.0000000000e+00 1.1000000000e+01
1.3253012048e-01 -3.6144578313e-02
2.4096385542e-02 8.4337349398e-02
1.0000000000e+00 0.0000000000e+00
3.5045023120e-02 1.0014602093e+00
Компилирую так: gcc -o2 test.c -lgsl -o test (и ещё возможно с -lgslcblas). Проверял на 64-х разрядной генте и в виртуалке с 32-х разрядным squeeze. gsl что там, что там 1.15.
Обратная матрица считается правильно. Не нравится абсолютная погрешность (0,035 вместо 0) и относительная (1,0015 вместо 1). Я работаю с double же, всё элементарно считается на калькуляторе. По идее, погрешности должны быть этак в миллион раз меньше. А сейчас у меня диф. уравнения не решаются, возможно из-за этого.
PS Программу писал по мотивам этого треда.