你好,游客 登录 注册 搜索
背景:
阅读新闻

gsl数据类型之向量 - Augusdi的专栏

[日期:2013-04-15] 来源:  作者: [字体: ]

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

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

gsl_blockdouble
gsl_block_floatfloat
gsl_block_long_doublelong double
gsl_block_intint
gsl_block_uintunsigned int
gsl_block_longlong
gsl_block_ulongunsigned long
gsl_block_shortshort
gsl_block_ushortunsigned short
gsl_block_charchar
gsl_block_ucharunsigned char
gsl_block_complexcomplex double
gsl_block_complex_floatcomplex float
gsl_block_complex_long_doublecomplex long double

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

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

  1. 【gsl_block_double.h】  
  2. ……  
  3. struct gsl_block_struct  
  4. {  
  5.   size_t size;  
  6.   double *data;  
  7. };  
  8. typedef struct gsl_block_struct gsl_block;  
  9. ……  

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

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

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

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

 

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

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

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

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

  1. 【gsl_vector_double.h】  
  2. ……  
  3. gsl_vector *gsl_vector_alloc (const size_t n);  
  4. gsl_vector *gsl_vector_calloc (const size_t n);  
  5. gsl_vector *gsl_vector_alloc_from_block (gsl_block * b,  
  6.                                                      const size_t offset,   
  7.                                                      const size_t n,   
  8.                                                      const size_t stride);  
  9. gsl_vector *gsl_vector_alloc_from_vector (gsl_vector * v,  
  10.                                                       const size_t offset,   
  11.                                                       const size_t n,   
  12.                                                       const size_t stride);  
  13. 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的源代码,便容易明白两者的区别:

 

  1. 【init_source.c】  
  2. ......  
  3. TYPE (gsl_vector) *  
  4. FUNCTION (gsl_vector, alloc) (const size_t n)  
  5. {  
  6.   TYPE (gsl_block) * block;  
  7.   TYPE (gsl_vector) * v;  
  8.   ......  
  9.   v = (TYPE (gsl_vector) *) malloc (sizeof (TYPE (gsl_vector)));  
  10.   ......  
  11.   block = FUNCTION (gsl_block,alloc) (n);//分配新的数据块。  
  12.   ......        
  13.   v->data = block->data ;  
  14.   v->size = n;  
  15.   v->stride = 1;  
  16.   v->block = block;  
  17.   v->owner = 1;//新创建的数据块的所有者即该向量本身。  
  18.   return v;  
  19. }  
  20. ......  
  21. TYPE (gsl_vector) *  
  22. FUNCTION (gsl_vector, alloc_from_block) (TYPE(gsl_block) * block,   
  23.                                          const size_t offset,   
  24.                                          const size_t n,   
  25.                                          const size_t stride)  
  26. {  
  27.   TYPE (gsl_vector) * v;  
  28.   ......  
  29.       
  30.   v = (TYPE (gsl_vector) *) malloc (sizeof (TYPE (gsl_vector)));  
  31.   ......  
  32.   v->data = block->data + MULTIPLICITY * offset ;//offset表示新向量的第一个元素的地址相对于传入的block->data的偏移量。  
  33.   v->size = n;//设置新向量的位数  
  34.   v->stride = stride;//设置新向量中每个元素的大小  
  35.   v->block = block;  
  36.   v->owner = 0;//block是由外界传入的,不为新向量所有  
  37.   return v;  
  38. }  
  39. ......  
 

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

[c-sharp] view plaincopy
  1. 【templates_on.h】  
  2. ......  
  3. #define MULTIPLICITY 1  
  4. ......  
  5. #define CONCAT2x(a,b) a ## _ ## b   
  6. #define CONCAT2(a,b) CONCAT2x(a,b)  
  7. ......  
  8. #define FUNCTION(dir,name) CONCAT2(dir,name)  
  9. #define TYPE(dir) dir  
  10. ......  
 

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

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

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

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

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

[c-sharp] view plaincopy
  1. 【gsl_vector_double.h】  
  2. ......  
  3. typedef struct  
  4. {  
  5.   gsl_vector vector;  
  6. } _gsl_vector_view;  
  7. typedef _gsl_vector_view gsl_vector_view;  
  8. ......  
 

 

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

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

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

  1. 【subvector_source.c】  
  2. ......  
  3. QUALIFIED_VIEW(_gsl_vector, view)  
  4. FUNCTION(gsl_vector, subvector) (QUALIFIED_TYPE(gsl_vector) * v, size_t offset, size_t n)  
  5. {  
  6.   QUALIFIED_VIEW(_gsl_vector,view) view = NULL_VECTOR_VIEW;  
  7.   ......  
  8.   {  
  9.     TYPE(gsl_vector) s = NULL_VECTOR;  
  10.     s.data = v->data +  MULTIPLICITY * v->stride * offset ;  
  11.     s.size = n;  
  12.     s.stride = v->stride;  
  13.     s.block = v->block;  
  14.     s.owner = 0;  
  15.     view.vector = s;  
  16.     return view;  
  17.   }  
  18. }  
  19. ......  

其中宏的定义如下:

  1. 【templates_on.h】  
  2. ......  
  3. #define TYPE(dir) dir  
  4. #define VIEW(dir,name) CONCAT2(dir,name)  
  5. #define QUALIFIED_TYPE(dir) TYPE(dir)  
  6. #define QUALIFIED_VIEW(dir,name) CONCAT2(dir,name)  
  7. ......  
  8. 【views.h】  
  9. #define NULL_VECTOR {0, 0, 0, 0, 0}  
 

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

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

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





收藏 推荐 打印 | 录入:admin | 阅读:
相关新闻