第29章 一个矩阵的设计(A Matrix Design)
目录
29.1 引言
29.1.1 基本的 Matrix 用法
29.1.2 Matrix 的要求
29.2 一个 Matrix 模板
29.2.1 构造和赋值(Construction and Assignment)
29.2.2 下标和分片(Subscripting and Slicing)
29.3 Matrix算术运算(Matrix Arithmetic Operations)
29.3.1 标量运算(Scalar Operations)
29.3.2 加法运算(Addition)
29.3.3 乘法(Multiplication)
29.4 Matrix实现(Matrix Implementation)
29.4.1 slice()
29.4.2 Matrix分片(Matrix Slices)
29.4.3 Matrix_ref
29.4.4 Matrix 列表初始化
29.4.5 Matrix 访问
29.4.6 0维Matrix
29.5 解线性方程(Solving Linear Equations)
29.5.1 经典Gauss 消元法
29.5.2 主元(Pivoting)
29.5.3 测试(Testing)
29.5.4 融合操作(Fused Operations)
29.6 建议(Advice)
29.1 引言
孤立的语言特性是枯燥无味且毫无用处的。本章将介绍如何组合使用特性来解决一项具有挑战性的设计任务：通用 N 维矩阵。
我从未见过完美的矩阵类。事实上，鉴于矩阵的用途十分广泛，是否存在这样的矩阵类值得怀疑。在这里，我将介绍编写简单的 N 维稠密矩阵所需的编程和设计技术。至少，这个Matrix比程序员有时间直接使用vector或内置数组编写的任何东西都更容易使用，并且紧凑而快速。用于Matrix的设计和编程技术具有广泛的适用性。
29.1.1 基本的 Matrix 用法
Matrix <T,N> 是某个值类型的 N 维矩阵。可以这样使用它：
Matrix<double ,0> m0 {1}; // 0 维: 一个标量
Matrix<double ,1> m1 {1,2,3,4}; // 一维:一个向量 (4 个元素)
Matrix<double ,2> m2 { // 二维 (4*3 个元素)
{00,01,02,03}, // 0 行
{10,11,12,13}, // 1 行
{20,21,22,23} // 2 行
};
Matrix<double ,3> m3(4,7,9); // 三维 (4*7*9 个元素), 所有初始化为 0
Matrix<complex<double>,17> m17; // 17 维(目前无元素)
元素类型必须是我们可以存储的。我们不需要每个元素类型都具有浮点数的所有属性。例如：
Matrix<double ,2> md; //OK
Matrix<string,2> ms; // OK: 只是别试图进行算术运算
Matrix<Matrix<int,2>,2> mm { // 2×2 矩阵的 3×2 矩阵
// 矩阵是一个合理的“数字”
{ // row 0
{{1, 2}, {3, 4}}, // col 0
{{4, 5}, {6, 7}}, // col 1
},
{ // row 1
{{8, 9}, {0, 1}}, // col 0
{{2, 3}, {4, 5}}, // col 1
},
{ // row 2
{{1, 2}, {3, 4}}, // col 0
{{4, 5}, {6, 7}}, // col 1
}
};
矩阵运算并不具备与整数或浮点运算完全相同的数学性质(例如，矩阵乘法不满足交换律)，因此我们必须小心使用这种矩阵。
至于向量，我们使用 () 来指定大小，使用 {} 来指定元素值(§17.3.2.1,§17.3.4.1)。行数必须与指定的维数匹配，并且每个维度(每列)中的元素数必须匹配。例如：
Matrix<char,2> mc1(2,3,4); // 错 : 维数过多
Matrix<char,2> mc2 {
{'1','2','3'} //错: 第二维丢失初始化器
};
Matrix<char,2> mc2 {
{'1','2','3'},
{'4','5'} //错: 第三列丢失元素
};
Matrix<T,N> 的维度数(其 order())被指定为模板参数(此处为 N )。每个维度都有从初始化列表推导出来的元素数(其 extreme())，或者使用 () 符号指定为 Matrix 构造函数参数。元素总数称为 size()。例如：
Matrix<double ,1> m1(100); // 一维: 一个向量 (100 个元素)
Matrix<double ,2> m2(50,6000); //二维: 50*6000 个元素
auto d1 = m1.order(); // 1
auto d2 = m2.order(); // 2
auto e1 = m1.extent(0); // 100
auto e1 = m1.extent(1); // 错: m1 是一个一维
auto e2 = m2.extent(0); // 50
auto e2 = m2.extent(1); // 6000
auto s1 = m1.size(); // 100
auto s2 = m2.size(); // 50*6000
我们可以通过几种下标方式来访问 Matrix 的元素。例如：
Matrix<double ,2> m { //二维 (4*3 个元素)
{00,01,02,03}, // 0 行
{10,11,12,13}, // 1 行
{20,21,22,23} // 2 行
};
double d1 = m(1,2); // d==12
double d2 = m[1][2]; // d==12
Matrix<double ,1> m1 = m[1]; // 1 行 : {10,11,12,13}
double d3 = m1[2]; // d==12
我们可以定义一个输出函数用于调试，如下所示：
template<typename M>
Enable_if<Matrix_type<M>(),ostream&>
operator<<(ostream& os, const M& m)
{
os << '{';
for (siz e_t i = 0; i!=rows(m); ++i) {
os << m[i];
if (i+1!=rows(m)) os << ',';
}
return os << '}';
}
这里，Matrix_type 是一个概念(§24.3)。Enable_if 是 enable_if 类型的别名(§28.4)，因此此 operator<<() 返回一个 ostream&。
鉴于此，cout<<m 打印：{{0,1,2,3},{10,11,12,13},{20,21,22,23}}。
29.1.2 Matrix 的要求
在继续实施之前，请考虑我们可能希望拥有哪些属性：
• N 维，其中 N 是一个参数，其值可以从 0 到多个，无需为每个维度编写专门的代码。
• N 维存储通常很有用，因此元素类型可以是任何我们可以存储的元素(例如向量元素)。
• 数学运算应适用于任何可以合理描述为数字的类型，包括矩阵。
• Fortran 样式的下标，每个维度使用一个索引，例如，3-D 矩阵的 m(1,2,3) 产生一个元素。
• C 样式的下标，例如，m[7] 产生一行(一行是 N-D 矩阵的 N-1-D 子矩阵)。
• 下标应尽可能快速，并可能进行范围检查。
• 移动赋值和移动构造函数以确保有效传递矩阵结果并消除昂贵的临时变量。
• 一些数学矩阵运算，例如 + 和 ∗=。
• 一种读取、写入和传递子矩阵引用的方法，Matrix_refs，用于读取和写入元素。
• 基本保证形式下不存在资源泄漏(§13.2)。
• 融合关键操作，例如，m∗v+v2 作为单个函数调用。
这是一个相对较长且雄心勃勃的清单，但它并不等于“为所有人提供一切”。例如，我没有列出：
• 更多数学矩阵运算
• 专用矩阵(例如对角矩阵和三角矩阵)
• 稀疏矩阵支持
• 支持并行执行矩阵操作
无论这些属性多么有价值，它们都超出了提供基本编程技术所需的范围。
为了实现这一点，我结合使用了几种语言特性和编程技术：
• 类(当然)
• 使用数字和类型进行参数化
• 移动构造函数和赋值(以最小化复制)
• RAII(依赖于构造函数和析构函数)
• 可变参数模板(用于指定范围和索引)
• 初始化列表
• 运算符重载(以获得常规符号)
• 函数对象(用于携带有关下标的信息)
• 一些简单的超越板元编程(例如，用于检查初始化列表和区分 Matrix_refs 的读写)
• 实现继承以最小化代码复制。
显然，像这样的矩阵可以是内置类型(就像在许多语言中一样)，但这里的重点恰恰是它不是内置的，而是为用户提供自己制作的工具。
29.2 一个 Matrix 模板
为了给出一个概述，这里是 Matrix 的声明及其最有趣的操作：
template<typename T, size_t N>
class Matrix {
public:
static constexpr size_t order = N;
using value_type = T;
using iterator = typename std::vector<T>::iterator;
using const_iterator = typename std::vector<T>::const_iterator;
Matrix() = default;
Matrix(Matrix&&) = default; // 移动构造
Matrix& operator=(Matrix&&) = default;
Matrix(Matrix const&) = default; // 复制构造
Matrix& operator=(Matrix const&) = default;
˜Matrix() = default;
template<typename U>
Matrix(const Matrix_ref<U,N>&); //从 Matrix_ref 构造
template<typename U>
Matrix& operator=(const Matrix_ref<U,N>&); // 作 Matrix_ref 赋值
template<typename...Exts> //指定 extents
explicit Matrix(Exts... exts);
Matrix(Matrix_initializer<T,N>); //从 list 初始化
Matrix& operator=(Matrix_initializer<T,N>); // 从 list 赋值
template<typename U>
Matrix(initializer_list<U>) =delete; //除 template<typename U> 元素不使用 {}
Matrix& operator=(initializer_list<U>) = delete;
static constexpr size_t order() { return N; } // 维数
size_t extent(siz e_t n) const { return desc.extents[n]; } // 第 n 维中的元素
size_t siz e() const { return elems.size(); } // 元素总数
const Matrix_slice<N>& descriptor() const { return desc; } // 定义下标的分片
T∗ data() { return elems.data(); } // ‘‘flat’’ 元素访问
const T∗ data() const { return elems.data(); }
// ...
private:
Matrix_slice<N> desc; // 在 N 维中定义 extent 分片
vector<T> elems; // 元素
};
使用 vector<T> 保存元素可以让我们免于内存管理和异常安全的担忧。Matrix_slice 保存了以 N 维矩阵形式访问元素所需的大小(§29.4.2)。可以将其视为专门用于我们的 Matrix 的 gslice(§40.5.6)。
Matrix_ref (§29.4.3) 的行为与 Matrix 类似，只不过它引用的是 Matrix，通常是子 Matrix(例如行或列)，而不是拥有自己的元素。可以将其视为对子 Matrix 的引用。
Matrix_initializer<T,N> 是 Matrix<T,N> (§29.4.4) 的适当嵌套初始化器列表。
29.2.1 构造和赋值(Construction and Assignment)
默认的复制和移动操作具有恰到好处的语义：逐个复制或移动 desc(定义下标的分片描述符)和元素。请注意，对于元素存储的管理，Matrix 可从 vector 获得所有好处。同样，默认的构造函数和析构函数也具有恰到好处的语义。
采用范围(维度中元素的数量)的构造函数是可变参数模板的一个相当简单的例子(§28.6)：
template<typename T, siz e_t N>
template<typename... Exts>
Matrix<T,N>::Matrix(Exts... exts)
:desc{exts...}, // 复制 extents
elems(desc.size) // 分配 desc.size 个元素空间并默认初始它们
{ }
采用初始化列表的构造函数需要做一些工作：
template<typename T, siz e_t N>
Matrix<T, N>::Matrix(Matrix_initializ er<T,N> init)
{
Matrix_impl::derive_extents(init,desc.extents); //从初始化列表中推断范围 (§29.4.4)
elems.reserve(desc.size); //为 slices 保留空间
Matrix_impl::insert_flat(init,elems); //从初始化列表初始化(§29.4.4)
assert(elems.size() == desc.size);
}
Matrix_initializer 是一个适当嵌套的初始化列表 (§29.4.4)。范围由 Matrix_slice 构造函数推导、检查并存储在 desc 中。然后，元素由来自 Matrix_impl 命名空间的 insert_flat() 存储在 elems 中。
为了确保 {} 初始化仅用于元素列表，我删除了简单的 initializer_list 构造函数。这是为了强制使用 () 初始化范围。例如：
enum class Piece { none, cross, naught };
Matrix<Piece,2> board1 {
{Piece::none, Piece::none , Piece::none},
{Piece::none, Piece::none , Piece::none},
{Piece::none, Piece::none , Piece::cross}
};
Matrix<Piece,2> board2(3,3); // OK
Matrix<Piece,2> board3 {3,3}; // 错: 构造函数已从 initializer_list<int> 中删除
如果没有 =delete，最后一个定义就会被接受。
最后，我们必须能够从 Matrix_ref 构建，也就是说，从对 Matrix 或 Matrix 的一部分(子矩阵)的引用构建：
template<typename T, siz e_t N>
template<typename U>
Matrix<T,N>::Matrix(const Matrix_ref<U,N>& x)
:desc{x.desc}, elems{x.begin(),x.end()} // copy desc and elements
{
static_assert(Conver tible<U,T>(),"Matrix constructor: incompatible element types");
}
使用模板使我们能够从具有兼容元素类型的 Matrix 构建。
像往常一样，赋值类似于构造函数。例如：
template<typename T, siz e_t N>
template<typename U>
Matrix<T,N>& Matrix<T,N>::operator=(const Matrix_ref<U,N>& x)
{
static_assert(Conver tible<U,T>(),"Matrix =: incompatible element types");
desc = x.desc;
elems.assign(x.begin(),x.end());
return ∗this;
}
即，我们复制 Matrix 的成员。
29.2.2 下标和分片(Subscripting and Slicing)
可以通过下标(元素或行)、行和列或分片(行或列的一部分)来访问矩阵。
| Matrix<T,N>访问 | |
| m.row(i) | m的i行;一个Matrix_ref<T,N−1> |
| m.column(i) | m的i列; 一个Matrix_ref<T,N−1> |
| m[i] | C语言风格下标；m.row(i) |
| m(i,j) | Fortran 风格元素访问：m[i][j]; 一个T& 类型；下标数必须为 N |
| m(slice(i,n),slice(j)) | 具有分片的子矩阵：一个 Matrix_ref<T,N> ；slice(i,n)是下标的维数的元素 [i:i+n); slice(j)是下标的维数的元素 [i:max); max 是维数的程度(extent);下标的数量必须为N 。 |
这些是所有的成员函数：
template<typename T, siz e_t N>
class Matrix {
public:
// ...
template<typename... Args> //m(i,j,k) 按整数取下标
Enable_if<Matrix_impl::Requesting_element<Args...>(), T&>
operator()(Args... args);
template<typename... Args>
Enable_if<Matrix_impl::Requesting_element<Args...>(), const T&>
operator()(Args... args) const;
template<typename... Args> //m(s1,s2,s3) 按分片取下标
Enable_if<Matrix_impl::Requesting_slice<Args...>(), Matrix_ref<T, N>>
operator()(const Args&... args);
template<typename... Args>
Enable_if<Matrix_impl::Requesting_slice<Args...>(), Matrix_ref<const T,N>>
operator()(const Args&... args) const;
Matrix_ref<T,N−1> operator[](size_t i) { return row(i); } // m[i] 访问行
Matrix_ref<const T,N−1> operator[](size_t i) const { return row(i); }
Matrix_ref<T,N−1> row(siz e_tn); //访问行
Matrix_ref<const T,N−1> row(siz e_t n) const;
Matrix_ref<T,N−1> col(size_t n); // 访问列
Matrix_ref<const T,N−1> col(size_t n) const;
// ...
};
C 风格下标是通过 m[i] 选择并返回第 i 行来完成的：
template<typename T, siz e_t N>
Matrix_ref<T,N−1> Matrix<T,N>::operator[](siz e_t n)
{
return row(n); // §29.4.5
}
将 Matrix_ref (§29.4.3) 视为对子Matrix的引用。
Matrix_ref<T,0> 是特殊的，因此它引用单个元素(§29.4.6)。
Fortran 风格的下标是通过列出每个维数的索引来完成的，例如 m(i,j,k)，从而产生一个标量：
Matrix<int,2> m2 {
{01,02,03},
{11,12,13}
};
m(1,2) = 99; // 按 1 行 2 列重写元素; 即第 13 个元素
auto d1 = m(1); // 错: 下标太少
auto d2 = m(1,2,3); // 错: 下标太少
除了使用整数下标外，我们还可以使用分片下标。分片描述维数元素的子集(§40.5.4)。具体而言，分 片 {i,n} 指的是它适用的维度的元素 [i:i+n)。例如：
Matrix<int> m2 {
{01,02,03},
{11,12,13},
{21,22,23}
};
auto m22 = m(slice{1,2},slice{0,3});
现在 m22 是一个 Matrix<int,2>，其值为
{
{11,12,13},
{21,22,23}
}
第一个(行)下标slice{1,2} 选择最后两行，第二个(列)下标slice{0,3} 选择列中的所有元素。
带有分片下标的 () 的返回类型是 Matrix_ref，因此我们可以将其用作赋值的目标。例如：
m(slice{1,2},slice{0,3}) = {
{111,112,113},
{121,122,123}
}
现在，m具有值：
{
{01,02,03},
{111,112,113},
{121,122,123}
}
从某个点开始选择所有元素非常常见，因此有一个简写：slice{i} 表示slice{i,max}，其中 max 大于维度中的最大下标。因此，我们可以将 m(slice{1,2},slice{0,3}) 简化为等价的 m(slice{1,2},slice{0})。
另一个简单的常见情况是选择单行或单列的所有元素，因此一组分片下标中的纯整数下标 i 解释为分片{ i,1}。例如：
Matrix<int> m3 {
{01,02,03},
{11,12,13},
{21,22,23}
};
auto m31 = m(slice{1,2},1); // m31 becomes {{12},{22}}
auto m32 = m(slice{1,2},0); // m33 becomes {{11},{21}}
auto x = m(1,2); // x == 13
基本上所有用于数字编程的语言都支持分片下标的概念，所以希望它不会太陌生。
row()，column() 和 operator()() 的实现在 §29.4.5 中介绍。这些函数的 const 版本的实现与它们的非 const 版本的实现基本相同。关键区别在于 const 版本返回带有 const 元素的结果。
29.3 Matrix算术运算(Matrix Arithmetic Operations)
因此，我们可以创建矩阵、复制它们、访问它们的元素和行。但是，我们通常希望有数学运算，这样我们就不必通过访问单个元素(标量)来表达我们的算法。例如：
Matrix<int,2> mi {{1,2,3}, {4,5,6 }}; // 2-乘以-3
Matrix<int,2> m2 {mi}; // 复制
mi∗=2; //scale: {{2,4,6},{8,10,12}}
Matrix<int,2> m3 = mi+m2; // add: {{3,6,9},{12,15,18}}
Matrix<int,2> m4 {{1,2}, {3,4}, {5,6}}; // 3-by-2
Matrix<int,1> v = mi∗m4; //乘以: {{18,24,30},{38,52,66},{58,80,102}}
数学运算定义如下：
template<typename T, siz e_t N>
class Matrix {
// ...
template<typename F>
Matrix& apply(F f); // f(x) for every element x
template<typename M, typename F>
Matrix& apply(const M& m, F f); // 相应元素的 f(x,mx)
Matrix& operator=(const T& value); //标量赋值
Matrix& operator+=(const T& value); // 标量加
Matrix& operator−=(const T& value); // 标量减
Matrix& operator∗=(const T& value); // 标量乘
Matrix& operator/=(const T& value); // 标量除
Matrix& operator%=(const T& value); // 标量取模
template<typename M> // 矩阵加
Matrix& operator+=(const M& x);
template<typename M> // 矩阵减
Matrix& operator−=(const M& x);
// ...
};
// 二进制 +, -, * 作为非数函数提供
29.3.1 标量运算(Scalar Operations)
标量算术运算只是将其运算和右侧操作数应用于每个元素。例如：
template<typename T, size_t N>
Matrix<T,N>& Matrix<T,N>::operator+=(const T& val)
{
return apply([&](T& a) { a+=val; } ); // using a lambda (§11.4)
}
此 apply() 将函数(或函数对象)应用于其矩阵的每个元素：
template<typename T, siz e_t N>
template<typename F>
Matrix<T,N>& Matrix<T,N>::apply(F f)
{
for (auto& x : elems) f(x); // 这个循环使用了步进迭代器
return ∗this;
}
像往常一样，返回 ∗this 可以实现链接。例如：
m.apply(abs).apply(sqr t); // 对于所有的 i m[i] = sqrt(abs(m[i]))
与往常一样(§3.2.1.1,§18.3)，我们可以使用赋值运算符(例如 +=)在类外部定义“普通运算符”，例如 +。例如：
template<typename T, siz e_t N>
Matrix<T,N> operator+(const Matrix<T,N>& m, const T& val)
{
Matrix<T,N> res = m;
res+=val;
return res;
}
如果没有移动构造函数，这种返回类型将是一个严重的性能错误。
29.3.2 加法运算(Addition)
两个矩阵的加法与标量版本的加法非常相似：
template<typename T, size_t N>
template<typename M>
Enable_if<Matrix_type<M>(),Matrix<T,N>&> Matrix<T,N>::operator+=(const M& m)
{
static_assert(m.order()==N,"+=: mismatched Matrix dimensions");
assert(same_extents(desc,m.descriptor())); // 确保大小匹配
return apply(m, [](T& a,Value_type<M>&b) { a+=b; });
}
Matrix::apply(m,f) 是 Matrix::apply(f) 的双参数版本。它将 f 应用于其两个矩阵(m 和 ∗this)：
template<typename T, siz e_t N>
template<typename M, typename F>
Enable_if<Matrix_type<M>(),Matrix<T,N>&> Matrix<T,N>::apply(M& m, F f)
{
assert(same_extents(desc,m.descriptor())); // 确保类型匹配
for (auto i = begin(), j = m.begin(); i!=end(); ++i, ++j)
f(∗i,∗j);
return ∗this;
}
现在，operator+() 很容易定义：
template<typename T, size_t N>
Matrix<T,N> operator+(const Matrix<T,N>& a, const Matrix<T,N>& b)
{
Matrix<T,N> res = a;
res+=b;
return res;
}
这为两个相同类型的 Matrix 定义了 +，产生该类型的结果。我们应该概括一下：
template<typename T, typename T2, size_t N,
typename RT = Matrix<Common_type<Value_type<T>,Value_type<T2>>,N>
Matrix<RT,N> operator+(const Matrix<T,N>& a, const Matrix<T2,N>& b)
{
Matrix<RT,N> res = a;
res+=b;
return res;
}
如果 T 和 T2 是同一类型，则 Common_type 就是该类型。Common_type 类型函数派生自 std::common_type (§35.4.2)。对于内置类型，它(如 ?:)给出一种最能保留算术运算值的类型。如果对于我们想要组合使用的一对类型没有定义 Common_type，我们可以定义它。例如：
template<>
struct common_type<Quad,long double> {
using type = Quad;
};
现在 Common_type<Quad,long double> 是 Quad。
我们还需要涉及 Matrix_refs (§29.4.3) 的操作。例如：
template<typename T, siz e_t N>
Matrix<T,N> operator+(const Matrix_ref<T,N>& x, const T& n)
{
Matrix<T,N> res = x;
res+=n;
return res;
}
此类操作与 Matrix 操作完全相同。Matrix 和 Matrix_ref 元素访问之间没有区别：Matrix 和 Matrix_ref 之间的区别在于元素的初始化和所有权。
标量的减法、乘法等以及 Matrix_ref 的处理只是加法技术的重复。
29.3.3 乘法(Multiplication)
矩阵乘法并不像加法那么简单：N 乘 M 矩阵和 M 乘 P 矩阵的乘积是 N 乘 P 矩阵。对于 M==1，我们得到两个向量的乘积是一个矩阵，而对于 P==1，我们得到矩阵和向量的乘积是一个向量。我们可以将矩阵乘法推广到更高维度，但要做到这一点，我们必须引入张量 [Kolecki,2002]，我不想将关于编程技术和如何使用语言特性的讨论转移到物理和工程数学课上。所以，我将坚持一维和二维。
将一个 Matrix<T,1> 视为 N×1 矩阵，将另一个 Matrix<T,1> 视为 1×M 矩阵，我们得到：
template<typename T>
Matrix<T,2> operator∗(const Matrix<T,1>& u, const Matrix<T,1>& v)
{
const size_t n = u.extent(0);
const size_t m = v.extent(0);
Matrix<T,2> res(n,m); // an n-by-m matrix
for (siz e_t i = 0; i!=n; ++i)
for (siz e_t j = 0; j!=m; ++j)
res(i,j) = u[i]∗v[j];
return res;
}
这是最简单的情况：矩阵元素 res(i,j) 是 u[i]∗v[j]。我没有尝试推广以处理向量元素类型不同的情况。如有必要，可以使用讨论过的加法技巧。
请注意，我对 res 的每个元素都进行了两次写入：一次初始化为 T{}，一次赋值 u[i]∗v[j]。这大约使乘法的成本翻倍。如果这让您感到困扰，请编写一个没有这种开销的乘法，看看您的程序中是否存在差异。
接下来，我们可以将一个 N×M 矩阵与一个视为 M×1 矩阵的向量相乘。结果是一个 N×1 矩阵：
template<typename T>
Matrix<T,1> operator∗(const Matrix<T,2>& m, const Matrix<T,1>& v)
{
assert(m.extent(1)==v.extent(0));
const size_t n = m.extent(0);
Matrix<T,1> res(n);
for (siz e_t i = 0; i!=n; ++i)
for (siz e_t j = 0; j!=n; ++j)
res(i) += m(i,j)∗v(j);
return res;
}
请注意，res 的声明将其元素初始化为 T{}，对于数字类型为零，因此 += 从零开始。
N×M 矩阵乘以 M×P 矩阵的处理方式类似：
template<typename T>
Matrix<T,2> operator∗(const Matrix<T,2>& m1, const Matrix<T,2>& m2)
{
const size_t n = m1.extent(0);
const size_t m = m1.extent(1);
assert(m==m2.extent(0)); // 列和行必须匹配
const size_t p = m2.extent(1);
Matrix<T,2> res(n,p);
for (siz e_t i = 0; i!=n; ++i)
for (siz e_t j = 0; j!=m; ++j)
for (siz e_t k = 0; k!=p; ++k)
res(i,j) = m1(i,k)∗m2(k,j);
return res;
}
有很多方法可以优化这一重要操作。
最内层的循环可以更优雅地表达为：
res(i,j) = dot_product(m1[i],m2.column(j))
这里，dot_product() 只是标准库 inner_product() (§40.6.2) 的一个接口：
template<typename T>
T dot_product(const Matrix_ref<T,1>& a, const Matrix_ref<T,1>& b)
{
return inner_product(a.begin(),a.end(),b.begin(),0.0);
}
29.4 Matrix实现(Matrix Implementation)
到目前为止，我推迟了对 Matrix 实现中最复杂(对某些程序员来说也是最有趣)的“机械”部分的介绍。例如：什么是 Matrix_ref？什么是 Matrix_slice？如何从一串初始化列表初始化一个 Matrix 并确保其维度合理？我们如何确保不使用不合适的元素类型来实例化 Matrix？
呈现此代码的最简单方法是将 Matrix 的所有内容放在头文件中。在这种情况下，添加内联每个非成员函数的定义。
不是 Matrix ,Matrix_ref ,Matrix_slice 成员或通用接口一部分的函数定义放在命名空间 Matrix_impl 中。
29.4.1 slice()
用于分片下标的简单分片根据三个值描述从整数(下标)到元素位置(索引)的映射：
struct slice {
slice() :start(−1), length(−1), stride(1) { }
explicit slice(size_t s) :start(s), length(−1), stride(1) { }
slice(size_t s, size_t l, size_t n = 1) :start(s), length(l), stride(n) { }
size_t operator()(siz e_t i) const { return start+i∗stride; }
static slice all;
size_t star t; // 第一个索引
size_t length; // 包含索引数(可用于范围验证)
size_t stride; // 元素在序列中的距离
};
有一个标准库版本的分片；有关更详细的讨论，请参阅 §40.5.4。此版本提供了符号便利(例如，构造函数提供的默认值)。
29.4.2 Matrix分片(Matrix Slices)
Matrix_slice 是 Matrix 实现的一部分，它将一组下标映射到元素的位置。它使用广义分片的思想(§40.5.6)：
template<size_t N>
struct Matrix_slice {
Matrix_slice() = default; // 一个空矩阵: 无元素
Matrix_slice(size_t s, initializer_list<siz e_t> exts); // extents
Matrix_slice(size_t s, initializer_list<siz e_t> exts, initializer_list<siz e_t> strs);// extents 和 strides
template<typename...Dims> //N extents
Matrix_slice(Dims... dims);
template<typename... Dims,
typename = Enable_if<All(Conver tible<Dims,siz e_t>()...)>>
size_t operator()(Dims... dims) const; // 基于下标集计算索引
size_t size; // 元素总数
size_t start; // 起始偏移
array<siz e_t,N> extents; // 每一个维数中元素数目
array<siz e_t,N> strides; // 每一个维数中元素之间的偏移
};
换言之，Matrix_slice 描述了内存区域中的行和列。在通常的 C/C++ 行主序矩阵布局中，行的元素是连续的，而列的元素由固定数量的元素(步长)分隔。Matrix_slice 是一个函数对象，其运算符 ()() 进行步长计算(§40.5.6)：
template<size_t N>
template<typename... Dims>
size_t Matrix_slice<N>::operator()(Dims... dims) const
{
static_assert(sizeof...(Dims) == N, "");
size_t args[N] { size_t(dims)... }; // 复制参数进数组
return inner_product(args,args+N,strides.begin(),siz e_t(0));
}
下标必须高效。这是一个需要优化的简化算法。如果没有其他选择，可以使用特化来消除可变参数模板的参数包中的下标简化副本。例如：
template<>
struct Matrix_slice<1> {
// ...
size_t operator()(siz e_t i) const
{
return i;
}
}
template<>
struct Matrix_slice<2> {
// ...
size_t operator()(size_t i, size_t j) const
{
return i∗stides[0]+j;
}
}
Matrix_slice 是定义矩阵形状(其范围)和实现 N 维下标的基础。但是，它对于定义子矩阵也很有用。
29.4.3 Matrix_ref
Matrix_ref 基本上是用于表示子矩阵的 Matrix 类的克隆。但是，Matrix_ref 并不拥有其元素。它由 Matrix_slice 和指向元素的指针构成：
template<typename T, siz e_t N>
class Matrix_ref {
public:
Matrix_ref(const Matrix_slice<N>& s, T∗ p) :desc{s}, ptr{p} {}
// ... mostly like Matr ix ...
private:
Matrix_slice<N> desc; // 矩阵形状
T∗ptr; //矩阵第一个元素
};
Matrix_ref 仅指向“其”Matrix 的元素。显然，Matrix_ref 的寿命不应长于其 Matrix。例如：
Matrix_ref<double ,1> user()
{
Matrix<double ,2> m = {{1,2}, {3,4}, {5,6}};
return m.row(1);
}
auto mr = user(); // 麻烦
Matrix 和 Matrix_ref 之间非常相似，因此容易产生重复。如果这造成麻烦，我们可以从一个共同的基础中得出两者：
template<typename T, size_t N>
class Matrix_base {
// ... common stuff ...
};
template<typename T, siz e_t N>
class Matrix : public Matrix_base<T,N> {
// ... Matrix 特有...
private:
Matrix_slice<N> desc; // 矩阵形状
vector<T> elements;
};
template<typename T, size_t N>
class Matrix_ref : public Matrix_base<T,N> {
// ... Matrix_ref 特有 ...
private:
Matrix_slice<N> desc; // 矩阵形状
T∗ ptr;
};
29.4.4 Matrix 列表初始化
从 initializer_list 构造的 Matrix 构造函数以别名 Matrix_initializer 作为其参数类型：
template<typename T, size_t N>
using Matrix_initializer = typename Matrix_impl::Matrix_init<T, N>::type;
Matrix_init 描述了嵌套的 initializer_list 的结构。
Matrix_init<T,N> 仅具有 Matrix_init<T,N−1> 作为其成员类型：
template<typename T, siz e_t N>
struct Matrix_init {
using type = initializer_list<typename Matrix_init<T,N−1>::type>;
};
N==1 比较特殊。这就是我们得到(嵌套最深的)initializer_list<T> 的地方：
template<typename T>
struct Matrix_init<T,1> {
using type = initializer_list<T>;
};
为了避免意外，我们将 N=0 定义为错误：
template<typename T>
struct Matrix_init<T,0>; // 未基于目的定义
现在我们可以完成采用 Matrix_initializer 的 Matrix 构造函数：
template<typename T, size_t N>
Matrix<T, N>::Matrix(Matrix_initializer<T,N> init)
{
Matrix_impl::derive_extents(init,desc.extents); //从initializer list推导 extents (§29.4.4)
elems.reserve(desc.size); // 为分片创建空间
Matrix_impl::insert_flat(init,elems); // 从初始化器列表初始化(§29.4.4)
assert(elems.size() == desc.size);
}
为此，我们需要两个操作来对 Matrix<T,N> 的初始化列表树进行递归：
• derive_extents() 确定矩阵的形状：
• 检查树是否真的有 N 深
• 检查每行(子初始化列表)是否具有相同数量的元素
• 设置每行的范围
• insert_flat() 将初始化列表<T>树的元素复制到矩阵的元素中。
从 Matrix 构造函数调用 derived_extents() 来初始化其 desc，如下所示：
template<size_t N, typename List>
array<size_t, N> derive_extents(const List& list)
{
array<siz e_t,N> a;
auto f = a.begin();
add_extents<N>(f,list); // 将 extents 从 list 放入 f[]
return a;
}
你给它一个initializer_list，它会返回一个范围数组。
递归从 N 到最后的 1，其中 initializer_list 是 initializer_list<T>。
template<size_t N, typename I, typename List>
Enable_if<(N>1),void> add_extents(I& first, const List& list)
{
assert(check_non_jagged(list));
∗first = list.size();
add_extents<N−1>(++first,∗list.begin());
}
template<size_t N, typename I, typename List>
Enable_if<(N==1),void> add_extents(I& first, const List& list)
{
∗first++ = list.size(); // 我们已达最深嵌套
}
check_non_jagged() 函数检查所有行是否具有相同数量的元素：
template<typename List>
bool check_non_jagged(const List& list)
{
auto i = list.begin();
for (auto j = i+1; j!=list.end(); ++j)
if (i−>size()!=j−>siz e())
return false;
return true;
}
我们需要 insert_flat() 来获取可能嵌套的初始化列表，并将其元素作为 vector<T> 呈现给 Matrix<T>。它将提供给 Matrix 的初始化列表作为 Matrix_initializer，并提供元素作为目标：
template<typename T, typename Vec>
void insert_flat(initializer_list<T> list, Vec& vec)
{
add_list(list.begin(),list.end(),vec);
}
遗憾的是，我们不能依赖在内存中连续分配的元素，所以我们需要通过一组递归调用来构建向量。如果我们有一个initializer_list列表，我们会递归遍历每个列表：
template<typename T, typename Vec> // 嵌入 initializer_list
void add_list(const initializer_list<T>∗ first, const initializer_list<T>∗ last, Vec& vec)
{
for (;first!=last;++first)
add_list(first−>begin(),first−>end(),vec);
}
当我们到达一个包含非 initializer_list 元素的列表时，我们将这些元素插入到我们的vector 中：
template<typename T, typename Vec>
void add_list(const T∗ first, const T∗ last, Vec& vec)
{
vec.insert(vec.end(),first,last);
}
我使用 vec.insert(vec.end(),first,last) 因为没有接受序列参数的 push_back()。
29.4.5 Matrix 访问
矩阵可通过行、列、分片 (§29.4.1) 和元素 (§29.4.3) 进行访问。row() 或 column() 操作返回 Matrix_ref<T,N−1>，带整数的 () 下标操作返回 T&，带分片的 () 下标操作返回 Matrix<T,N>。
只要 1 < N ，Matrix<T,N> 的行就是 Matrix_ref<T,N−1>：
template<typename T, siz e_t N>
Matrix_ref<T,N−1> Matrix<T,N>::row(siz e_t n)
{
assert(n<rows());
Matrix_slice<N−1> row;
Matrix_impl::slice_dim<0>(n,desc,row);
return {row,data()};
}
我们需要对 N==1 和 N==0 进行特化：
template<typename T>
T& Matrix<T,1>::row(siz e_t i)
{
return &elems[i];
}
template<typename T>
T& Matrix<T,0>::row(siz e_t n) = delete;
选择列 () 本质上与选择行()相同。区别仅在于 Matrix_slice 的构造：
template<typename T, siz e_t N>
Matrix_ref<T,N−1> Matrix<T,N>::column(siz e_t n)
{
assert(n<cols());
Matrix_slice<N−1> col;
Matrix_impl::slice_dim<1>(n,desc,col);
return {col,data()};
}
const 版本是等价的。
Requesting_element() 和 Requesting_slice() 是一组整数的概念，分别用于用一组整数取下标和用分片取下标(§29.4.5)。它们检查一系列访问函数参数是否具有适合用作下标的类型。
用整数取下标像这样定义：
template<typename T, siz e_tN> //subscripting with integers
template<typename... Args>
Enable_if<Matrix_impl::Requesting_element<Args...>(),T&>
Matrix<T,N>::operator()(Args... args)
{
assert(Matrix_impl::check_bounds(desc, args...));
return ∗(data() + desc(args...));
}
check_bounds() 谓词检查下标的数量是否等于维数，以及下标是否在界限内：
template<size_t N, typename... Dims>
bool check_bounds(const Matrix_slice<N>& slice, Dims... dims)
{
size_t indexes[N] {size_t(dims)...};
return equal(indexes, indexes+N, slice.extents, less<size_t> {});
}
矩阵中元素的实际位置是通过调用矩阵的 Matrix_slice 的广义分片计算(以函数对象的形式呈现)来计算的：desc(args...)。将其添加到数据(data())的开头，我们就得到了我们的位置：
return ∗(data() + desc(args...));
这样，声明中最神秘的部分就留到最后了。operator()() 的返回类型的规范如下所示：
Enable_if<Matrix_impl::Requesting_element<Args...>(),T&>
因此返回类型为 T&，前提是
Matrix_impl::Requesting_element<Args...>()
为 true(§28.4)。此谓词仅检查每个下标是否都可以通过使用标准库谓词 is_convertible(§35.4.1)的概念版本转换为所需的 size_t：
template<typename... Args>
constexpr bool Requesting_element()
{
return All(Conver tible<Args,siz e_t>()...);
}
All() 只是将其谓词应用于可变模板的每个元素：
constexpr bool All() { return true; }
template<typename... Args>
constexpr bool All(bool b, Args... args)
{
return b && All(args...);
}
在 Request 中使用谓词 (Requesting_element) 和 Enable_if() ‘‘hidden’’ 的原因是为了在元素和分片下标运算符之间进行选择。分片下标运算符使用的谓词如下所示：
template<typename... Args>
constexpr bool Requesting_slice()
{
return All((Conver tible<Argssiz e_t>() || Same<Args,slice>())...)
&& Some(Same<Args,slice>()...);
}
也就是说，如果至少有一个切片参数，并且所有参数都可以转换为切片或 size_t，那么我们就有一些东西可以用来描述 Matrix<T,N> ：
template<typename T, siz e_t N> // 分片取下标
template<typename... Args>
Enable_if<Matrix_impl::Requesting_slice<Args...>(), Matrix_ref<T,N>>
Matrix<T,N>::operator()(const Args&... args)
{
matrix_slice<N> d;
d.start = matrix_impl::do_slice(desc,d,args...);
return {d,data()};
}
以 Matrix_slice 中的范围和步幅表示并用于分片下标的分片按以下方式计算：
template<size_t N, typename T, typename ... Args>
size_t do_slice(const Matrix_slice<N>& os, Matrix_slice<N>& ns, const T& s, const Args&... args)
{
size_t m = do_slice_dim<siz eof...(Args)+1>(os,ns,s);
size_t n = do_slice(os,ns,args...);
return m+n;
}
与往常一样，递归通过一个简单的函数终止：
template<size_t N>
size_t do_slice(const Matrix_slice<N>& os, Matrix_slice<N>& ns)
{
return 0;
}
do_slice_dim() 是一个比较棘手的计算(为了得到正确的分片值)，但是并没有展示出新的编程技术。
29.4.6 0维Matrix
矩阵代码包含大量 N-1 的出现，其中 N 是维数。因此，N==0 很容易成为一个令人讨厌的特殊情况(对于编程和数学而言)。在这里，我们通过定义一个特化来解决问题：
template<typename T>
class Matrix<T,0> {
public:
static constexpr size_t order = 0;
using value_type = T;
Matrix(const T& x) : elem(x) { }
Matrix& operator=(const T& value) { elem = value; return ∗this; }
T& operator()() { return elem; }
const T& operator()() const { return elem; }
operator T&() { return elem; }
operator const T&() { return elem; }
private:
T elem;
};
Matrix<T,0> 并不是真正的矩阵。它存储单个 T 类型的元素，并且只能转换为该类型的引用。
29.5 解线性方程(Solving Linear Equations)
如果你理解要解的问题和用于表达解的数学知识，那么数值计算的代码就有意义，如果你不理解，则代码往往看起来完全是胡说八道。如果你学过基本的线性代数，这里使用的示例应该相当简单；如果没有，只需将其视为将教科书求解转录为代码的示例即可，而无需进行太多改写。
这里选择的示例是为了展示 Matrixes 的合理现实和重要用途。我们将求解一组(任意一组)这种形式的线性方程：
在这里，x 表示 n 个未知数；a 和 b 是给定的常数。为简单起见，我们假设未知数和常数都是浮点值。目标是求得同时满足 n 个方程的未知数的值。这些方程可以用一个矩阵和两个向量简洁地表示：
Ax = b
这里，A 是由以下系数定义的 n×n 方阵：
。
向量 x 和 b 分别是未知数和常数向量：
和
。
该系统可能有零个、一个或无数个解，具体取决于矩阵 A 和向量 b 的系数。有各种方法可以解线性方程组。我们使用一种经典方案，称为 Gauss 消元法 [Freeman,1992]、[Stewart,1998]、[Wood,1999]。首先，我们变换 A 和 b，使 A 成为上三角矩阵。所谓“上三角”，是指 A 对角线下方的所有系数均为零。换句话说，系统如下所示：
这很容易做到。位置 a(i,j) 的零点是通过将行 i 的方程乘以一个常数而得到的，这样 a(i,j) 等于列 j 中的另一个元素，比如 a(k,j)。完成后，我们只需减去这两个方程，a(i,j)==0，行 i 中的其他值就会相应改变。
如果我们能让所有对角线系数都非零，那么该方程组就有唯一的解，可以通过“反向代入”求得。最后一个方程很容易解：
。
显然，x[n] 是 b[n]/a(n,n)。完成后，从系统中消除第 n 行，继续查找 x[n−1] 的值，依此类推，直到计算出 x[1] 的值。对于每个 n，我们除以 a(n,n)，因此对角线值必须为非零。如果不成立，则反向代换方法失败，这意味着系统有零个或无限多个解。
29.5.1 经典Gauss 消元法
现在让我们看一下 C++ 代码来表达这一点。首先，我们通过常规命名我们将使用的两种 Matrix 类型来简化我们的符号：
using Mat2d = Matrix<double ,2>;
using Vec = Matrix<double ,1>;
接下来，我们将表达我们期望的计算：
Vec classical_gaussian_elimination(Mat2d A, Vec b)
{
classical_elimination(A, b);
return back_substitution(A, b);
}
也就是说，我们复制输入 A 和 b (使用按值调用)，调用一个函数来解决系统，然后通过反向替换计算返回结果。关键是我们对问题的分解和解决方案的符号都来自教科书。为了完成我们的解决方案，我们必须实现 classic_elimination() 和 back_substitution()。同样，解决方案在教科书中：
void classical_elimination(Mat2d& A, Vec& b)
{
const size_t n = A.dim1();
// 从第一列到倒数第二列遍历，将对角线下的所有元素填充为零:
for (siz e_t j = 0; j!=n−1; ++j) {
const double pivot = A(j, j);
if (pivot==0) throw Elim_failure(j);
// 将第 i 行对角线下的每个元素填零:
for (siz e_t i = j+1; i!=n; ++i) {
const double mult = A(i,j) / pivot;
A[i](slice(j)) = scale_and_add(A[j](slice(j)), −mult,A[i](slice(j)));
b(i) −= mult∗b(j); // 对b进行相应的修改
}
}
}
主元(pivot)是位于我们当前正在处理的行对角线上的元素。它必须非零，因为我们需要除以它；如果它为零，我们通过抛出异常来放弃：
Vec back_substitution(const Mat2d& A, const Vec& b)
{
const size_t n = A.dim1();
Vec x(n);
for (siz e_t i = n−1; i>=0; −−i) {
double s = b(i)−dot_product(A[i](slice(i+1)),x(slice(i+1)));
if (double m = A(i,i))
x(i) = s/m;
else
throw Back_subst_failure(i);
}
return x;
}
29.5.2 主元(Pivoting)
我们可以通过对行进行排序，使零和小值远离对角线，从而避免除以零的问题，并实现更稳健的解决方案。“更稳健”是指对舍入误差不那么敏感。但是，随着我们将零放在对角线下，值会发生变化，因此我们还必须重新排序，以使小值远离对角线(也就是说，我们不能只对矩阵进行重新排序，然后使用经典算法)：
void elim_with_partial_pivot(Mat2d& A, Vec& b)
{
const size_t n = A.dim1();
for (siz e_t j = 0; j!=n; ++j) {
size_t pivot_row = j;
// look for a suitable pivot:
for (siz e_t k = j+1; k!=n; ++k)
if (abs(A(k,j)) > abs(A(pivot_row,j)))
pivot_row = k;
// swap the rows if we found a better pivot:
if (pivot_row!=j) {
A.swap_rows(j,pivot_row);
std::swap(b(j),b(pivot_row));
}
// elimination:
for (siz e_t i = j+1; i!=n; ++i) {
const double pivot = A(j,j);
if (pivot==0) error("can't solve: pivot==0");
const double mult = A(i,j)/pivot;
A[i].slice(j) = scale_and_add(A[j].slice(j), −mult, A[i].slice(j));
b(i) −= mult∗b(j);
}
}
}
我们使用 swap_rows() 和 scale_and_multiply() 来使代码更加常规，并避免编写显式循环。
29.5.3 测试(Testing)
显然，我们必须测试我们的代码。幸运的是，有一个简单的方法可以做到这一点：
void solve_random_system(size_t n)
{
Mat2d A = random_matrix(n); // 产生随机 Mat2d
Vec b = random_vector(n); // 产生随机 Vec
cout << "A = " << A << endl;
cout << "b = " << b << endl;
try {
Vec x = classical_gaussian_elimination(A, b);
cout << "classical elim solution is x = " << x << endl;
Vec v = A ∗ x;
cout << " A ∗ x = " << v << endl;
}
catch(const exception& e) {
cerr << e.what() << endl;
}
}
我们可以通过三种方式来获得 catch 子句：
• 代码中存在错误(但作为乐观主义者，我们认为不存在任何错误)
• 导致 classic_elimination() 出错的输入(使用 elim_with_partial_pivot() 可最大程度地降低发生这种情况的可能性)
• 舍入错误
然而，我们的测试并不像我们所希望的那样现实，因为真正的随机矩阵不太可能给classical_elimination()带来问题。
为了验证我们的解，我们打印出 A∗x，它最好等于 b(或者考虑到舍入误差，足够接近我们的目的)。舍入误差的可能性是我们不这样做的原因：
if (A∗x!=b) error("substitution failed");
因为浮点数只是实数的近似值，所以我们必须接受近似正确的答案。一般来说，最好避免在浮点计算的结果上使用 == 和 !=：浮点本质上是一种近似值。如果我觉得需要进行机器检查，我会定义一个 equal() 函数，并考虑哪些错误范围是可接受的，然后写成：
if (equal(A∗x,b)) error("substitution failed");
random_matrix() 和 random_vector() 是随机数的简单用法，留给读者作为简单的练习。
29.5.4 融合操作(Fused Operations)
除了提供高效的原始操作之外，通用矩阵类还必须处理三个相关问题以满足注重性能的用户：
[1] 必须尽量减少临时变量的数量。
[2] 必须尽量减少矩阵的复制。
[3] 必须尽量减少复合操作中对同一数据的多次循环。
考虑 U=M∗V+W，其中 U、V 和 W 是向量 (Matrix<T,1>)，M 是 Matrix<T,2>。一个简单的实现为 M∗V 和 M∗V+W 引入了临时向量，并复制了 M∗V 和 M∗V+W 的结果。一个聪明的实现调用一个函数 mul_add_and_assign(&U,&M,&V,&W)，该函数不引入任何临时向量。
移动构造函数有帮助：用于 M∗V 的临时变量用于 (M∗V)+W。如果我们写成
Matrix<double ,1> U=M∗V+W;
我们将消除所有元素副本：在 M∗V 中的局部变量中分配的元素是最终出现在 U 中的元素。
剩下的就是合并循环的问题：循环融合。对于多种表达式，这种程度的优化很少是必要的，因此解决效率问题的一个简单方法是提供诸如 mul_add_and_assign() 之类的函数，并让用户在重要的地方调用它们。但是，可以设计一个矩阵，以便自动将此类优化应用于正确形式的表达式。也就是说，我们可以将 U=M∗V+W 视为具有四个操作数的单个运算符的使用。基本技术已针对 ostream 操纵器进行了演示(§38.4.5.2)。一般来说，它可用于使 n 个二进制运算符的组合像 (n+1) 元运算符一样工作。处理 U=M∗V+W 需要引入两个辅助类。但是，通过启用更强大的优化技术，该技术可以在某些系统上带来令人印象深刻的加速(比如 30 倍)。首先，为简单起见，我们将自己限制为双精度浮点数的二维矩阵：
using Mat2d = Matrix<double ,2>;
using Vec = Matrix<double ,1>;
我们定义 Mat2d 与 Vec 相乘的结果：
struct MVmul {
const Mat2d& m;
const Vec& v;
MVmul(const Mat2d& mm, const Vec &vv) :m{mm}, v{vv} { }
operator Vec(); // 计算并返回结果
};
inline MVmul operator∗(const Mat2d& mm, const Vec& vv)
{
return MVmul(mm,vv);
}
这种“乘法”应该取代 §29.3 中的乘法，除了存储对其操作数的引用外，什么也不做；M∗V 的估算被推迟。∗ 产生的对象与许多技术社区中所谓的闭包密切相关。同样，我们可以处理如果添加 Vec 会发生什么：
struct MVmulVadd {
const Mat2d& m;
const Vec& v;
const Vec& v2;
MVmulVadd(const MVmul& mv, const Vec& vv) :m(mv.m), v(mv.v), v2(vv) { }
operator Vec(); // 计算并返回结果
};
inline MVmulVadd operator+(const MVmul& mv, const Vec& vv)
{
return MVmulVadd(mv,vv);
}
这推迟了 M∗V+W 的估算。现在我们必须确保在将其分配给 Vec 时使用良好的算法对其进行估算：
template<>
class Matrix<double ,1> { // 特化(仅为此例)
// ...
public:
Matrix(const MVmulVadd& m) // 通过 m 的结果进行初始化
{
// 分配元素内存, 等
mul_add_and_assign(this,&m.m,&m.v,&m.v2);
}
Matrix& operator=(const MVmulVadd& m) // 将 m 的结果赋给 *this
{
mul_add_and_assign(this,&m.m,&m.v,&m.v2);
return ∗this;
}
// ...
};
现在 U=M∗V+W 自动扩展为
U.operator=(MVmulVadd(MVmul(M,V),W))
由于内联，它解析为所需的简单调用
mul_add_and_assign(&U,&M,&V,&W)
显然，这消除了复制和临时变量。此外，我们可能会以优化的方式编写mul_add_and_assign()。但是，如果我们只是以一种相当简单和未优化的方式编写它，它仍然会以一种为优化器提供巨大机会的形式。
这种技术的重要性在于，大多数真正时间紧迫的向量和矩阵计算都是使用一些相对简单的语法形式完成的。通常，优化六个运算符的表达式并没有真正的好处。为此，我们通常无论如何都想编写一个函数。
该技术基于使用编译时分析和闭包对象将子表达式的求值转移到表示复合操作的对象的思想。它可以应用于各种问题，这些问题具有共同的属性，即在进行求值之前需要将几条信息收集到一个函数中。我将生成的用于延迟求值的对象称为组合闭包对象，或简称为合成器(compositors)。
如果使用这种组合技术来延迟所有操作的执行，则它被称为表达式模板 [Vandevoorde,2002] [Veldhuizen,1995]。表达式模板系统地使用函数对象将表达式表示为抽象语法树 (AST)。
29.6 建议(Advice)
[1] 列出基本用例；§29.1.1。
[2] 始终提供输入和输出操作以简化简单测试(例如，单元测试)；§29.1.1。
[3] 仔细列出程序、类或库理想情况下应具有的属性；§29.1.2。
[4] 列出程序、类或库中被认为超出项目范围的属性；§29.1.2。
[5] 设计容器模板时，请仔细考虑元素类型的要求；§29.1.2。
[6] 考虑设计如何适应运行时检查(例如，用于调试)；§29.1.2。
[7] 如果可能，设计一个类来模仿现有的专业符号和语义；§29.1.2。
[8] 确保设计不会泄漏资源（例如，每个资源都有唯一的所有者并使用 RAII）；§29.2。
[9] 考虑如何构造和复制类；§29.1.1。
[10] 提供完整、灵活、高效且语义上有意义的元素访问；§29.2.2、§29.3。
[11] 将实现细节放在自己的 _impl 命名空间中；§29.4。
[12] 提供不需要直接访问表示的通用操作作为辅助函数；§29.3.2、§29.3.3。
[13] 为了快速访问，请保持数据紧凑并使用访问器对象提供必要的非平凡访问操作；§29.4.1、§29.4.2、§29.4.3。
[14] 数据的结构通常可以表示为嵌套的初始化列表；§29.4.4。
[15] 处理数字时，始终要考虑“最终情况”，例如零和“许多”；§29.4.6。
[16] 除了单元测试和测试代码是否符合要求外，还要通过实际使用示例测试设计；§29.5。
[17] 考虑设计如何适应异常严格的性能要求；§29.5.4
内容来源：
<<The C++ Programming Language >> 第4版，作者 Bjarne Stroustrup
