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的庐山真面目：[cpp] view plaincopy【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提供了三个函数：[c-sharp] view plaincopy【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为例进行分析，它的定义如下：[cpp] view plaincopy【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个关于向量内存管理的函数：[cpp] view plaincopy【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的源代码，便容易明白两者的区别： [cpp] view plaincopy【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，这些宏的代码如下：[c-sharp] view plaincopy【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中包含一些操作向量元素的函数，经常用到的有以下几个：[cpp] view plaincopy【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的一个简单封装：[c-sharp] view plaincopy【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个元素构建一个新的子向量。这一功能可以从它的源码中看出来：[cpp] view plaincopy【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;
}
}
......

其中宏的定义如下：[cpp] view plaincopy【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