现代C++25-数字计算：介绍线性代数和数值计算库

科学计算在今天已经完全可以使用 C++ 了。我不是从事科学计算这一领域的工作的，不过，在工作中也多多少少接触到了一些计算相关的库。今天，我就给你介绍几个有用的计算库。

Armadillo

说到计算，你可能首先会想到矩阵、矢量这些东西吧？这些计算，确实就是科学计算中的常见内容了。这些领域的标准，即是一些 Fortran 库定下的，如：

BLAS [1]

LAPACK [2]

ARPACK [3]

它们的实现倒不一定用 Fortran，尤其是 BLAS：

OpenBLAS [4] 是用汇编和 C 语言写的

Intel MKL [5] 有针对 Intel 的特定 CPU 指令集进行优化的汇编代码

Mir GLAS [6] 是用 D 语言写的

不管实现的方法是哪一种，暴露出来的函数名字是这个样子的：

ddot

dgemv

dsyrk

sgemm

……

这个接口的唯一好处，应该就是，它是跨语言并且跨实现的😅。所以，使用这些函数时，你可以切换不同的实现，而不需要更改代码。唯一需要修改的，通常就是链接库的名字或位置而已。

假设我们需要做一个简单的矩阵运算，对一个矢量进行旋转：

PRP′​=[1 0​]=[cos(θ)sin(θ)​−sin(θ)cos(θ)​]=R⋅P​

这么一个简单的操作，用纯 C 接口的 BLAS 来表达，有点痛苦：你需要使用的大概是 dgemv_ 函数，而这个函数需要 11 个参数！我查阅了一下资料之后，也就放弃了给你展示一下如何调用 dgemv_ 的企图，我们还是老老实实地看一下在现代 C++ 里的写法吧：

#include <armadillo>

#include <cmath>

#include <iostream>

using namespace std;

int main()

{

// 代表位置的向量

arma::vec pos{1.0, 0.0};

// 旋转矩阵

auto& pi = arma::datum::pi;

double angle = pi / 2;

arma::mat rot = {

{cos(angle), -sin(angle)},

{sin(angle), cos(angle)}};

cout << "Current position:\n"

<< pos;

cout << "Rotating "

<< angle * 180 / pi

<< " deg\n";

arma::vec new_pos = rot * pos;

cout << "New position:\n"

<< new_pos;

}

这就是使用 Armadillo [7] 库来实现矢量旋转的代码。这个代码，基本就是上面的数学公式的一一对应了。代码相当直白，我只需要稍稍说明一下：

所有的 Armadillo 的类型和函数都定义在 arma 名空间下。

Armadillo 在 arma::datum 下定义了包括 pi 和 e 在内的一些数学常量。

vec 是矢量类型，mat 是矩阵类型，这两个类型实际上是 Col<double> 和 Mat<double> 的缩写别名。

Armadillo 支持使用 C++11 的列表初始化语法来初始化对象。

Armadillo 支持使用流来输出对象。

上面代码的输出为：

Current position:

1.0000

0

Rotating 90 deg

New position:

6.1232e-17

1.0000e+00

输出里面的 6.1232e-17 是浮点数表示不精确的后果，把它理解成 0 就对了。

我们上面已经提到了 vec 实际上是 Col<