Dec 23, 2008

GSL 学习笔记

1.    In linux, the sizes of int, float, double are 4 4 8
2.    Macro 全大写,如GSL_SIGN, 其他小写
3.    函数返回是否成功 int, 0代表成功
GSL_EFAILED 失败  (=-1)
4.    gsl 常数:和exp 相关的,以及 nan, inf
5.    gsl 基本函数,有些是为了避免 overflow
gsl_hypot:
hypotenuse
(直角三角形的)斜边, 弦
gsl_log1p 计算 log(1+x) for small x

6.    Small integer pow,
— Function: double gsl_pow_2 (const double x)
— Function: double gsl_pow_3 (const double x)
— Function: double gsl_pow_4 (const double x)
— Function: double gsl_pow_5 (const double x)
— Function: double gsl_pow_6 (const double x)
— Function: double gsl_pow_7 (const double x)
— Function: double gsl_pow_8 (const double x)
— Function: double gsl_pow_9 (const double x)

7.    Special function: gsl_sf.h
Mode 调整计算精度和速度,参见 section 7.3
typedef struct
{
double val;
double err;
} gsl_sf_result;

typedef struct
{
double val;
double err;
int
e10;
} gsl_sf_result_e10;

函数名 _e 表示有error 输出,需采用gsl_sf_result
同理 _e10,采用gsl_sf_result_e10

7.16 exponential function,不错,提防 overflow

gsl_sf_gamma
gsl_sf_lngamma
gsl_sf_beta
gsl_sf_psi
gsl_sf_psi_1   (trigamma)
gsl_sf_psi_n (int n, double x)   即 psi n 阶导数


8.    Blocks
我们通常不用 block, 而用 vector 和 matrix
Bolck 的作用是和文件交互,省去在程序中初始化的麻烦
typedef struct
{
size_t size;
double * data;
} gsl_block;

gsl_block * gsl_block_alloc (size_t n)   不初始化
gsl_block * gsl_block_calloc (size_t n)  初始化为 0
void gsl_block_free (gsl_block * b)
基本上 gsl 的数据结构都手动需要 allocate 和 free

相应于 block 有
gsl_block_float 等等

gsl_vector * gsl_vector_alloc (size_t n)
gsl_vector * gsl_vector_calloc (size_t n)   initializes all the elements of the vector to zero.
void gsl_vector_free (gsl_vector * v)

gsl_vector_get(v,i)     v->data[i*v->stride]           v->stride 代表 element 之间在 memory block 中的距离
gsl_vector_set(v,i,x)   v->data[i*v->stride]=x.


void gsl_vector_set_all (gsl_vector * v, double x)
void gsl_vector_set_zero (gsl_vector * v)


int gsl_vector_fwrite (FILE * stream, const gsl_vector * v)
int gsl_vector_fread (FILE * stream, gsl_vector * v)
int gsl_vector_fprintf (FILE * stream, const gsl_vector * v, const char * format)
int gsl_vector_fscanf (FILE * stream, gsl_vector * v)

typedef struct
{
gsl_vector vector;
} _gsl_vector_const_view;

typedef const _gsl_vector_const_view gsl_vector_const_view;

从 double * 获得view
gsl_vector_view gsl_vector_view_array (double * base, size_t n)
gsl_vector_view gsl_vector_view_array_with_stride (double * base, size_t stride, size_t n)

View 可以处理 sub vector,以及带间隔的 vector,但是没办法像 matlab 一样,index 可以是一向量

9.    Matrix

原则上,Vector 和 matrix 在内存中是连续的
但是 vector 有 stride,可以间隔固定
Matrix 有 tda,但是row 是连续的,即row 没有stride

同样,对于vector 和 matrix 有 float, int 版本
一般就加后缀
Gsl_vector_float

Index 顺序和 matlab 一样,(ith row, jth col)

typedef struct
{
size_t size1;
size_t size2;
size_t tda; //physical row-length of the matrix.
double * data;
gsl_block * block;
int owner;
} gsl_matrix;
Eg. (3,4) however, tda = 8
00 01 02 03 XX XX XX XX
10 11 12 13 XX XX XX XX
20 21 22 23 XX XX XX XX
Matrices are stored in row-major order, meaning that each row of elements forms a contiguous block in memory.
m->data[i * m->tda + j]

gsl_matrix_view gsl_matrix_submatrix (gsl_matrix * m, size_t k1, size_t k2, size_t n1, size_t n2)
gsl_matrix_const_view gsl_matrix_const_submatrix (const gsl_matrix * m, size_t k1, size_t k2, size_t n1, size_t n2)
-->
m'(i,j) = m->data[(k1*m->tda + k2) + i*m->tda + j]
同理,可以由 matrix 中获得 view
使用
&View.matrix 用于 gsl_matrix *

从 matrix 中获得
Col 和 row 的view
gsl_vector_const_view gsl_matrix_row (gsl_matrix * m, size_t i)
gsl_vector_const_view gsl_matrix_const_row (const gsl_matrix * m, size_t i)
gsl_vector_view gsl_matrix_column (gsl_matrix * m, size_t j)
gsl_vector_const_view gsl_matrix_const_column (const gsl_matrix * m, size_t j)

gsl_vector_view gsl_matrix_diagonal (gsl_matrix * m)

10.    BLAS Support
Level 1
Vector operations, e.g. y = \alpha x + y
Level 2
Matrix-vector operations, e.g. y = \alpha A x + \beta y
Level 3
Matrix-matrix operations, e.g. C = \alpha A B + C

11.    Random Number Generation  & Random Number Distributions

Random number generation 指 [0 1) 之间的随机数

而distribution 计算各种 distribution 的pdf 和随机数产生
Distribution 需要
const gsl_rng * r

12.    再接再厉

暂时学习到这里

找到一个 Gaussian multivariate
gsl_mvg.rar


0 comments: