gsl数据类型之向量

最新推荐文章于 2024-09-13 08:23:41 发布

Augusdi

最新推荐文章于 2024-09-13 08:23:41 发布

阅读量1.6k

点赞数

分类专栏： GSL-GNU scientific Library

GSL-GNU scientific Library 专栏收录该内容

94 篇文章 4 订阅

订阅专栏

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}  