2008年5月23日星期五

GSL笔记20080523

最近的一个程序要用到一些矩阵和向量的操作,主程序是用c写的,自己实现一些矩阵的操作倒也不算困难,但是有些运算就不太好实现了,而且质量也难以保证。在google上搜索了一番,于是就发现了GSL(GNU Scientific Library),GNU的好东西还真不少。

很快的翻了一下gsl的手册,写的比较详细,再配合它的例子,一般的使用应该不成问题。
照葫芦画瓢,写了一个求逆矩阵的程序,当然都是数值运算的,好在这次的程序没有符号运算。

#include
#include

int
main (void)
{
double a_data[] = { 0.3, 0, 0, 0,
0, 0.3, 0, 0,
0, 0, 0.3, 0,
0, 0, 0, 0.3 };

double b_data[] = { 1.0, 2.0, 3.0, 4.0 };

gsl_matrix_view m
= gsl_matrix_view_array(a_data, 4, 4);

gsl_vector_view b
= gsl_vector_view_array(b_data, 4);

gsl_vector *x = gsl_vector_alloc (4);
gsl_matrix *mi= gsl_matrix_alloc (4,4);

int s;

gsl_permutation * p = gsl_permutation_alloc (4);

gsl_linalg_LU_decomp (&m.matrix, p, &s);

gsl_linalg_LU_invert (&m.matrix, p, mi);

printf ("inv=\n");
gsl_matrix_fprintf (stdout, mi, "%f");

gsl_permutation_free (p);
return 0;
}


permutation没有查到准确的意思,从手册和例子中的用法猜测,意思大概是行列式前面的符号。
gsl中的线性代数部分在操作矩阵的时候,都是把矩阵分解成某种乘积形式,然后再做运算,上面程序用的是LU分解把m分解成一个下三角阵和上三角阵的乘积,之后再求逆。

1 条评论:

丽佳 说...
此评论已被作者删除。