LINUX.ORG.RU

gsl - матрицы перемножаются с большой погрешностью


0

2

Вот набросал код, иллюстрирующий проблему. Перемножаю две матрицы 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 Программу писал по мотивам этого треда.


Разобрался. Быдлоалгоритм в

InvertMatrix

по ссылке переписывал исходную матрицу.

gogi ()

А вот в фортране перемножение матриц искаропки.

aedeph ()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.