Boost源码笔记：boost::multi_array
谢轩 /文
动机
C++是一门自由的语言，允许你自由的表达自己的意图，对不对? 所以我们既然可以new一个一维数组，也应该可以new出多维数组，对不对?先来看一个例子：
int* pOneDimArr = new int[10]; //新建一个10个元素的一维数组
pOneDimArr[0] = 0; //访问
int** pTwoDimArr = new int[10][20]; //错误！
pTwoDimArr[0][0] = 0; //访问
但是，很可惜，三四两行代码的行为并非如你所想象的那样——虽然从语法上它们看起来是那么“自然”。
这里的问题在于，new int[10][20]返回的并非int**类型的指针，而是int (*)[20]类型的指针（这种指针被称为行指针，对它“+1”相当于在数值上加上一行的大小（本例为20），也就是说，让它指向下一行），所以我们的代码应该像这样：
int (*pTwoDimArr)[20] = new int[i][20]; //正确
pTwoDimArr[1][2] = 0; //访问
注意pTwoDimArr的类型——int(*)[20]是个很特殊的类型，它不能转化为int**，虽然两者索引元素的语法形式一样，都是“p[i][j]”的形式，但是访问内存的次数却不一样，语义也不一样。
最关键的问题还是：以上面这种朴素的方式来创建多维数组，有一个最大的限制，就是：除了第一维，其它维的大小都必须是编译期确定的。例如：
int (*pNdimArr)[N2][N3][N4] = new int[n1][N2][N3][N4];
这里N2,N3,N4必须都是编译期常量，只有n1可以是变量，这个限制与多维数组的索引方式有关——无论多少维的数组都是线性存储在内存中的，所以：
pTwoDimArr[i][j] = 0;
被编译器生成的代码类似于：
*( (int*)pTwoDimArr+i*20+j ) = 0;
20就是二维数组的行宽，问题在于，如果允许二维数组的行宽也是动态的，这里编译器就无法生成代码（20所在的地方应该放什么呢？）。基于这个原因，C++只允许多维数组的第一维是动态的。
不幸的是，正由于这个限制，C++中的多维数组就在大多数情况下变成了有名无实的无用之物。我们经常可以在论坛上看到关于多维数组的问题，一般这类问题的核心都在于：如何模仿一个“完全动态的”多维数组。这里“完全动态”的意思是，所有维的大小都可以是动态的变量，而不仅是第一维。论坛上给出的答案不一而足，有的已经相当不错，但是要么缺乏可扩展性（即扩展到N维的情况），要么在访问元素的形式上远远脱离了内建的多维数组的访问形式，要么消耗了过多额外的空间。归根到底，我们需要的是一个类似这样的多维数组实现：
//创建一个int型的DIMS维数组，dim_sizes用于指定各维的大小，即n1*n2*n3
multi_array<int,DIMS> ma ( dim_sizes[n1][n2][n3] );
ma[i][j][k] = value; //为第i页j行k列的元素赋值
ma[i][j] = value; //编译错！
ma[i] = value; //编译错！
ma[i][j][k][l] = value;//编译错！
这样一个multi_array，能够自动管理内存，拥有和内建多维数组一致的界面，并且各维的大小都可以是变量——正符合我们的要求。看起来，实现这个multi_array并非难事，但事实总是出乎意料，下面就是对boost中已有的一个multi_array实现的剖析——你几乎肯定会发现一些出乎意料的（甚至是令人惊奇的）地方。
Boost中的多维数组实现——boost::multi_array
在Boost库中就有一个用于描述多维数组的功能强大的MultiArray库。它实现了一个通用、与标准库的容器一致的接口，并且具有与C++中内建的多维数组一样的界面和行为。正是基于这种通用性的设计，MultiArray库与标准库组件甚至用户自定义的泛型组件之间可以具有很好的兼容性，并能够很好的协同工作。
除此之外，MultiArray还提供了诸如改变大小、重塑（reshaping）以及对多维数组的视图访问等极为有用的特性，从而使MultiArray比其它描述多维数组的方式（譬如：std::vector< std::vector<...> > ）更为便捷、高效。
下面我们就逐层深入，去揭开boost::multi_array的神秘面纱——对示例程序进行调试、跟踪是分析库源代码最有效的手段之一。我们就从MultiArray文档中的示例程序入手：
// 略去头文件包含
int main () {
// 创建一个尺寸为3×4×2的三维数组
#define DIMS 3 //数组是几维的
typedef boost::multi_array<double,DIMS> array_type; // (1-1)
array_type A(boost::extents[3][4][2]); // (1-2)
// 为数组中元素赋值
A[1][2][0] = 120; // (1-3)
... ...
return 0;
}
在上述代码中，boost::multi_array的两个模板参数分别代表数组元素的类型和数组的维度。（1-2）处是三维数组对象的构造语句。boost::extents[3][4][2]则用于指明该数组各维的大小，这里的含义是：定义一个3*4*2的三维数组。你肯定对boost::extents心存疑惑——为什么对它可以连用“[]”？如果用多了一个或用少了一个“[]”又会如何？下面我就为你层层剥开boost::extents的所有奥秘——
extents——与内建数组一致的方式
boost::extents是一个全局对象，在base.hpp中：
typedef detail::multi_array::extent_gen<0> extent_gen;
... ...
multi_array_types::extent_gen extents; //注意它的类型！
可见extents的类型为extent_gen，后者位于extent_gen.hpp中：
// extent_gen.hpp
template <std::size_t NumRanges>
class extent_gen {
range_list ranges_; // （2-1）
... ...
extent_gen(const extent_gen<NumRanges-1>& rhs, const range& a_range)
{ // （2-2）
... //
}
extent_gen<NumRanges+1> operator[](index idx) //(2-3)
{ return extent_gen<NumRanges+1>(*this,range(0,idx)); }
//返回一个extent_gen，不过第二个模板参数增加了1
};
boost::extent_gen重载了operator[]操作符，但是，既然这里只有一个“[]”，为什么我们可以写“extents[n1][n2][n3][...]呢？继续看——
如果把boost::extents[n1][n2][n3]展开为操作符调用的方式就相当于：
boost::extents.operator[](n1).operator[](n2).operator[](n3);
boost::extents对象的类型是extent_gen<0> ——其NumRanges模板参数为0。extents.operator[](n1)调用（2-3），返回一个extent_gen<NumRanges+1>，也就是extent_gen<1>，由于这个类型也是extent_gen，所以仍然可以对它调用operator[](n2)，返回extent_gen<2>，然后再对该返回值调用operator[](n3)，最终返回一个extent_gen<3>类型的对象。
这里，每一重operator[]调用都转到（2-3），在那里将参数idx以range包装一下再转发给构造函数（2-2），注意此时调用的是extent_gen<NumRange+1>类型的构造函数。至于range(0,idx)则表示一个[0,idx)的下标区间。
总结一下extents的基本工作方式——每对它调用一次operator[]，都会返回一个extent_gen<NumRange+1>类型的对象，所以，对于boost::extents[n1][n2][n3]，依次返回的类型是：
extent_gen<1> => extent_gen<2> => extent_gen<3>
最后一个也是最终的返回类型——extent_gen<3>。其模板参数3表示这是用于三维数组的下标指定。其成员ranges_中，共有[0，n1）、[0，n2）、[0，n3）三组区间。这三组区间指定了我们定义的multi_array对象的三个维度的下标区间，值得注意的是这些区间都是前闭后开的（这和STL的迭代器表示的区间一样）。当boost::extents准备完毕后，就被传入multi_array的构造函数，用于指定各维的下标区间：
// multi_array.hpp
explicit multi_array(const extent_gen<NumDims>& ranges):
super_type((T*)initial_base_,ranges){
allocate_space(); // （2-5）
}
这里，multi_array接受了ranges参数中的信息，取出其中各维的下标区间，然后保存，最后调用allocate_space()来分配底层内存。
使用extent_gen的好处——do things right!
使用boost::extents作参数的构造过程和内建多维数组的方式一致，简练直观，语义清晰。首先，boost::extents使用“[]”，能让人很容易想到内建多维数组的声明，也很清晰地表达了每个方括号中数值的含义——表明各维度的下标区间；最关键的还是，使用boost::extents，可以防止用户一不小心写出错误的代码，例如：
multi_array<int,3> A(boost::extents[3][4][2][5]);
//错！不能用四维的下标指定来创建一个三维数组!
上面的语句是完全错误的，因为mult_array是个三维数组，而boost::extents后面却跟了四个“[]”，这显然是个错误。这个错误被禁止在编译期，这是由于在语法层面，multi_array<int,3>的构造函数只能接受extent_gen<3>类型的参数，而根据我们前面对extents的分析，boost::extents[3][4][2][5]返回的却是extent_gen<4>类型的对象，于是就会产生编译错误。这种编译期的强制措施阻止了用户一不小心犯下的错误（如果你正在打瞌睡呢？），也很清晰明了地表达（强制）了语义的需求。
另一种替代方案及其缺点
另外，还有一种声明各维大小的替代方式，就是使用所谓的Collection Concept，例如：
// 声明一个shape（“形状”），即各个维度的size
boost::array<int,3> shape = {{ 3, 4, 2 }};
array_type B(shape); //3*4*2的三维数组
这种构造方式将调用multi_array的第二个构造函数：
// multi_array.hpp
template <class ExtentList>
explicit multi_array( ExtentList const& extents ) :
super_type((T*)initial_base_,extents) {
boost::function_requires< // （2-4）
detail::multi_array::CollectionConcept<ExtentList> >();
allocate_space(); // （2-6）
}
这个构造函数的形参extents只要是符合collection concept就可以了——shape的类型为boost::array，当然符合这个concept。这个构造函数的行为与接受extents_gen的构造函数是一样的——仍然是先取出各维的range保存下来，然后分配底层内存。至于（2-4）处的代码，则是为了在编译期静态检查模板参数ExtentList是否符合Collection concept，实现细节在此不再赘述。
把这种方式与使用extent_gen的方式作一个简单的比较，很容易就看出优劣：采用这种方式，就不能保证编译期能够进行正确性的检查了，例如：
boost::array<int,4> shape = {{3,4,2,5}}; //一个四维数组的shape
multi_array<int,3> A(shape); // 竟然可以通过编译！！
这里，用一个四维的shape来指定一个三维multi_array是不恰当的，但是却通过了编译，这是由于该构造函数对它的参数extents没什么特殊要求，只是把它作为一个普通的collection来对待，构造函数会根据自己的需要从extents中取出它所需要的各维下标区间——A是三维数组，于是构造函数从shape中取出前三个数值作为A三个维度的下标区间，而不管shape究竟包含了几个数值。这样的语句在语义上是不清晰甚至错误的。但是既然这样的构造函数存在，设计者自然有他的道理，文档中就明确的表明，这个构造函数最大的用处就是编写维度无关（dimension-independent）的代码，除此之外multi_array库默认为前一种构造函数。
multi_array的架构
无论采用哪一种构造函数，构造的流程都是相似的——将一系列下标区间传入基类的构造函数中去，基类构造完成之后就调用相同的allocate_space()函数（见（2-5）和（2-6）处），allocate_space，顾名思义，用于为多维数组的元素分配空间。
然而，在这个看似简单的构造过程背后，却隐藏了两三层层次结构，每层的结构和功能都有所不同，担任了不同的角色，有些层次则能够单独拉出来复用于其他地方。
下面我们就来看看multi_array的架构到底长什么样J首先，multi_array继承自multi_array_ref：
// multi_array_ref.hpp
template <typename T, std::size_t NumDims>
class multi_array_ref : //multi_array的基类！！
public const_multi_array_ref<T,NumDims,T*>
{
typedef const_multi_array_ref<T,NumDims,T*> super_type;
... ...
explicit multi_array_ref(T* base, //指向数组存储空间的指针
const extent_gen<NumDims>& ranges): //下标区间
super_type(base,ranges) //把初始化的任务转发给基类（3-1）
{ }
... ...
};
而multi_array_ref又以const_multi_array_ref为基类：
// multi_array_ref.hpp
class const_multi_array_ref : //multi_array_ref的基类!管理底层存储！
public multi_array_impl_base<T,NumDims>
{
... ...
explicit const_multi_array_ref(TPtr base,
const extent_gen<NumDims>& ranges) :
base_(base), storage_(c_storage_order()) // （3-2）
{ init_from_extent_gen(ranges); }
... ...
storage_order_type storage_;//支持多种存储策略！（3-3）
};
multi_array的构造之路途经（3-1）处（multi_array_ref的构造函数），延伸至（3-2）处（const_multi_array_ref的构造函数）——这里看似一个终结，因为再没有参数传递给 const_multi_array_ref的基类multi_array_impl_base了。但是心中还是疑惑：为什么会有如此多层的继承结构？这样的类层次结构设计究竟有什么玄机呢？
多层继承的奥秘——访问界面&&复用性
转到基类const_multi_array_ref的声明，似乎可以看出一些端倪：
template< ... >
class const_multi_array_ref {
... ...
//和所有的STL容器一致的迭代器界面！！
const_iterator begin() const;
const_iterator end() const;
... ...
//和std::vector一致的元素访问界面！！
const_reference operator[](index i) const;
... ...
};
看到上面这些声明，是不是有些面熟？STL！对，这些成员函数的声明是与STL中container concept完全一致的。也就是说，这里提供了与STL容器一致的访问界面，所谓与STL的兼容性正是在这里体现出来了。而const_multi_array_ref更是“类如其名”，const_multi_array_ref中所有访问元素、查询数组信息等成员函数都返回const的reference或iterator。而反观multi_array_ref的声明，其中只比const_multi_array_ref多了访问元素、查询数组信息的对应的non-const版本成员函数。
那么const_multi_array_ref的基类multi_array_impl_base的职责又是什么呢?multi_array_impl_base是属于实现细节的，它的作用只是根据数组信息（const_multi_array_ref中的成员变量）计算偏移量、步长等，也就是把多维的下标最终转化为一维偏移量。而multi_array_impl_base的基类——value_accessor_n或者value_accessor_one——的功能则是提供一个对原始数据的访问。这种访问方式是大家所熟悉的把多维索引转化为一维偏移量的方式——凡是写过用一维数组模拟多维数组的人都应该清楚。至于为什么要有value_accessor_n和value_accessor_one两个版本，后文会详细阐释。
至此，multi_array的大致架构就已经浮现出来了：
multi_array -> multi_array_ref -> const_multi_array_ref -> multi_array_impl_base -> value_accessor_n/value_accessor_one
其中每一层都担任各自的角色：
¨ multi_array : 为数组元素分配空间，将各种操作转发至基类。
¨ multi_array_ref : 提供与STL容器一致的数据访问界面。也可以独立出来作为一个adapter使用。
¨ const_multi_array_ref : 提供const的STL数据访问界面。也可以作为一个const adapter使用。
¨ multi_array_impl_base及其基类 :最底层实现，提供一组对原始数据的基本操作。
这种架构看似复杂，却提供了很高的复用性，其中的(const_)multi_array_ref类都可以独立出来作为一个adapter使用——例如：
int a[24]; //一维的10个元素数组，位于栈上！
//把一维数组a“看成”一个3*4*2的三维数组：
multi_array_ref<int,3> arr_ref(a,boost::extents[3][4][2]);
arr_ref[i][j][k] = value; //和multi_array一样的使用界面
很简单吧！从此你就不用辛辛苦苦的去模拟多维数组了。即使是位于栈上的一维数组，也可以把它“看成”是一个多维数组。倘若你不想让multi_array来自动分配内存的话，你可以自行分配数组（可以位于栈上或堆上）然后用multi_array_ref把它包装成一个多维的数组。
multi_array的存储策略
接下来，就来看看multi_array的存储策略，例如：C风格的多维数组存储方式是按行存储，而fortran恰恰相反，是按列存储，甚至，用户可能有自己的存储策略要求。那么，如何支持多种风格的存储策略呢？秘密就在于代码（3-3）处，const_multi_array_ref的成员storage_——其类型为storage_order_type，下面的声明指出了storage_order_type的“本来面目”——general_storage_order<NumDims>：
// multi_array_ref.hpp
... ...
typedef general_storage_order<NumDims> storage_order_type;
... ...
// storage_order.hpp
template <std::size_t NumDims>
class general_storage_order {
general_storage_order(const c_storage_order&){ //（4-1）
for (size_type i=0; i != NumDims; ++i)
{ ordering_[i] = NumDims - 1 - i; }
ascending_.assign(true); //缺省为升序排列
}
... ...
boost::array<size_type,NumDims> ordering_; //各维的优先顺序
boost::array<bool,NumDims> ascending_;//升序排列还是降序排列
};
在（4-1）处的构造函数中，ordering_和ascending_是两个数组，缺省情况下，当函数（4-1）执行完毕后，ordering_中的元素应当是{NumDims-1, NumDims-2,...1,0}，也就是说，NumDims-1维在先，然后是NumDims-2维，...，如果将这些元素作为各维度存储顺序的标识——具有较小ordering_值的维度优先——那么就这和C语言中的存储方式就完全一致了，而ascending_勿庸置疑就是用来表明各维度是否升序存储。其实general_storage_order还有一个模板构造函数，它是为了支持更为一般化的存储策略（例如fortran的按列存储或用户自定义的存储策略）。这里不作详述。
除了存储策略，const_multi_array_ref的构造还通过调用init_from_extent_gen函数，将extents中的内容取出来进行处理，并以此设定其它若干表述多维数组的变量（（3-3）处其它一些变量），具体细节不再赘述。
现在关于一个多维数组的所有信息都已经准备齐备，可谓“万事具备，只欠‘空间’”。multi_array下面要做的就是调用前面提到的allocate_space来为数组中的元素分配空间了。
// multi_array.hpp
void allocate_space() {
... ...
base_ = allocator_.allocate(this->num_elements(),no_hint);
... ... std::uninitialized_fill_n(base_,allocated_elements_,T());
}
原来，在底层，存储仍然是退化为一维数组的存储： allocate_space使用allocator_分配一块连续的空间用以存储元素，其中num_elements()返回的就是数组各维度的大小的乘积，即数组的总元素个数。分配完之后，就将首指针赋给表述数组基地址的成员base_，然后std::uninitialized_fill_n负责把该数组进行缺省初始化。至此multi_array的构造工作终于大功告成了。
一致性界面——GP的灵魂
multi_array的另一重要特性就是支持与内建多维数组相同的访问方式，也就是说，multi_array支持以连续的“[]”来访问数组元素。就以文章开头给出的示例代码（1-3）处的赋值语句为例，让我们看看multi_array是如何支持这种与内建数组兼容的访问方式的。
// multi_array_ref.hpp
// 使用operator[]来访问元素,返回类型reference是什么呢？不是T&!
reference operator[](index idx) {
return super_type::access(boost::type<reference>(),
idx,origin(),this->shape(),this->strides(),
this->index_bases());
}
这个调用转入了value_accessor_n::access(...)之中：
// base.hpp
// in class value_accessor_n
template <typename Reference, typename TPtr>
Reference access(boost::type<Reference>,
index idx,TPtr base,const size_type* extents,
const index* strides,const index* index_base)
{
TPtr newbase = base + idx * strides[0];
return Reference(newbase,extents+1,
strides+1,index_base+1);
}
这个连续调用operator[]的过程和extend_gen是很类似的——每调用一层就返回一个“proxy(替身)”，这个替身并非实际元素，因为这时候后面还跟有“[]”，所以对这个替身也应该能调用operator[]，从而使这个过程能够连续下去直到后面没有“[]”为止”。
举个例子，如果以A[x1][x2][x3]方式访问A中的元素，就相当于
A.operator[x1].operator[x2].operator[x3] //连续调用“[]”
这三次operator[]调用返回的类型依次为：
sub_array<T,2> -> sub_array<T,1> -> T&
注意sub_array的第二个模板参数（维度）递减的过程，当递减到0时，表明已经递归到了最后一个维度，真正的访问开始了，所以最后一次调用“[]”返回的恰好是对元素的引用（这就刚好印证了前面所说的——当且仅当“[]”的个数和数组的维数相同的时候，才能够取出元素，否则你得到的返回值要么还是sub_array<...>（而不是元素），要么会由于试图在T&上继续调用“[]”而编译失败）那么，这一切究竟是如何做到的呢？
从上面的递归过程，我们可以轻易看出：真正对元素进行访问并返回T&的任务交给了sub_array<T,1>。那么这个sub_array到底是个什么东东？
sub_array的秘密
sub_array的定义在sub_array.hpp中：
// sub_array.hpp
template <typename T, std::size_t NumDims>
class sub_array : public const_sub_array<T,NumDims,T*>;
template <typename T, std::size_t NumDims, typename TPtr>
class const_sub_array :
public multi_array_impl_base<T,NumDims>;
//base.hpp
template <typename T, std::size_t NumDims>
class multi_array_impl_base:public
value_accessor_generator<T,mpl::size_t<NumDims> >::type ;
唔，原来sub_array最终继承自multi_array_impl_base，后者的基类型由value_accessor_generator来生成——它会根据NumDims的不同而生成不同的基类型：
// base.hpp
template <typename T, typename NumDims>
struct value_accessor_generator {
... ...
typedef typename //如果NumDims为1，则类型为value_accessor_one
mpl::apply_if_c<(dimensionality == 1),
choose_value_accessor_one<T>,
choose_value_accessor_n<T,dimensionality>
>::type type; //把这个类型作为multi_array_impl_base的基类！
};
很显然，如果dimensionality == 1，那么“::type”就是value_accessor_one<T>，而只有对value_accessor_one使用“[]”才能返回T&，否则，“::type”被推导为value_accessor_n，这只是个“proxy”而已，对它运用“[]”只会返回sub_array<T,NumDims-1>，从而可以对它继续这个连续调用“[]”的过程，直到dimensionality 降为1，sub_array的基类才变成了value_accessor_one，这时候再使用“[]”就会返回T&！
取出元素
根据上面的分析，取元素的任务最终交给value_accessor_one，其成员函数access如下：
template <typename Reference, typename TPtr>
Reference access(boost::type<Reference>,index idx,TPtr base,
const size_type*,const index* strides,
const index*) const {
return *(base + idx * strides[0]); //终于取出了数据！
}
这里，access的返回类型Reference就是T&，即数组中元素类型的引用，从而可以将指定元素的引用返回，达到访问数组元素的目的。看到这里，multi_array以内建数组访问方式访问数组元素的过程基本已经弄清楚了，至于其中一些细节，尤其是计算地址的细节，譬如：偏移量的计算、步长的使用等，皆已略去了。
现在也许你会有这样的疑惑：以内建多维数组的“[][][]...”访问方式来访问multi_array的元素的能力真的如此重要吗？费这么大力气、写这么多代码还不如以多参数的方式重载operator[]呢（“[i,j,...]”形式）！这么大代价真的值得吗？值得！这是勿庸置疑的。以内建数组访问方式访问multi_array元素的能力最重要的表现就是：可以把multi_array作为一个内建多维数组来对待，与内建类型的一致性是你的类能够和泛型算法合作的关键。举个例子：用户编写了一个函数模板，
template <typename ReturnType, typename _3DArray>
ReturnType someAlgo(_3Array& mda){//可以作用于内建多维数组
... ...
mda[x][y][z] = mda[z][x][y]; //对内建多维数组的访问形式！
... ...
}
因为有了以内建数组访问方式访问multi_array元素的能力，这个someAlgo算法就可以同时应用在内建数组和multi_array上（否则用户就得为multi_array提供一个单独的重载版本），如此一来，代码的可重用性、可扩展性都大大提高了。
效率
效率是C++永恒的主题，MultiArray库也不例外。执行时间效率上，纵观MultiArray库对数组元素的访问代码，虽然函数调用嵌套层数甚多，但多数调用都是简单的转发或者编译期的递归，在一个成熟的现代C++编译器下，这些转发的函数调用代码应该可以很轻易地被优化掉，所以在效率上几乎没有什么损失。在空间方面，由于大量运用模板技术，基本能够在编译期决定的内容都已决定，没有为运行期带来不必要的空间上的负担。总的看来，Boost.MultiArray库的确是难得的高效又通用的多维数组的实现。
结语
本文只是将multi_array最基本的功能代码做了一个扼要的分析，正如文章开始所说，multi_array还有许多很有用的特性，如果读者想充分了解multi_array的运作机制与实现技巧，就请深入完整地分析multi_array的代码吧，相信一定会大有收获的！
