Boost C++ 库

...世界上最受尊敬和专业设计的 C++ 库项目之一。 — Herb Sutter 和 Andrei Alexandrescu, C++ 编码标准

uBLAS 概述

理由

如果每种数值软件都可以用 C++ 编写而不会损失效率，那就太好了，但除非能找到在不损害 C++ 类型系统的情况下实现这一目标的方法，否则最好依赖 Fortran、汇编语言或特定于架构的扩展（Bjarne Stroustrup）。

这个 C++ 库的目标是科学计算，主要关注基本的线性代数结构，包括矩阵和向量，以及它们相应的抽象运算。主要的设计目标是

另一个目的是评估使用这种矩阵和向量类所导致的抽象开销是否可以接受。

资源

这个库的开发受到了几个类似项目的指导

BLAS 似乎是最广泛使用的基本线性代数结构库，因此可以称之为事实标准。它的接口是过程式的，各个函数在某种程度上是从简单的线性代数运算中抽象出来的。由于它使用了 Fortran 及其优化实现，它似乎也是最快的库之一。由于我们决定以面向对象的方式设计和实现我们的库，因此技术方法截然不同。然而，任何人都应该能够用我们的库运算符表达 BLAS 抽象，并比较实现的效率。

Blitz++ 是一个令人印象深刻的 C++ 库。它的主要设计似乎面向多维数组及其相关的运算符，包括张量。Blitz++ 的作者声明，由于他使用表达式模板和模板元编程的实现技术，该库的性能与相应的 Fortran 代码相当甚至更好。但是，我们看到了一些理由来开发我们自己的设计和实现方法。我们不知道是否有人尝试使用 Blitz++ 实现传统的线性代数和其他数值算法。我们还推测，即使在今天，Blitz++ 也需要最先进的 C++ 编译器技术，因为它使用了大量的语言特性。另一方面，Blitz++ 使我们确信，使用表达式模板对于将抽象开销降低到可接受的限度是强制性的。

POOMA 的设计目标在许多方面似乎与 Blitz++ 的并行。它扩展了 Blitz++ 的概念，加入了来自偏微分方程和理论物理领域的类。该实现甚至支持并行架构。

MTL 是另一种支持 C++ 中基本线性代数运算的方法。它的设计主要似乎受到 BLAS 和 C++ 标准模板库的影响。我们认同线性代数库必须提供与 BLAS 相当的功能的观点。另一方面，我们认为，C++ 标准库的概念尚未被证明能够支持所需的数值计算。另一个不同之处是，MTL 目前似乎没有使用表达式模板。这可能会导致以下两种后果之一：可能损失表达能力或可能损失性能。

概念

数学符号

数学符号的使用可以简化科学算法的开发。因此，一个仔细实现基本线性代数概念的 C++ 库应该重载矩阵和向量类上选定的 C++ 运算符。

我们决定为以下原语使用运算符重载

描述 运算符
向量和矩阵的索引 vector::operator(size_t i);
matrix::operator(size_t i, size_t j);
向量和矩阵的赋值 vector::operator = (const vector_expression &);
vector::operator += (const vector_expression &);
vector::operator -= (const vector_expression &);
vector::operator *= (const scalar_expression &);
matrix::operator = (const matrix_expression &);
matrix::operator += (const matrix_expression &);
matrix::operator -= (const matrix_expression &);
matrix::operator *= (const scalar_expression &);
向量和矩阵上的一元运算 vector_expression operator - (const vector_expression &);
matrix_expression operator - (const matrix_expression &);
向量和矩阵上的二元运算 vector_expression operator + (const vector_expression &, const vector_expression &);
vector_expression operator - (const vector_expression &, const vector_expression &);
matrix_expression operator + (const matrix_expression &, const matrix_expression &);
matrix_expression operator - (const matrix_expression &, const matrix_expression &);
向量和矩阵与标量的乘法 vector_expression operator * (const scalar_expression &, const vector_expression &);
vector_expression operator * (const vector_expression &, const scalar_expression &);
matrix_expression operator * (const scalar_expression &, const matrix_expression &);
matrix_expression operator * (const matrix_expression &, const scalar_expression &);

我们决定不为以下其他原语使用运算符重载

描述 函数
向量与矩阵的左乘 vector_expression prod<vector_type > (const matrix_expression &, const vector_expression &);
vector_expression prod (const matrix_expression &, const vector_expression &);
向量与矩阵的右乘 vector_expression prod<vector_type > (const vector_expression &, const matrix_expression &);
vector_expression prod (const vector_expression &, const matrix_expression &);
矩阵的乘法 matrix_expression prod<matrix_type > (const matrix_expression &, const matrix_expression &);
matrix_expression prod (const matrix_expression &, const matrix_expression &);
向量的内积 scalar_expression inner_prod (const vector_expression &, const vector_expression &);
向量的外积 matrix_expression outer_prod (const vector_expression &, const vector_expression &);
矩阵的转置 matrix_expression trans (const matrix_expression &);

效率

为了实现数值计算的效率目标，必须克服在 C++ 中制定抽象概念时的两个难题，即临时对象和虚函数调用。表达式模板解决了这些问题，但往往会减慢编译时间。

消除临时对象

向量和矩阵上的抽象公式通常由几个一元和二元运算组成。评估这种公式的传统方法是首先将组合的每个叶子运算评估为临时对象，然后评估组合结果为另一个临时对象。这种方法在时间方面成本很高，特别是对于小型向量和矩阵，在空间方面成本很高，特别是对于大型向量和矩阵。解决这个问题的方法是使用惰性求值，这在现代函数式编程语言中是已知的。这种方法的原则是逐元素评估复杂表达式，并将其直接分配给目标。

导致两个有趣且危险的事实

别名

在向量或矩阵上使用逐元素求值可能会产生严重的副作用。考虑矩阵向量积 x = A x。评估 A1x 并赋值给 x1 会更改右侧，从而导致 A2x 的评估返回错误的结果。在这种情况下，赋值的左侧和右侧都有元素 xn 的别名。

我们针对此问题的解决方案是将赋值的右侧评估为临时对象，然后将此临时对象赋值给左侧。为了允许进一步优化，我们为每个赋值运算符提供相应的成员函数，并提供 noalias 语法。 通过使用此语法，程序员可以确认赋值的左侧和右侧是独立的，因此逐元素求值和直接赋值给目标是安全的。

复杂度

在某些情况下，计算复杂度可能会出乎意料地高。考虑链式矩阵向量积 A (B x)。A (B x) 的传统评估是二次的。B xi 的延迟评估是线性的。由于每个元素 B xi 的需求线性地依赖于大小，因此链式矩阵向量积 A (B x) 的完全延迟评估是三次的。在这种情况下，需要在表达式中重新引入临时对象。

消除虚函数调用

惰性表达式求值通常会导致术语类层次结构的定义。这导致使用动态多态性来访问向量和矩阵的单个元素，这也已知在时间方面成本很高。几年前，David Vandervoorde 和 Todd Veldhuizen 独立地找到了一个解决方案，通常称为表达式模板。表达式模板包含惰性求值，并将动态多态性替换为静态多态性，即编译时多态性。表达式模板严重依赖著名的 Barton-Nackman 技巧，Jim Coplien 也将其称为“奇特的递归模板定义”。

表达式模板构成了我们实现的基础。

编译时间

众所周知，表达式模板对当前可用的编译器提出了挑战。我们通过始终如一地使用 Barton-Nackman 技巧，显著减少了所需表达式模板的数量。

我们还决定支持双重传统实现（即不使用表达式模板），其中包含向量和矩阵操作的广泛边界和类型检查，以支持开发周期。从调试模式切换到发布模式由 <cassert> 的 NDEBUG 预处理器符号控制。

功能性

每个支持线性代数的 C++ 库都将与长期存在的 Fortran 包 BLAS 进行比较。我们现在描述如何将 BLAS 调用映射到我们的类。

页面矩阵和向量运算概述简要总结了向量和矩阵上最常用的运算。

Blas Level 1

BLAS 调用 映射的库表达式 数学描述 注释
sasum 或 dasum norm_1 (x) sum |xi| 计算实向量的 l1（和）范数。
scasum 或 dzasum real (sum (v)) + imag (sum (v)) sum re(xi) + sum im(xi) 计算复向量元素的和。
_nrm2 norm_2 (x) sqrt (sum |xi|2 ) 计算向量的 l2（欧几里得）范数。
i_amax norm_inf (x)
index_norm_inf (x)		 max |xi| 计算向量的 linf（最大）范数。
BLAS 计算具有此值的第一个元素的索引。
_dot
_dotu
_dotc		 inner_prod (x, y)或
inner_prod (conj (x), y)		 xT y 或
xH y		 计算两个向量的内积。
BLAS 实现了某些循环展开。
dsdot
sdsdot		 a + prec_inner_prod (x, y) a + xT y 以双精度计算内积。
_copy x = y
y.assign (x)		 x <- y 将一个向量复制到另一个向量。
BLAS 实现了某些循环展开。
_swap swap (x, y) x <-> y 交换两个向量。
BLAS 实现了某些循环展开。
_scal
csscal
zdscal		 x *= a x <- a x 缩放一个向量。
BLAS 实现了某些循环展开。
_axpy y += a * x y <- a x + y 添加一个缩放的向量。
BLAS 实现了某些循环展开。
_rot
_rotm
csrot
zdrot		 t.assign (a * x + b * y),
y.assign (- b * x + a * y),
x.assign (t)		 (x, y) <- (a x + b y, -b x + a y) 应用平面旋转。
_rotg
_rotmg		   (a, b) <-
  (? a / sqrt (a2 + b2),
    ? b / sqrt (a2 + b2)) 或
(1, 0) <- (0, 0)		 构造平面旋转。

Blas Level 2

BLAS 调用 映射的库表达式 数学描述 注释
_t_mv x = prod (A, x) 或
x = prod (trans (A), x) 或
x = prod (herm (A), x)		 x <- A x 或
x <- AT x 或
x <- AH x		 计算矩阵与向量的乘积。
_t_sv y = solve (A, x, tag) 或
inplace_solve (A, x, tag) 或
y = solve (trans (A), x, tag) 或
inplace_solve (trans (A), x, tag) 或
y = solve (herm (A), x, tag)或
inplace_solve (herm (A), x, tag)		 y <- A-1 x 或
x <- A-1 x 或
y <- AT-1 x 或
x <- AT-1 x 或
y <- AH-1 x 或
x <- AH-1 x		 求解三角形式的线性方程组，即 A 是三角形的。
_g_mv
_s_mv
_h_mv		 y = a * prod (A, x) + b * y 或
y = a * prod (trans (A), x) + b * y 或
y = a * prod (herm (A), x) + b * y		 y <- a A x + b y 或
y <- a AT x + b y
y <- a AH x + b y		 添加矩阵与向量的缩放乘积。
_g_r
_g_ru
_g_rc		 A += a * outer_prod (x, y) 或
A += a * outer_prod (x, conj (y))		 A <- a x yT + A 或
A <- a x yH + A		 执行秩 1 更新。
_s_r
_h_r		 A += a * outer_prod (x, x) 或
A += a * outer_prod (x, conj (x))		 A <- a x xT + A 或
A <- a x xH + A		 执行对称或 Hermitian 秩 1 更新。
_s_r2
_h_r2		 A += a * outer_prod (x, y) +
 a * outer_prod (y, x)) 或
A += a * outer_prod (x, conj (y)) +
 conj (a) * outer_prod (y, conj (x)))		 A <- a x yT + a y xT + A 或
A <- a x yH + a- y xH + A		 执行对称或 Hermitian 秩 2 更新。

Blas Level 3

BLAS 调用 映射的库表达式 数学描述 注释
_t_mm B = a * prod (A, B) 或
B = a * prod (trans (A), B) 或
B = a * prod (A, trans (B)) 或
B = a * prod (trans (A), trans (B)) 或
B = a * prod (herm (A), B) 或
B = a * prod (A, herm (B)) 或
B = a * prod (herm (A), trans (B)) 或
B = a * prod (trans (A), herm (B)) 或
B = a * prod (herm (A), herm (B))		 B <- a op (A) op (B)，其中
  op (X) = X 或
  op (X) = XT 或
  op (X) = XH		 计算两个矩阵的缩放乘积。
_t_sm C = solve (A, B, tag) 或
inplace_solve (A, B, tag) 或
C = solve (trans (A), B, tag) 或
inplace_solve (trans (A), B, tag) 或
C = solve (herm (A), B, tag) 或
inplace_solve (herm (A), B, tag)		 C <- A-1 B 或
B <- A-1 B 或
C <- AT-1 B 或
B <- A-1 B 或
C <- AH-1 B 或
B <- AH-1 B		 求解三角形式的线性方程组，即 A 是三角形的。
_g_mm
_s_mm
_h_mm		 C = a * prod (A, B) + b * C 或
C = a * prod (trans (A), B) + b * C 或
C = a * prod (A, trans (B)) + b * C 或
C = a * prod (trans (A), trans (B)) + b * C 或
C = a * prod (herm (A), B) + b * C 或
C = a * prod (A, herm (B)) + b * C 或
C = a * prod (herm (A), trans (B)) + b * C 或
C = a * prod (trans (A), herm (B)) + b * C 或
C = a * prod (herm (A), herm (B)) + b * C		 C <- a op (A) op (B) + b C，其中
  op (X) = X 或
  op (X) = XT 或
  op (X) = XH		 添加两个矩阵的缩放乘积。
_s_rk
_h_rk		 B = a * prod (A, trans (A)) + b * B 或
B = a * prod (trans (A), A) + b * B 或
B = a * prod (A, herm (A)) + b * B 或
B = a * prod (herm (A), A) + b * B		 B <- a A AT + b B 或
B <- a AT A + b B 或
B <- a A AH + b B 或
B <- a AH A + b B		 执行对称或 Hermitian 秩 k 更新。
_s_r2k
_h_r2k		 C = a * prod (A, trans (B)) +
 a * prod (B, trans (A)) + b * C 或
C = a * prod (trans (A), B) +
 a * prod (trans (B), A) + b * C 或
C = a * prod (A, herm (B)) +
 conj (a) * prod (B, herm (A)) + b * C 或
C = a * prod (herm (A), B) +
 conj (a) * prod (herm (B), A) + b * C		 C <- a A BT + a B AT + b C 或
C <- a AT B + a BT A + b C 或
C <- a A BH + a- B AH + b C 或
C <- a AH B + a- BH A + b C		 执行对称或 Hermitian 秩 2 k 更新。

存储布局

uBLAS 支持多种不同的存储布局。完整详细信息可以在 类型概述 中找到。大多数类型，如 vector<double> 和 matrix<double>，默认情况下与 C 数组兼容，但也可以配置为包含 FORTAN 兼容的数据。

兼容性

出于兼容性原因，我们为向量和矩阵提供了类似数组的索引。对于某些类型（Hermitian、稀疏等），对于矩阵来说，这可能会很昂贵，因为需要临时的代理对象。

uBLAS 使用 STL 兼容的分配器来分配其容器所需的存储空间。

基准测试结果

下表包含我们其中一项基准测试的结果。此基准测试比较了原生 C 实现（“C 数组”）和一些基于库的实现。基于库的安全变体假设别名，快速变体不使用临时对象，并且在功能上等同于原生 C 实现。除了通用向量和矩阵类之外，基准测试还利用了特殊的类 c_vector 和 c_matrix，这些类旨在避免通过泛型产生的任何开销。

基准测试程序 bench1 使用 GCC 4.0 编译，并在 Athlon 64 3000+ 上运行。时间按比例缩放以获得合理的精度，通过运行 bench1 100。

首先，我们评论维度分别为 3 和 3 x 3 的双精度向量和矩阵的结果。

注释
inner_prod C 数组 0.61 782 一些抽象开销
c_vector 0.86 554
vector<unbounded_array> 1.02 467
vector + vector C 数组 0.51 1122 抽象开销：因子 2
c_vector fast 1.17 489
vector<unbounded_array> fast 1.32 433
c_vector safe 2.02 283
vector<unbounded_array> safe 6.95 82
outer_prod C 数组 0.59 872 一些抽象开销
c_matrix, c_vector fast 0.88 585
matrix<unbounded_array>, vector<unbounded_array> fast 0.90 572
c_matrix, c_vector safe 1.66 310
matrix<unbounded_array>, vector<unbounded_array> safe 2.95 175
prod (matrix, vector) C 数组 0.64 671 无明显的抽象开销
c_matrix, c_vector fast 0.70 613
matrix<unbounded_array>, vector<unbounded_array> fast 0.79 543
c_matrix, c_vector safe 0.95 452
matrix<unbounded_array>, vector<unbounded_array> safe 2.61 164
matrix + matrix C 数组 0.75 686 无明显的抽象开销
c_matrix fast 0.99 520
matrix<unbounded_array> fast 1.29 399
c_matrix safe 1.7 303
matrix<unbounded_array> safe 3.14 164
prod (matrix, matrix) C 数组 0.94 457 无明显的抽象开销
c_matrix fast 1.17 367
matrix<unbounded_array> fast 1.34 320
c_matrix safe 1.56 275
matrix<unbounded_array> safe 2.06 208

我们注意到小型向量和矩阵的性能损失是双重的：首先是使用类的通用抽象开销，然后是使用通用向量和矩阵类时的一小部分损失。关于别名假设的差异也很显著。

接下来，我们评论维度分别为 100 和 100 x 100 的双精度向量和矩阵的结果。

操作 实现 耗时 [秒] MFLOP/秒 注释
inner_prod C 数组 0.64 889 无明显的抽象开销
c_vector 0.66 862
vector<unbounded_array> 0.66 862
vector + vector C 数组 0.64 894 无明显的抽象开销
c_vector fast 0.66 867
vector<unbounded_array> fast 0.66 867
c_vector safe 1.14 501
vector<unbounded_array> safe 1.23 465
outer_prod C 数组 0.50 1144 无明显的抽象开销
c_matrix, c_vector fast 0.71 806
matrix<unbounded_array>, vector<unbounded_array> fast 0.57 1004
c_matrix, c_vector safe 1.91 300
matrix<unbounded_array>, vector<unbounded_array> safe 0.89 643
prod (matrix, vector) C 数组 0.65 876 无明显的抽象开销
c_matrix, c_vector fast 0.65 876
matrix<unbounded_array>, vector<unbounded_array> fast 0.66 863
c_matrix, c_vector safe 0.66 863
matrix<unbounded_array>, vector<unbounded_array> safe 0.66 863
matrix + matrix C 数组 0.96 596 无明显的抽象开销
c_matrix fast 1.21 473
matrix<unbounded_array> fast 1.00 572
c_matrix safe 2.44 235
matrix<unbounded_array> safe 1.30 440
prod (matrix, matrix) C 数组 0.70 813 无明显的抽象开销
c_matrix fast 0.73 780
matrix<unbounded_array> fast 0.76 749
c_matrix safe 0.75 759
matrix<unbounded_array> safe 0.76 749

对于更大的向量和矩阵，使用类的通用抽象开销似乎有所减少，使用通用向量和矩阵类时的一小部分损失似乎仍然存在。关于别名假设的差异仍然可见。

版权 (©) 2000-2002 Joerg Walter, Mathias Koch
使用、修改和分发受 Boost 软件许可协议 1.0 版约束。(请参阅随附文件 LICENSE_1_0.txt 或副本位于 https://boost.ac.cn/LICENSE_1_0.txt )。