gsl数据类型之向量

    gsl中的向量与矩阵是基于一种底层的数据类型，即数据块（block）实现的。因此我们的分析，就从数据块开始。

    与复数一样，gsl中也有各种不同数据类型的数据块，详细列表如下：

gsl_block	double
gsl_block_float	float
gsl_block_long_double	long double
gsl_block_int	int
gsl_block_uint	unsigned int
gsl_block_long	long
gsl_block_ulong	unsigned long
gsl_block_short	short
gsl_block_ushort	unsigned short
gsl_block_char	char
gsl_block_uchar	unsigned char
gsl_block_complex	complex double
gsl_block_complex_float	complex float
gsl_block_complex_long_double	complex long double

    与上一节的分析一样，我们将着眼于double类型对应的gsl_block类型展开分析。

    让我们来看看gsl_block的庐山真面目：

【gsl_block_double.h】 …… struct gsl_block_struct { size_t size; double *data; }; typedef struct gsl_block_struct gsl_block; ……
      可见数据块本质上仍然是一个结构体。其中size_t在32位机上实际上就是unsigned int，在这里gsl_block->size表示该数据块中存储的双精度数的个数。而gsl_block->data则指向一块size*sizeof(double)大小的连续内存，这块内存就是数据块的主体。

      类似于c语言中动态内存管理，为了创建与销毁这样的数据块，gsl提供了三个函数：

【gsl_block_double.h】 …… gsl_block *gsl_block_alloc (const size_t n);//创建包含n个double元素的一块连续空间，并返回指向该内存空间头部的指针。 gsl_block *gsl_block_calloc (const size_t n);//功能与gsl_block_alloc类似，但申请到得空间中所以double元素都将初始化为零。 void gsl_block_free (gsl_block * b);//释放b所指向的数据块 ……

      正如你所预料的，这里的内存分配函数只是对c语言中malloc函数的简单封装，只不过其中包含了一些异常处理的判断语句。

 

      而gsl中向量与矩阵构建在数据块的基础上。因此与数据块一样，也有与前文列表中各种数据类型相对应的各种向量与矩阵。

      向量相关的结构体与函数的标识符以gsl_vector打头，后缀规则与数据块完全一样。这里仍以double型的gsl_vector为例进行分析，它的定义如下：

【gsl_vector_double.h】 …… typedef struct { size_t size; //表示向量的维数，即向量中元素的个数。 size_t stride; //向量中每个元素所占的内存大小。 double *data; //为指向向量中第一个数据的数据指针，实际上就是下面的block中的data指针。 gsl_block *block; //指向向量引用的数据块。 int owner; //所有者标识符。 } gsl_vector; ……

      值得注意的是，不同的向量可以指向相同的一块数据块。因此结构体里需要一个所有者标识符。如果owner=1，则表示block归该向量所有，当向量内存被释放时，该block也将被释放；而owner=0时，则表示所有者为其他向量block，因此在本向量被释放时，该block不会被释放。
      由于考虑到block复用的问题，gsl一共提供了5个关于向量内存管理的函数：

【gsl_vector_double.h】 …… gsl_vector *gsl_vector_alloc (const size_t n); gsl_vector *gsl_vector_calloc (const size_t n); gsl_vector *gsl_vector_alloc_from_block (gsl_block * b, const size_t offset, const size_t n, const size_t stride); gsl_vector *gsl_vector_alloc_from_vector (gsl_vector * v, const size_t offset, const size_t n, const size_t stride); void gsl_vector_free (gsl_vector * v);

      其中gsl_vector_alloc、gsl_vector_calloc的功能与gsl_block_alloc、gsl_block_calloc的功能完全类似，调用这两个函数在创建新的向量同时也将创建新的block。而gsl_vector_alloc_from_block与gsl_vector_alloc_from_vector则将从已经存在的block创建新的向量。

      比较一下alloc与alloc_from_block的源代码，便容易明白两者的区别：

 

【init_source.c】 ...... TYPE (gsl_vector) * FUNCTION (gsl_vector, alloc) (const size_t n) { TYPE (gsl_block) * block; TYPE (gsl_vector) * v; ...... v = (TYPE (gsl_vector) *) malloc (sizeof (TYPE (gsl_vector))); ...... block = FUNCTION (gsl_block,alloc) (n);//分配新的数据块。 ...... v->data = block->data ; v->size = n; v->stride = 1; v->block = block; v->owner = 1;//新创建的数据块的所有者即该向量本身。 return v; } ...... TYPE (gsl_vector) * FUNCTION (gsl_vector, alloc_from_block) (TYPE(gsl_block) * block, const size_t offset, const size_t n, const size_t stride) { TYPE (gsl_vector) * v; ...... v = (TYPE (gsl_vector) *) malloc (sizeof (TYPE (gsl_vector))); ...... v->data = block->data + MULTIPLICITY * offset ;//offset表示新向量的第一个元素的地址相对于传入的block->data的偏移量。 v->size = n;//设置新向量的位数 v->stride = stride;//设置新向量中每个元素的大小 v->block = block; v->owner = 0;//block是由外界传入的，不为新向量所有 return v; } ...... 

代码中涉及到了一些宏，对于gsl_vector，这些宏的代码如下：

【templates_on.h】 ...... #define MULTIPLICITY 1 ...... #define CONCAT2x(a,b) a ## _ ## b #define CONCAT2(a,b) CONCAT2x(a,b) ...... #define FUNCTION(dir,name) CONCAT2(dir,name) #define TYPE(dir) dir ...... 

注意，其中FUNCTION(dir,name)宏将产生字符串dir_name。

      从源码中我们容易看出调用gsl_vector_alloc_from_block与gsl_vector_alloc_from_vector并重用block中的数据时，允许设置新向量第一个元素的起始地址、元素大小与元素个数。使用非常灵活。

      gsl中包含一些操作向量元素的函数，经常用到的有以下几个：

【gsl_vector_double.h】 ...... double gsl_vector_get (const gsl_vector * v, const size_t i); //返回向量v中第i+1个元素的值 void gsl_vector_set (gsl_vector * v, const size_t i, double x); //将向量v中第i+1个元素置为x ...... void gsl_vector_set_zero (gsl_vector * v);//将v的所有元素均置为0 void gsl_vector_set_all (gsl_vector * v, double x);//将v的所有元素均设置为x int gsl_vector_set_basis (gsl_vector * v, size_t i);//将v的第i+1个元素设置为1，其余元素为0 ...... 

      gsl中还提供了一种vector views的数据结构以便更灵活的描述向量。实际上这种神秘的view就是gsl_vector的一个简单封装：

【gsl_vector_double.h】 ...... typedef struct { gsl_vector vector; } _gsl_vector_view; typedef _gsl_vector_view gsl_vector_view; ...... 

 

      基于vector views，gsl实现了从向量中提取子向量的一系列函数，其中最简单的版本是 

gsl_vector_view gsl_vector_subvector (gsl_vector * v, size_t offset, size_t n)。

我们可以利用该函数从v向量中的第offset+1个元素开始提取n个元素构建一个新的子向量。这一功能可以从它的源码中看出来：

【subvector_source.c】 ...... QUALIFIED_VIEW(_gsl_vector, view) FUNCTION(gsl_vector, subvector) (QUALIFIED_TYPE(gsl_vector) * v, size_t offset, size_t n) { QUALIFIED_VIEW(_gsl_vector,view) view = NULL_VECTOR_VIEW; ...... { TYPE(gsl_vector) s = NULL_VECTOR; s.data = v->data + MULTIPLICITY * v->stride * offset ; s.size = n; s.stride = v->stride; s.block = v->block; s.owner = 0; view.vector = s; return view; } } ......

其中宏的定义如下：

【templates_on.h】 ...... #define TYPE(dir) dir #define VIEW(dir,name) CONCAT2(dir,name) #define QUALIFIED_TYPE(dir) TYPE(dir) #define QUALIFIED_VIEW(dir,name) CONCAT2(dir,name) ...... 【views.h】 #define NULL_VECTOR {0, 0, 0, 0, 0} 

可见由该函数创建的新向量只是对原始向量的一种引用，两者共享一块block。对新向量中数据的改变会同步的体现在原始向量中——这正是gsl中vector views数据结构的含义：它只是对原有数据的一种动态的重新描述。

      除此之外，gsl还用一系列丰富的操作向量的函数，更详细的信息可以参见

http://www.gnu.org/software/gsl/manual/html_node/Vectors.html