Перемножение матриц, gsl (Неверный результат)

Модератор: Модераторы разделов

Аватара пользователя
ZyX
Сообщения: 355
ОС: Gentoo

Перемножение матриц, gsl

Сообщение ZyX »

Код:

/* test.c */ /*{{{ Директивы #include, внешние и глобальные переменные */ #include <stdio.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_linalg.h> #include <gsl/gsl_blas.h> /*{{{ _cplusplus */ #ifdef _cplusplus extern "C" { #endif /*}}}*/ gsl_matrix *MatrixMultiply(gsl_matrix *mat1, gsl_matrix *mat2) { gsl_matrix *r=gsl_matrix_calloc(mat1->size1, mat2->size2); gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., mat1, mat2, 0.0, r); /* * size_t i=r->size1; * while(i--) * { * size_t j=r->size2; * while(j--) * { * size_t k=mat1->size2; * double *cur=gsl_matrix_ptr(r, i, j); * while(k--) * [> gsl_matrix_set(r, i, j, <] * [> gsl_matrix_get(r, i, j)+ <] * *cur+=( * gsl_matrix_get(mat1, i, k)* * gsl_matrix_get(mat2, k, j)); * } * } */ return ®; } void MatrixPrint(gsl_matrix *mat) { /* 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->size1!=mat->size2) return (0); 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 ®; } /*MAIN: {{{*/ int main(int argc, char **argv) { gsl_matrix *T=gsl_matrix_calloc(3, 3); long double nu=0.27; long double beta=0.; gsl_matrix_set(T, 0, 0, 0.346-0.259*nu-0.551*beta); gsl_matrix_set(T, 0, 1, -0.344-0.268*nu-0.571*beta); gsl_matrix_set(T, 0, 2, 0.793*nu-0.373*beta); gsl_matrix_set(T, 1, 0, 0.021-0.624*nu+0.306*beta); gsl_matrix_set(T, 1, 1, -0.020-0.646*nu+0.316*beta); gsl_matrix_set(T, 1, 2, -0.440*nu-0.898*beta); gsl_matrix_set(T, 2, 0, 0.630+0.163*nu+0.292*beta); gsl_matrix_set(T, 2, 1, -0.609+0.169*nu+0.302*beta); gsl_matrix_set(T, 2, 2, -0.421*nu+0.235*beta); MatrixPrint(T); puts(""); MatrixPrint(InvertMatrix(T)); puts(""); MatrixPrint(MatrixMultiply(T, InvertMatrix(T))); return (0); } /*}}}*/ /*{{{ _cplusplus */ #ifdef _cplusplus } #endif /*}}}*/ /* vim: fdm=marker:fenc=utf-8:ts=4:tw=80:ft=c */

zyx@zyx-desktop

(zyx:*) % gcc -lgsl test.c (zyx:*) % ./a.out 2.7607000000e-01 -4.1636000000e-01 2.1411000000e-01 -1.4748000000e-01 -1.9442000000e-01 -1.1880000000e-01 6.7401000000e-01 -5.6337000000e-01 -1.1367000000e-01 -6.0751906991e-01 -2.2760744287e+00 1.2344658579e+00 -1.3123304835e+00 -2.3809992182e+00 1.6535825688e-02 2.9018535777e+00 -1.6953760543e+00 -1.5595447808e+00 1.0000000000e+00 0.0000000000e+00 0.0000000000e+00 -2.8605144439e-02 -2.8961857607e+00 -1.6449198815e-02 2.0708673491e-02 8.1675344277e+00 3.5365321655e+00

maxima@zyx-desktop

(%i31) T; [ 0.27607 - 0.41636 0.21411 ] [ ] (%o31) [ - 0.14748 - 0.19442 - 0.1188 ] [ ] [ 0.67401 - 0.56337 - 0.11367 ] (%i32) iT; /* Эта матрица просто тупо скопирована из вывода ./a.out. */ [ - 0.60751906991 - 2.2760744287 1.2344658579 ] [ ] (%o32) [ - 1.3123304835 - 2.3809992182 0.016535825688 ] [ ] [ 2.9018535777 - 1.6953760543 - 1.5595447808 ] (%i33) T.iT; (%o33) [ 1.000000000001353 - 2.763000939154381E-11 - 1.009070604851559E-11 ] [ ] [ 1.6368573163561E-12 0.99999999999796 5.687006421339902E-12 ] [ ] [ 2.19690932112826E-12 - 3.847200336082324E-11 .9999999999988665 ]

Суть претензий следующая: как видно, gsl правильно считает обратные матрицы, но
почему-то проваливается на умножении. Никто не может сказать, в чём дело?

PS: Про погрешность округления я зная, но не настолько же она огромна.
Спасибо сказали:
Аватара пользователя
Фантом
Сообщения: 453
ОС: openSUSE

Re: Перемножение матриц, gsl

Сообщение Фантом »

ZyX писал(а):
23.04.2009 18:41
Суть претензий следующая: как видно, gsl правильно считает обратные матрицы, но
почему-то проваливается на умножении. Никто не может сказать, в чём дело?

PS: Про погрешность округления я зная, но не настолько же она огромна.


Во-первых, при перемножении прямой и обратной матриц неизбежно вылезают классические вычислительные грабли: сильнейший рост относительной погрешности при вычитании двух близких положительных чисел. Т.е. проблема не столько в точности умножения матриц вообще, сколько в точности вычислений в данной конкретной задаче. Попробуйте умножить две произвольных матрицы друг на друга - "провал" будет не таким заметным.

Во-вторых, погрешность округления не так уж и мала. Насколько я помню, там используется double, т.е. абсолютная погрешность результата не может быть лучше 10^{-15}, а упомянутое первым обстоятельство вполне способно поднять ее еще на два-три порядка.
Спасибо сказали:
Аватара пользователя
Vital86
Сообщения: 79
ОС: Debian SID, Mandriva 2010

Re: Перемножение матриц, gsl

Сообщение Vital86 »

С точностью все норм, я сам часто программировал инженерные задачи, так такая точность зачастую и не требуется, а чтоб получались 1 и 0 - просто округли до 10 разряда при выводе.
Спасибо сказали: