C 和 C++ 的开源矩阵库Meschach、Cooperware 和 Blitz

评估和比较 Meschach、Cooperware 矩阵和 Blitz

本文将介绍一些目前可在 Linux 环境中使用的开放源代码 C/C++ 矩阵库。在这里具
体讨论的库有三个：Meschach 库为 C 代码编写的项目提供例程，用于矩阵和向量的
运算；Cooperware Matrix（CwMtx）库可用于 C++ 代码编写；Blitz 库为 C++ 提供
可用整数、浮点数、复数和规范的用户定义的类型的 n 维数组类。Andrew Blais 是
Gnosis，Inc. 的研究员和作家，他在 developerWorks 已投稿多次，从事神经网络
方面的工作。
本文假设读者对 C/C++ 有一定的了解并非常关注 C/C++ 本身没有矩阵功能。您可能
在分析计量经济学的数据或模拟雨林。对于我来说，我正从事神经网络的研究而一两
个矩阵可以大大简化神经网络的实现。虽然 C/C++ 包括可被看作矩阵（例如数组和
标准库中的向量、列表和图）的容器，但是，真正是矩阵的容器将使手中的任务变得
很容易。所以我们将介绍三个开放源代码库供您选择，它们不要求您从头开始构建矩
阵，但的确让您使用矩阵库。如果您期望用从未见过的方法来使用它，那就太好了。

Meschach：用 C 的选择
对于用 C 编写代码的项目，Meschach（读作 me-shark）提供了例程，用于矩阵和向
量的运算。它的优点是能在 Linux 和多数其它操作系统下编译，并且在版权保护下
可公开获得，只要您作出例行的承认并报告错误。Meschach 可以解稠密或稀疏线性
方程组、计算特征值和特征向量和解最小平方问题，另外还有其它功能。它为双精度
数和复数提供了近 400 个函数。它提供的教程以说明性的小案例研究的形式介绍了
这些函数。David Stewart 和 Zbigniew Leyk 通过一些主题的讨论来介绍 Meschach
，这些主题包括超定方程组的广义的最小平方方程解答器（generalized least square
equation solver for over-determined equations）和涉及稀疏矩阵的问题。他们
的教程还包括三维矩阵和错误报告等稍稍高级的主题。

对象和类函数往往与代码关联，C 结构可能看起来有点神秘，所以 C 库往往不被作
为可行的解决方案。但是作为对此的反击，这个库的组织很合理，所以不应不加思索
地就舍弃它。在下载 Meschach 后过了一刻钟，我就可以制造、填充和显示矩阵了（
在概念上等同于创建 Hello World! 程序）。您可以参考一本名为“Meschach: Matrix
Computations in C”的便宜的印刷品手册（请参阅本文末尾的参考资料）。特别是
测试程序“torture”，其中包括不少有益的线索。

矩阵可被容易地发送到文件或标准输出。Meschach 能计算快速傅立叶变换（Fast Fourier
Transform）、提取列和行以及计算对称矩阵的特征值。您可以在矩阵中填充随机整
数和复数。信不信由您，该库甚至还有矩阵相加的工具。Meschach 有一个返回在 [
0,1) 之间的随机双精度数的函数，它是 Meschach 的一个特色，可用于简化明显的
踏脚石程序的编写。虽然 Meschach 有一个用 1.0 填充矩阵的函数，但不幸的是，
它没有用一个任意的双精度数填充矩阵的函数，也没有用随机双精度数填充矩阵的函
数。不过，添加它们是容易的。

清单 1. 用一个任意的双精度数或随机双精度数来填充矩阵
MAT *m_fill( MAT *A, double x)
/* MAT *m_random_fill( MAT *A ) */
{
int i, j;
for ( i = 0; i < A->m; i++ ) for ( j = 0; j < A->n; j++ )
{ A->me[j] = x ; }
/* { A->me[j] = m_random ; } */
return A;
}

所以，正如我们所看到的，Meschach 代码可被容易地扩展，尽管在理想的世界中应
该更加详细地注释代码。但是如果您正好在用 C 做计算工作并且要用到矩阵，那么
这是很有用的库。

Cooperware 矩阵：基础知识
如果您要用 C++ 编写面向对象的代码并且您认为概念的清晰比速度更重要，那么 Harry
Kuiper 的 Cooperware Matrix（CwMtx）能很好地运行。在这里讨论的三个矩阵库
中，我觉得它的概念性体系结构最容易理解。在构造矩阵时，您直接使用：

清单 2. 构造矩阵
CWMatrix A( rows, columns ) ;

在这里考虑的三个库中，以三个评估性任务为标准，CwMtx 的性能最差，我们将在速
度部分详细地讨论。但是如果清晰比性能更重要（例如当您要确信您的数据被正确地
处理），那么 CwMtx 是很好的选择。先使它正确，再使它快速。

CwMtx 中的矩阵包括向量和方阵，其中向量包括空间向量和四元数。一个矩阵可被映
射到另一个矩阵、用某个元素来填充、转置和进行常见的数学运算。Kuiper 原先用
CwMtx 来模拟用离散的交互式状态机器构建的系统。除了必须的矩阵类，还有四元
数类。对于明显问题的回答是，仅当 q = r + xi + yi + zi 时 q 才是四元数，其
中 r 是实数，i 是 -1 的平方根，x、y 和 z 是复数。四元数可能把三维旋转的概
念推广到四维（请参阅参考资料，其中有四元数的参考资料的链接）。

CwMtx 没有内置的随机数生成器，也没有用随机元素填充矩阵的类函数。但是，它是
免费的且它的发布受 GNU LGPL 许可证的保护，所以如果您愿意的话，您有创建这些
（生成器和类函数）的自由。如果仅仅是用随机元素填充矩阵，那么以下代码是不错
的和容易的选择。

清单 3. 用随机元素填充矩阵
#include <stdlib.h>
...
void random_fill( CWMatrix &M )
{
int SIZE = M.GetRows() ;
for ( int r=0; r<SIZE; r++ ) for ( int c=0; c<SIZE; c++ ) { M[r][c]= drand48
(); }
}

这里的文档并不详细，但文档是清楚的并且组织得也很好。您能容易地找到类层次结
构、构造函数、成员函数选项等信息。虽然没有教程，但您不会想念它；由于有了文
档和测试程序，您并不需要它。

Blitz：和 Fortran 一样快？
Blitz 是另一个 C++ 库，它的发行受 GNU GPL 的保护，您可以用它来免费地创建对
象。它支持 KAI、Intel、gcc、Metroworks 和 Cray 3.0.0.0 C++ 编译器，它还提
供 n 维数组类，这些类可包括整数、浮点数、复数和用户定义的表现良好的类型。
它的构造函数比 CwMtx 的构造函数更复杂，下面的示例证明了这一点：

清单 4. Blitz 构造函数
Array<double,2> A(4,7) ;

示例中创建了包含双精度数的秩为 2 的 4x7 数组。但是这样有点不清楚，所以 Blitz
使您把数组看作矩阵。还有，它没有实现许多矩阵函数。例如，它没有返回矩阵的
特征值的函数。它也没有用随机双精度数来填充数组或矩阵的函数。但是，Blitz 确
实有两个基本优点。

一个优点是它的广度。通过使用它自身的功能可以容易地实现和构造随机双精度数填
充函数，如下面的示例所示：

清单 5. Blitz 中的随机双精度数填充函数

template <class Array, class Uniform>
void fill( Array &a, int size, Uniform &u )
{
for ( int i=0; i<size; i++ )
for ( int j=0; j<size; j++ )
{ a(i,j) = u.random() ; }
}

Blitz 的 Uniform 类提供返回在 [0,1) 间的双精度数的成员函数。它还提供三种访
问数组元素的方法：标准索引、创建子数组和切片（slicing），切片能产生维数更
小的数组片段的视图。Blitz 还处理标准计算器类型函数，所以数组可在标准输出上
显示，而且可从文件中读入或发送到文件。Laplacian、坡度（gradient）和 Jacobian
运算符只是 Blitz 的模版（stencil）函数的三个示例。

Blitz 的另一个优点是它的速度。根据所用的编译器，C++ 的性能可以赶上 Fortran
，而 Fortran 在科学和工程计算方面所表现的高性能是出名的。看一看下表中的比
较，但是请阅读后面的速度部分，其中分析了这些数据和基于这些数据的性能表现。

表 1. 在不同平台上的 Blitz 性能
平台        编译器     高速缓存外 高速缓存内
HPC-160     KAI C++    100.2%     97.5%
Pentium II  egcs       98.4%      79.6%
Cray T3E    KAI C++    95.7%      98.1%
Origin 2000 KAI C++    88.1%      79.8%

Blitz 带有一本手册，格式是 HTML 和 Postscript，但不幸的是没有教程。然而，
有不少说明性的代码，从中可以了解 Blitz 语法的细微差别。类的参考资料提供通
常的信息。还有几个有用的邮件列表，已被归档，可供搜索（请参阅参考资料）。

速度
库的评价标准有库的功能性资源、文档、库的教程质量、库的扩展难度等。库的评价
标准还有性能和／或速度。但是，有时候比较它们是困难的，因为（我们的示例就是
这种情况）它们并不是都用相同的语言来编写的，而且它们相同的自身功能也不多。
在我们这里，三个库有足够的重复部分，这使我们能用三个比较简单但能说明问题的
任务来比较它们的速度，我们将在下面三个示例中显示和讨论这些任务。

清单 6. 第 1 个任务
For ( d=2; d<7; d++)
Construct 3 dxd matrices: A, B, C
Start Clock: Do the following one million times:
Fill A and B with 1.0s
Let C = A + B
Stop clock: Report elapsed time in seconds.

使用我们的库来实现并执行这个算法产生以下结果：

表 2. 第 1 个任务的结果
库        2x2   3x3   4x5   5x5   6x6
Blitz     0.40  0.48  0.62  0.75  0.91
CwMtx     2.64  3.57  4.58  5.60  6.60
Meschach  0.17  0.27  0.45  0.60  0.79

尽管 CwMtx 的体系结构易于理解，但不幸的是，它在这里的表现并不好。虽然 Blitz
的表现不如 Meschach，但是值得注意的是 Blitz 的性能大大超过了它的面向对象
的对手 CwMtx。

如前所述，Meschach 和 Blitz 有提供随机双精度数的函数（随机数生成器），而
CwMtx 自身没有产生随机数的功能。考虑到随机化在某些基于矩阵的模拟中的关键作
用，研究这些库在调用随机化的情况下的表现是有益的。

清单 7. 第 2 个任务
For ( d=2; d<7; d++)
Construct 3 dxd matrices: A, B, C
Start Clock: Do the following one million times:
Fill A and B with random doubles (using library RNG, if any)
Let C = A + B
Stop clock: Report elapsed time in seconds.

我们的库的表现如下：

表 3. 第 2 个任务的结果
库        2x2   3x3   4x5   5x5   6x6
Blitz     0.87  1.71  2.83  4.34  6.13
CwMtx     3.67  5.87  8.59  11.93 15.61
Meschach  0.42  0.80  1.32  1.86  2.52

Blitz 的表现再次比 Meschach 差，但它的表现超过了它的面向对象的对手 CwMtx，
而且这种差距令人吃惊。我们来看看第三个评价任务，以免您认为这是因为随机数生
成器的性能有所不同。

清单 8. 第 3 个任务
For ( d=2; d<7; d++)
Construct 3 dxd matrices: A, B, C
Start Clock: Do the following one million times:
Fill A and B with random doubles (using shared RNG)
Let C = A + B
Stop clock: Report elapsed time in seconds.

我们的库的性能排行榜与前面两个任务的相同：

表 4. 第 3 个任务的结果
库        2x2   3x3   4x5   5x5   6x6
Blitz     1.31  2.62  4.50  6.85  9.71
CwMtx     3.67  5.87  8.59  11.93 15.61
Meschach  1.17  2.45  4.28  6.63  9.40

如您所预料的那样，CwMtx 的性能排名没有改变。而且，Blitz 和 Meschach 的排名
先后顺序也没有变。如果原始速度是决定性因素，那么这些库的排名现在已经很清楚
。

在 Red Hat Linux 7.1 上安装和编译
为了您的方便，以下列出的是安装和编译这三个库的注解。下载的链接可在参考资料
中找到。

清单 9. Blitz
tar zxf blitz-0.5beta3.tar.gz
cd blitz-0.5beta3
./configure --with-cxx=gcc
make all
cp -a blitz/ /usr/local/include/ ( or whatever you wish )
cp -a lib/ /usr/local/include/blitz/
Compile with: g++ -O2 pgm.cpp -o pgm

清单 10. Meschach
unzip -q mesch12b.zip -d mesch12b
cd mesch12b
./configure
make basic
mkdir /usr/local/include/meschach ( or whatever you wish )
cp *.h meschach.a /usr/local/include/meschach/
Compile with: gcc -O2 pgm.c /usr/local/include/meschach/meschach.a -o pgm

清单 11. CwMtx
tar zxf cwmtx-0.3.0.tar.gz
cd cwmtx-0.3.0
Open Makefile, and set INSTALL_DIR to /usr/local/include/cwmtx ( or whatever
you wish )
Execute: make install
Compile with: g++ -O2 pgm.cpp -o pgm

总结
我们已经了解了三个矩阵库的特色和它们本身的功能的详细信息。我们也看到了它们
的一些功能缺点以及克服这些缺点的方法。我还为您提供了一些简单测试，这些测试
所提供的原始的量化数据可能在您代码编写选择时对您有所帮助，但最终的决定要由
您来作出，而且应该依据您的项目的自身特点以及在给出的环境中库的速度来作出决
定。