F# 与数学：第一部分 - BLAS 和 LAPACK 入门

• 下载基础知识 - 1.99 MB

引言

对线性代数感兴趣的原因有很多，从基础研究到工程学，再到计量经济学和金融分析。无论出于何种原因，吸引 .NET 开发人员进入这个迷人的数学分支，都有许多可能的路径可供选择。不幸的是，通常必须在旅程开始时做出决定，而此时可能还没有经验来做出适当的选择。虽然许多线性代数库都是用纯 .NET 代码开发的，但我们将研究如何使用著名的 BLAS 和 LINPACK 库。这些库历史悠久，可靠性和性能都有良好的记录。不幸的是，这些库很大，它们的命名法对于初学者来说可能显得晦涩难懂，从而望而却步。我们将尝试在开发 .NET 接口库时，总体上阐明这些库的组织和使用。

尽管我们将在示例代码中使用 F#，但我们的目标是创建一个 BLAS 和 LINPACK 接口，该接口可供 C#、Visual Basic 或任何其他您可能喜欢的 .NET 语言使用。

背景

BLAS，即基本线性代数子程序，于 1979 年首次开发并用 Fortran 编程。今天，“BLAS”这个缩写有时指此库，有时指用其他语言开发的类似库，有时指替代库试图实现的临时“标准”。在本文中，我们将使用 BLAS 的原始含义，指代仍然是一个活跃项目并可在 www.NetLib.org 访问的 Fortran 库。最初，BLAS 经常与 LINPACK 结合使用。其组织方式简单明了；核心功能，如矩阵和向量的乘法，在 BLAS 中找到。更复杂的功能，如计算特征值，在 LINPACK 中实现。

您可能已经注意到，自 1979 年以来，计算机硬件已经得到了改进。LAPACK 的开发旨在利用新硬件，而 LINPACK 没有，尤其是在并行性方面。如今，BLAS/LAPACK 设定了数值线性代数的标准，其他软件包对其进行补充和优化，但并未取代它们。通常，无论是说“BLAS”还是“LAPACK”，人们通常都在谈论这一对。正如其名称所示，BLAS 是基础；LAPACK 建立在 BLAS 的功能之上。

该项目受到 F# Powerpack 源代码（可从 www.codeplex.com 下载）提供的单元测试中一些 F# 代码的启发。不出所料，测试的组织方式不利于一个将不断增长并需要维护的库。不幸的是，一些测试代码似乎也有错误，但我假设这没有引起注意，因为这些代码用于测试原生指针功能，而不是其本身。此处呈现的部分代码，特别是 `NativeUtilities.fs` 中找到的包装器，几乎未作更改，并标有 Microsoft 版权声明，表明它们属于 Microsoft 的 CodePlex 许可策略。

Using the Code

在第一部分中，我们将实现实矩阵和向量的乘法，这意味着我们只需要 BLAS。这将使我们能够专注于库的组织结构以及从 .NET 调用 Fortran 库所需的结构。在第二部分中，我们将对库进行扩展，以包含复数矩阵和向量，以及一些更复杂的计算，包括特征值和特征向量的计算，以及线性方程组的解。库的组织结构将允许我们在将来根据需要切换库。例如，如果我们希望将 BLAS 函数移动到 CUBLAS，这是 NVidia 提供的库，用于在其显卡或其高度并行的浮点处理器设备（如 Tesla）上执行超快的矩阵运算。

当然，我们需要 BLAS 和 LINPACK 库的副本。Netlib 提供可下载的 Windows 可执行文件。如果您愿意，可以自己编译源代码，但这涉及在 Windows 上安装和使用 gnu 编译器套件，不适合胆小者。本项目 zip 文件中包含 LAPACK 3.3.0 的已编译 Windows DLL（32 位）。为了方便（我的方便，不是您的），我没有在此编译中包含最新的 LAPACK 扩展，但到您需要它们时，您无论如何都会成为专家。

我们的库将很大程度上遵循 BLAS 和 LAPACK 的命名约定。大多数 .NET 开发人员宁愿调用函数“DoubleMatrixMultiply”而不是 DGEMM，这是理所当然的。但是，您需要一个用于乘三对角矩阵的程序，以及另一个用于三角矩阵的程序。讽刺的是，如果您需要考虑一千多个函数名，那么好的函数名可能会变成不好的函数名。BLAS 和 LINPACK 有一个经过深思熟虑的约定，但只有在一番努力之后，它才会开始显得自然。我们将使用标准名称来指示哪些 F# 函数最终依赖于哪些 BLAS/LAPACK 函数。在本文第一部分的代码中，我们将只使用 DGEMM 和 DGEMV。“D”表示双精度浮点数。“G”表示通用。换句话说，不假定矩阵具有任何特殊特性，例如三对角矩阵。“MM”表示矩阵乘矩阵，“MV”表示矩阵乘向量，但您已经明白了。

在 Fortran 中，所有函数参数都通过引用传递。这意味着我们将大量使用指针，即使对于整数也是如此。另一个复杂之处在于，与 .NET Framework 和 C 语言不同，Fortran 数组是列主序的。这意味着当二维数组在 RAM 中布局时，第 0 行第 0 列的值后面是第 1 行第 0 列。在 .NET 中，RAM 中的数组按第 0 行第 0 列、第 0 行第 1 列等顺序排列。我们需要解决这个问题。当我们为 DGEMM 和 DGEMV 创建声明时，我们注意到似乎有大量潜在冗余的参数。这些参数旨在为函数提供极大的灵活性，但我们将尝试在这里保持简单。好消息是，一旦您熟悉声明参数的含义，同样的参数也会出现在许多其他函数中，工作就会变得更容易。

// dgemm ****************************************************************************
// LAPACK/BLAS primitive matrix/matrix multiply routine 
[<DllImport(@"blas.dll", EntryPoint="dgemm_", 
		CallingConvention=CallingConvention.Cdecl)>] 
extern void DGEMM(char* transa, char* transb,
                  int* m, int* n, int* k, 
                  double* alpha, double* A, int* lda,double* B, 
			int* ldb, double* beta, 
                  double* C, int* ldc);

GNU 工具链在编译非托管 Windows DLL 时喜欢在函数名末尾加上下划线，我们的声明必须完全匹配。函数名当然区分大小写。在此声明中，我们指定了调用约定，但由于 `Cdecl` 是默认值，这更多是为了自文档化而非强制要求。在 64 位 Windows 中，其他调用约定不再存在，尽管这没有什么用处。

前两个参数，`transa` 和 `transb` 是字符，取值“`t`”或“`n`”，取决于您是否想使用矩阵 A 和 B 的转置。虽然不适用于此处，“`c`”也是一个可能的值，可用于请求复数矩阵的共轭转置。这些值不区分大小写。

“`m`”是矩阵 A 的行数，因此也是结果 C 的行数。“`n`”是矩阵 B 的列数，因此也是 C 的列数。“`k`”是矩阵 A 的列数，因此也必须是 B 的行数。

我们稍后会再回到 alpha 和 beta。“`double* A`”是第一个双精度浮点矩阵，`lda` 是 A 的主维。在当前的工作中，`lda` 始终等于 `m`。类似地，“`double* B`”是第二个矩阵，`ldb` 是 B 的主维，同样在我们的工作中等于 k。

'double* C' 是我们存储结果矩阵的矩阵，'ldc' 是它的主维度。然而，C 也可以提供输入。这就是 alpha 和 beta 的作用。DGEMM 计算公式 (alpha * A * B + beta * C)。换句话说，它将 A 和 B 的乘积乘以标量 alpha。通常，我们希望 alpha 为 1.0。此结果将添加到标量 beta 乘以作为输入提供的 C 的任何内容中……通常，我们希望 beta 为 0.0。结果存储在 C 中。

现在我们已经有了声明，我们可以从 F# 调用它了。我们的库的组织结构会与您习惯的有所不同。会有一个层，由 `LapackAPI.fs` 表示，其中包含所有声明。另一方面，我们的库接口将有两个独立的文件。一个 `fsLapack.fs` 将包含方便 F# 程序员使用的包装函数。它的姊妹文件 `clsLapack.fs` 将包含一个可以在不使用任何 F# 特定类型的情况下调用的接口。这可以直接供 C# 或 Visual Basic 程序员使用。注意不幸的前缀。尽管很多人使用“`cls`”作为“class”的匈牙利前缀，但我在这里用它来表示任何 CLS 语言都可以访问的函数。

在这两层之间，每个 BLAS/LAPACK 函数都会有一个单独的文件。尽管这可能看起来有些多余，但随着添加的函数越来越多，从接口到 API 调用的调用序列就变得越来越难以追踪。为每个 BLAS/LAPACK 函数设置一个单独的文件极大地简化了调试和维护。

以下代码展示了项目的标准模式。函数 `DGEMM` 直接调用 `LapackAPI` 中同名的函数声明。请注意，此函数将列主序的 Fortran 矩阵作为输入。第二个函数，用 _C 修饰，接受 F# C 风格的矩阵，并将其准备好导出到 Fortran。由于所有这些矩阵都同时充当输入和输出，因此在我们将结果用于 .NET 代码之前，封送机制将一切设置正确。

注意 `DGEMM` 中大量的“mutable”声明。我们说过所有内容都通过引用传递给 Fortran，F# 不会让你获取指向类型的指针，除非该类型是可变的。即使你不需要“静音”它。

let DGEMM trans alpha (A: FortranMatrix<double>) 
	(B: FortranMatrix<double>) beta (C: FortranMatrix<double>) = 
  let mutable trans = trans   
  Debug.Assert(trans.ToString().ToLower() = "n") 
  Debug.Assert(A.NumRows = B.NumCols)
  let mutable beta = beta 
  let mutable alpha = alpha 
  let mutable m = A.NumRows 
  let mutable n = B.NumCols 
  let mutable k = A.NumCols 
  let mutable lda = m
  let mutable ldb = k 
  let mutable ldc = m
  LapackAPI.DGEMM(&&trans, &&trans, &&m, &&n, &&k, &&alpha, 
	A.Ptr, &&lda, B.Ptr, &&ldb, &&beta, C.Ptr, &&ldc)

   /// Here we build a version that operates on row-major data
let DGEMM_C alpha (A: CMatrix<double>) (B: CMatrix<double>) 
	beta (C: CMatrix<double>) = 
   Debug.Assert(A.NumCols = B.NumRows);
   DGEMM 'n' alpha B.NativeTranspose A.NativeTranspose beta C.NativeTranspose    

如果您是 F# 新手，当我们将值类型按引用发送时，您可能会惊讶地看到两个而不是一个“&”符号。在 F# 中，一个“&”符号用于 .NET 内部，根据您是说 C# 还是 VB，它表示“`ref`”或“`ByRef`”。第二个“&”符号将 .NET 指针转换为原生指针以供 pInvoke 使用。

我们现在可以查看 `fsLapack.fs`，这个文件提供了 BLAS/LAPACK 代码与 F# 世界其他部分之间的接口。这些函数都带有“`fs`”前缀以示区分。显然，如果您愿意，可以更改它。这些函数的关键特性是它们将 .NET 矩阵固定在 RAM 中，这样垃圾回收器就不会在不合适的时候移动它们。这是将指针传递给非托管代码的绝对要求。类似的固定也适用于 .NET 数组。（矩阵在 F# 中定义；数组在公共语言规范中定义。）

固定函数，位于 `NativeUtilities.fs` 中，只是 F# 类型 `PinnedArray` 和 `PinnedArray2` 的包装器，它们本身只是 `GCHandle` 的包装器。如果您想掌握指针和封送处理，您将熟悉 `GCHandle`。`NativeUtitilies.fs` 包含的代码几乎与 Microsoft 的 CodePlex 测试和示例代码未作更改，并相应地拥有版权。

let fsDGEMM (A: matrix) (B: matrix) = 
     let C = Matrix.zero A.NumRows B.NumCols
     let (pA,pB,pC) as ptrs = pinMMM A B C
     try dgemm.DGEMM_C 1.0 pA.NativeArray pB.NativeArray 0.0 pC.NativeArray
     finally
         freeMMM ptrs
     C 

显然，`pinMMM` 是方便的包装器，它可以在单个调用中固定三个矩阵，而 `freeMMM` 则在我们完成后释放这些矩阵的三个原生指针。

同样的函数可以提供数组作为参数，而不是矩阵，这样库就可以直接从 C# 或 Visual Basic 使用，而无需额外的工作。

let clsDGEMM (A:double[,]) (B:double[,]) = 
    let C = Array2D.zeroCreate (Array2D.length1 A) (Array2D.length2 B)
    let (pA,pB,pC) as ptrs = pinA2A2A2 A B C
    try dgemm.DGEMM_C 1.0 pA.NativeArray pB.NativeArray 0.0 pC.NativeArray
    finally
        freeMMM ptrs
    C

请注意，为数组类型的对象提供了第二组固定包装器。

一个警告，或许是预兆

在 Visual Studio 设计环境中运行代码时，Visual Studio 会警告堆栈不平衡。这样的警告应该非常重视。如果您跳过警告，代码会正常运行（至少表面上如此）。如果您在 Visual Studio 外部运行，在 fsi 中运行的脚本代码和编译的代码都会在没有警告或错误的情况下执行。我一直无法确定 Visual Studio 如何或为何看到堆栈不平衡。我曾尝试在循环中使用这段代码乘以一千万对矩阵，以期显现出一些堆栈问题，但未能发现任何问题。

基础源代码项目文件

与本文相关的 zip 文件包含 Visual Studio 2008 和 Visual Studio 2010 的 Visual Studio 解决方案文件。还包含一些批处理文件，它们将启动 F# 解释器环境 fsi 并加载 F# 测试脚本。此外，还包含一些 PDF 文件，描述了您可能需要注意的一些细节。

历史

• 2011年3月24日：首次发布