C 和 C++ 的矩阵库

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

Andrew L. Blais（[email protected]）
研究员和作家

本文将介绍一些目前可在 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
...
void random_fill( CWMatrix &M )
{
int SIZE = M.GetRows() ;
for ( int r=0; r

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

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

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

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

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

清单 5. Blitz 中的随机双精度数填充函数
template
void fill( Array &a, int size, Uniform &u )
{
for ( int i=0; i
for ( int j=0; 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

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

参考资料

请在 David Stewart 的站点下载 Meschach，获得关于 Meschach 的信息，或把电子邮件发送到 [email protected]，在信的正文部分写“send all from c/meschach”（别加引号）。该站点有怎样购买 Meschach 手册“Meschach: Matrix Computations in C”的信息。随机数生成器是基于 Knuth 的滞后的基于 Fibonacci 的生成器。请参阅“Seminumerical Algorithms: The Art of Computer Programming”章节 3.2 到 3.3。

请参阅 The C++ Scalar, Vector, Matrix and Tensor Class Library Standard Page，了解定义矩阵的提议。

请下载 Cooperware Matrix 版本 0.4.2。Harry Kuiper 的网页包含更多关于 Cooperware Matrix 的信息。

请阅读Matrix and Quaternions FAQ，了解什么是四元数。

如果您想了解更多关于四元数和四维旋转的信息，请参阅Quaternion Powers。在为 CwMtx 添加 random_fill 的时候，我用的是 drand48；如果您想了解这方面的更多信息，请参阅手册页 drand48.3。

请下载 Blitz，如果您想了解更多信息，请参阅 Blitz 主页。

请看一个用 Blitz 实现矩阵的示例。

如果您想了解关于邮件列表的更多信息，请参阅Blitz 邮件列表。

如果您想了解更多关于 Blitz 中实现的随机数生成器的信息，请阅读“Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator”，ACM Transactions on Modeling and Computer Simulation，第 8 卷，第 1 号，1998 年 1 月，第 3 到 30 页。

请看一看Matrix TCL Lite（免费）或 Matrix TCL Pro（收费）。前者能在我喜欢的编译器 g++ 2.96 上运行，但由于库本身没有优化，所以我的分析不适用它。后者是优化的，但不能在我的编译器上运行，所以我的分析也不适用它。

获得面向对象库的列表，用于科学计算。

请阅读 Andrew 的文章“An introduction to neural networks”（developerWorks，2001 年 7 月）。

Improving the memory-system performance of sparse-matrix vector multiplication 讲述在超标量 RISC 处理器上用大于二的因子来提高矩阵相乘的速度的方法。

如果您想了解更多关于提高稀疏矩阵算法的信息，请阅读Fast and effective algorithms for graph partitioning and sparse-matrix ordering。

请下载WebSphere Studio Application Developer Linux 版的基本版的试用版。

请在 developerWorks Linux 专区查阅更多 Linux 的文章和教程。

作者简介
Andrew L. Blais 把他的时间分配在四件事上：对他的儿子 A. Van Blais 进行家教；当 K.J. Krawczyk-Blais 的丈夫；为 Gnosis Software，Inc. 研究和写作；在 Anna Maria College 教哲学和宗教。有时候，他会休息一下。您可以通过 [email protected] 与 Andrew 联系。