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分解成一个下三角阵和上三角阵的乘积，之后再求逆。