网络浏览器中的 LAPACK

这篇文章最初发布在 Quansight Labs 博客上，经 Quansight 许可修改并重新发布于此处。

Web 应用程序正在迅速崛起，成为高性能科学计算和支持 AI 的最终用户体验的新领域。支撑机器学习/人工智能革命的是​​线性代数，它是关于线性方程及其在向量空间和矩阵中的表示的数学分支。 LAPACK（“Linear Algebra Package”）是数值线性代数的基础软件库，提供强大的、经过实战检验的常见矩阵运算实现。尽管 LAPACK 是大多数数值计算编程语言和库的基础组件，但针对网络独特限制量身定制的全面、高质量的 LAPACK 实现尚未实现。那是……到现在为止。

今年早些时候，我有幸成为 Quansight Labs 的暑期实习生，该实验室是 Quansight 的公益部门，也是科学 Python 生态系统的领导者。在实习期间，我致力于为 stdlib 添加初始 LAPACK 支持，stdlib 是一个用 C 和 JavaScript 编写的科学计算基础库，并针对在 Web 浏览器和其他 Web 原生环境（例如 Node.js 和 Deno）中使用进行了优化。在这篇博文中，我将讨论我的旅程、一些预期和意外（！）的挑战以及未来的道路。我希望，如果运气好的话，这项工作能够为使网络浏览器成为数值计算和机器学习的一流环境提供关键的构建模块，并预示着更强大的人工智能网络应用程序的未来。

听起来很有趣吗？我们走吧！

什么是 stdlib？

熟悉 LAPACK 的本博客读者可能不太熟悉 Web 技术的狂野世界。对于那些来自数值和科学计算领域并熟悉科学 Python 生态系统的人来说，最简单的方式是将 stdlib 视为 NumPy 和 SciPy 模型中的开源科学计算库。它提供多维数组数据结构以及数学、统计和线性代数的相关例程，但使用 JavaScript 而不是 Python 作为其主要脚本语言。因此，stdlib 专注于 Web 生态系统及其应用程序开发范例。这种关注需要一些有趣的设计和项目架构决策，这使得 stdlib 与为数值计算设计的更传统的库相比相当独特。

以 NumPy 为例，NumPy 是一个单一的整体库，其中所有组件（除了可选的第三方依赖项（例如 OpenBLAS）之外）形成一个不可分割的单元。如果不安装所有 NumPy，就无法简单地安装 NumPy 例程来进行数组操作。如果您正在部署一个仅需要 NumPy 的 ndarray 对象及其一些操作例程的应用程序，那么安装和捆绑所有 NumPy 意味着包含大量“死代码”。用 Web 开发的术语来说，我们会说 NumPy 不是“tree shakeable”。对于正常的 NumPy 安装，这意味着至少需要 30MB 的磁盘空间，对于排除所有调试语句的自定义构建至少需要 15MB 的磁盘空间。对于 SciPy，这些数字可以分别增加到 130MB 和 50MB。不用说，在 Web 应用程序中仅提供几个功能的 15MB 库是不可能的，特别是对于需要将 Web 应用程序部署到网络连接较差或内存受限的设备的开发人员而言。

考虑到 Web 应用程序开发的独特限制，stdlib 采用自下而上的设计方法，其中每个功能单元都可以独立于代码库中不相关和未使用的部分进行安装和使用。通过采用可分解的软件架构和彻底的模块化，stdlib 使用户能够安装和使用他们所需要的内容，除了所需的 API 集及其显式依赖项之外几乎没有多余的代码，从而确保更小的内存占用、捆绑尺寸和更快的部署。

举个例子，假设您正在使用两个矩阵堆栈（即三维立方体的二维切片），并且您想要选择每隔一个切片并执行常见的 BLAS 运算 y = a * x，其中 x 和 y 是 ndarray，a 是标量常量。要使用 NumPy 执行此操作，您首先需要安装所有 NumPy

pip install numpy

然后执行各种操作

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

使用 stdlib，除了能够将项目安装为整体库之外，您还可以将各种功能单元安装为单独的包

npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy

然后执行各种操作

// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';

// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;

// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');

// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);

重要的是，您不仅可以独立安装 stdlib 超过 4,000 个软件包中的任何一个，还可以通过分叉关联的 GitHub 存储库来修复、改进和重新混合其中任何一个软件包（例如，请参阅 @stdlib/ndarray-fancy） ）。通过定义显式的抽象层和依赖树，stdlib 使您可以自由地为您的应用程序选择正确的抽象层。在某些方面，这是一个简单的想法，如果您习惯了传统的科学软件库设计，也许是非正统的想法，但是，当与 Web 平台紧密集成时，它会产生强大的后果并创造令人兴奋的新可能性！

WebAssembly 怎么样？

好吧，也许你的兴趣已经激起； stdlib 似乎很有趣。但这与网络浏览器中的 LAPACK 有什么关系呢？嗯，去年夏天我们的目标之一是应用 stdlib 精神——小型、范围狭窄的软件包，只做一件事，并做好一件事——将 LAPACK 引入网络。

但是等等，你说！这是一项极端的任务。 LAPACK 规模庞大，约有 1,700 个例程，在合理的时间范围内实施其中的 10% 都是一项重大挑战。将 LAPACK 编译为 WebAssembly 不是更好吗？WebAssembly 是 C、Go 和 Rust 等编程语言的可移植编译目标，可以在网络上部署，然后就到此为止了？

不幸的是，这种方法存在几个问题。

1. 将 Fortran 编译为 WebAssembly 仍然是一个正在积极开发的领域（请参阅 1、2、3、4 和 5）。在撰写本文时，常见的方法是使用 f2c 将 Fortran 编译为 C，然后执行单独的编译步骤将 C 转换为 WebAssembly。然而，这种方法是有问题的，因为 f2c 仅完全支持 Fortran 77，并且生成的代码需要大量修补。开发基于 LLVM 的 Fortran 编译器的工作正在进行中，但仍然存在差距和复杂的工具链。
2. 正如上面关于 Web 应用程序中的整体库的讨论中提到的，LAPACK 的庞大是问题的一部分。即使编译问题得到解决，在仅需要使用一两个 LAPACK 例程的 Web 应用程序中包含包含所有 LAPACK 的单个 WebAssembly 二进制文件也意味着大量死代码，导致加载时间变慢并增加内存消耗。
3. 虽然人们可以尝试将单个 LAPACK 例程编译为独立的 WebAssembly 二进制文件，但这样做可能会导致二进制膨胀，因为多个独立二进制文件可能包含来自公共依赖项的重复代码。为了减轻二进制膨胀，可以尝试执行模块拆分。在这种情况下，首先将公共依赖项分解为包含共享代码的独立二进制文件，然后为各个 API 生成单独的二进制文件。虽然在某些情况下适用，但这种方法很快就会变得笨拙，因为这种方法需要在加载时通​​过将一个或多个模块的导出与一个或多个其他模块的导入缝合在一起来链接各个 WebAssembly 模块。这不仅很乏味，而且这种方法还会带来性能损失，因为当 WebAssembly 例程调用导入的导出时，它们现在必须跨入 JavaScript，而不是保留在 WebAssembly 中。听起来很复杂？是的！
4. 除了专门对标量输入参数进行操作的 WebAssembly 模块（例如，计算单个数字的正弦值）之外，每个 WebAssembly 模块实例都必须与 WebAssembly 内存关联，该内存以 64KiB 的固定增量分配（即“页面”） ”）。重要的是，截至这篇博文，WebAssembly 内存只能增长而永远不会收缩。由于当前没有向主机释放内存的机制，因此 WebAssembly 应用程序的内存占用只会增加。这两个方面的结合增加了分配从未使用过的内存的可能性以及内存泄漏的普遍性。
5. 最后，虽然 WebAssembly 功能强大，但它需要更陡峭的学习曲线和一组更复杂且经常快速发展的工具链。在最终用户应用程序中，JavaScript（一种 Web 原生动态编译编程语言）和 WebAssembly 之间的接口进一步增加了复杂性，尤其是在必须执行手动内存管理时。

为了帮助说明最后一点，让我们回到 BLAS 例程 daxpy，它执行运算 y = a*x y，其中 x 和 y 是跨步向量，a 是标量常量。如果用 C 实现，基本实现可能类似于以下代码片段。

pip install numpy

编译到 WebAssembly 并将 WebAssembly 二进制文件加载到我们的 Web 应用程序中后，我们需要执行一系列步骤，然后才能从 JavaScript 调用 c_daxpy 例程。首先，我们需要实例化一个新的 WebAssembly 模块。

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

接下来，我们需要定义模块内存并创建一个新的 WebAssembly 模块实例。

npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy

创建模块实例后，我们现在可以调用导出的 BLAS 例程。但是，如果数据是在模块内存之外定义的，我们首先需要将该数据复制到内存实例，并且始终以小端字节顺序执行此操作。

// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';

// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;

// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');

// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);

现在数据已写入模块内存，我们可以调用 c_daxpy 例程。

void c_daxpy(const int N, const double alpha, const double *X, const int strideX, double *Y, const int strideY) {
    int ix;
    int iy;
    int i;
    if (N <= 0) {
        return;
    }
    if (alpha == 0.0) {
        return;
    }
    if (strideX < 0) {
        ix = (1-N) * strideX;
    } else {
        ix = 0;
    }
    if (strideY < 0) {
        iy = (1-N) * strideY;
    } else {
        iy = 0;
    }
    for (i = 0; i < N; i++) {
        Y[iy] += alpha * X[ix];
        ix += strideX;
        iy += strideY;
    }
    return;
}

最后，如果我们需要将结果传递到不支持 WebAssembly 内存“指针”（即字节偏移）的下游库（例如 D3）以进行可视化或进一步分析，我们需要从模块复制数据内存回原始输出数组。

const binary = new UintArray([...]);

const mod = new WebAssembly.Module(binary);

仅仅计算 y = a*x y 就需要大量工作。相反，与纯 JavaScript 实现相比，它可能类似于以下代码片段。

// Initialize 10 pages of memory and allow growth to 100 pages:
const mem = new WebAssembly.Memory({
    'initial': 10,  // 640KiB, where each page is 64KiB
    'maximum': 100  // 6.4MiB
});

// Create a new module instance:
const instance = new WebAssembly.Instance(mod, {
    'env': {
        'memory': mem
    }
});

通过 JavaScript 实现，我们可以使用外部定义的数据直接调用 daxpy，而无需上面 WebAssembly 示例中所需的数据移动。

// External data:
const xdata = new Float64Array([...]);
const ydata = new Float64Array([...]);

// Specify a vector length:
const N = 5;

// Specify vector strides (in units of elements):
const strideX = 2;
const strideY = 4;

// Define pointers (i.e., byte offsets) for storing two vectors:
const xptr = 0;
const yptr = N * 8; // 8 bytes per double

// Create a DataView over module memory:
const view = new DataView(mem.buffer);

// Resolve the first indexed elements in both `xdata` and `ydata`:
let offsetX = 0;
if (strideX < 0) {
    offsetX = (1-N) * strideX;
}
let offsetY = 0;
if (strideY < 0) {
    offsetY = (1-N) * strideY;
}

// Write data to the memory instance:
for (let i = 0; i < N; i++) {
    view.setFloat64(xptr+(i*8), xdata[offsetX+(i*strideX)], true);
    view.setFloat64(yptr+(i*8), ydata[offsetY+(i*strideY)], true);
}

至少在这种情况下，WebAssembly 方法不仅不太符合人体工程学，而且正如所预期的，考虑到所需的数据移动，还会对性能产生负面影响，如下图所示。

图 1：对于增加数组长度（x 轴）的 BLAS 例程 daxpy，stdlib 的 C、JavaScript 和 WebAssembly (Wasm) 实现的性能比较。在Wasm（复制）基准测试中，输入和输出数据被复制到Wasm内存中或从Wasm内存中复制，导致性能较差。

在上图中，我显示了 BLAS 例程 daxpy 的 stdlib 的 C、JavaScript 和 WebAssembly (Wasm) 实现的性能比较，以增加数组长度（沿 x 轴枚举）。 y 轴显示相对于基线 C 实现的标准化速率。在 Wasm 基准测试中，输入和输出数据直接在 WebAssembly 模块内存中分配和操作，并且在 Wasm（复制）基准测试中，输入和输出数据复制到 WebAssembly 模块内存中或从 WebAssembly 模块内存中复制，如上所述。从图表中，我们可以观察到以下几点：

1. 一般来说，得益于高度优化的即时 (JIT) 编译器，JavaScript 代码在精心编写的情况下，执行速度仅比本机代码慢 2 到 3 倍。对于松散类型、动态编译的编程语言来说，这个结果令人印象深刻，并且至少对于 daxpy 来说，在不同的数组长度上保持一致。
2. 随着数据大小以及 WebAssembly 模块中花费的时间量的增加，WebAssembly 可以接近原生（~1.5 倍）速度。此结果更符合预期的 WebAssembly 性能。
3. 虽然 WebAssembly 可以实现接近本机的速度，但数据移动要求可能会对性能产生不利影响，正如在 daxpy 中观察到的那样。在这种情况下，精心设计的 JavaScript 实现可以避免此类要求，即使不是更好，也可以实现相同的性能，就像 daxpy 的情况一样。

总体而言，WebAssembly 可以提供性能改进；然而，该技术并不是灵丹妙药，需要谨慎使用才能实现预期收益。即使提供卓越的性能，这种收益也必须与复杂性增加、潜在的更大的捆绑包大小和更复杂的工具链的成本相平衡。对于许多应用程序来说，简单的 JavaScript 实现就可以了。

彻底的模块化

现在我已经起诉了反对将整个 LAPACK 编译为 WebAssembly 的案件并就此结束，那我们该怎么办？好吧，如果我们要拥抱 stdlib 精神，那么我们就需要彻底的模块化。

拥抱彻底的模块化意味着要认识到最好的东西是高度上下文相关的，并且根据用户应用程序的需求和限制，开发人员需要灵活地选择正确的抽象。如果开发人员正在编写 Node.js 应用程序，这可能意味着绑定到硬件优化库，例如 OpenBLAS、Intel MKL 或 Apple Accelerate，以获得卓越的性能。如果开发人员正在部署需要一小组数字例程的 Web 应用程序，那么 JavaScript 可能是适合该工作的工具。如果开发人员正在开发大型的资源密集型 WebAssembly 应用程序（例如，用于图像编辑或游戏引擎），那么能够轻松编译单个例程作为大型应用程序的一部分将至关重要。简而言之，我们想要一个完全模块化的 LAPACK。

我的任务是为这样的努力奠定基础，解决问题并找到差距，并希望让我们更接近网络上的高性能线性代数。但是彻底的模块化是什么样的呢？这一切都始于基本的功能单元，包​​。

stdlib 中的每个包都是独立的，包含共同本地化的测试、基准测试、示例、文档、构建文件和关联的元数据（包括任何依赖项的枚举），并与外界定义清晰的 API 界面。为了向 stdlib 添加 LAPACK 支持，这意味着为每个 LAPACK 例程创建一个单独的独立包，其结构如下：

pip install numpy

简单来说，

• benchmark：包含微基准的文件夹，用于评估相对于参考实现（即参考 LAPACK）的性能。
• docs：包含辅助文档的文件夹，其中包括 REPL 帮助文本和定义类型化 API 签名的 TypeScript 声明。
• examples：包含可执行演示代码的文件夹，除了充当文档之外，还可以帮助开发人员健全性检查实现行为。
• include：包含 C 头文件的文件夹。
• lib：包含 JavaScript 源实现的文件夹，其中 index.js 作为包入口点，其他 *.js 文件定义内部实现模块。
• src：包含 C 和 Fortran 源实现的文件夹。每个模块化 LAPACK 包都应包含稍作修改的 Fortran 参考实现（F77 到自由格式的 Fortran）。 C 文件包括遵循 Fortran 参考实现的纯 C 实现、用于调用 Fortran 参考实现的包装器、用于在服务器端应用程序中调用硬件优化库（例如 OpenBLAS）的包装器，以及用于调用已编译的本机绑定来自 Node.js 中的 JavaScript 或兼容的服务器端 JavaScript 运行时的 C。
• test：包含单元测试的文件夹，用于测试 JavaScript 和本机实现中的预期行为。本机实现的测试是用 JavaScript 编写的，并利用本机绑定实现 JavaScript 和 C/Fortran 之间的互操作。
• binding.gyp/include.gypi：用于编译 Node.js 原生插件的构建文件，它提供了 JavaScript 和原生代码之间的桥梁。
• manifest.json：stdlib内部C和Fortran编译源文件包管理的配置文件。
• package.json：包含包元数据的文件，包括外部包依赖项的枚举以及用于基于浏览器的 Web 应用程序的纯 JavaScript 实现的路径。
• README.md：包含包的主要文档的文件，其中包括 API 签名以及 JavaScript 和 C 接口的示例。

考虑到 stdlib 苛刻的文档和测试要求，为每个例程添加支持是一项相当大的工作量，但最终结果是健壮的、高质量的，最重要的是，适合作为科学计算基础的模块化代码在现代网络上。但预赛就够了！让我们言归正传吧！

多阶段方法

基于之前为 stdlib 添加 BLAS 支持的努力，我们决定在添加 LAPACK 支持时遵循类似的多阶段方法，其中我们首先优先考虑 JavaScript 实现及其相关的测试和文档，然后，一旦存在测试和文档、回填 C 和 Fortran 实现以及任何与硬件优化库相关的本机绑定。这种方法使我们能够在董事会上提出一些早期观点，可以这么说，快速将 API 提供给用户，建立强大的测试程序和基准，并在深入研究构建工具链和性能之前调查工具和自动化的潜在途径优化。但从哪里开始呢？

为了确定首先定位哪些 LAPACK 例程，我解析了 LAPACK 的 Fortran 源代码以生成调用图。这使我能够推断出每个 LAPACK 例程的依赖关系树。有了这张图，我就进行了拓扑排序，从而帮助我识别没有依赖性的例程，并且这些例程将不可避免地成为其他例程的构建块。虽然我选择一个特定的高级例程并向后工作的深度优先方法将使我能够实现特定的功能，但这种方法可能会导致我在尝试实现日益复杂的例程时陷入困境。通过关注图表的“叶子”，我可以优先考虑常用例程（即具有高入度的例程），从而通过解锁稍后提供多个更高级别例程的能力来最大化我的影响我的努力或其他贡献者的努力。

有了我的计划，我很高兴能够开始工作。对于我的第一个例程，我选择了 dlaswp，它根据提供的主元索引列表在一般矩形矩阵上执行一系列行交换，并且它是 LAPACK 的 LU 分解例程的关键构建块。那就是我的挑战开始的时候......

挑战

旧版 Fortran

在 Quansight Labs 实习之前，我是（现在仍然是！）LFortran 的定期贡献者，LFortran 是一个基于 LLVM 构建的现代交互式 Fortran 编译器，我对自己的 Fortran 技能相当有信心。然而，我面临的首要挑战之一就是理解现在被认为是“遗留”的 Fortran 代码。我在下面强调了三个最初的障碍。

格式化

LAPACK 最初是用 FORTRAN 77 (F77) 编写的。虽然该库在版本 3.2 (2008) 中转移到 Fortran 90，但旧约定仍然保留在参考实现中。这些约定中最明显的一种是格式。

编写 F77 程序的开发人员使用从穿孔卡继承的固定表单布局来实现此目的。这种布局对字符列的使用有严格的要求：

• 占据整行的注释必须以第一列中的特殊字符（例如 *、! 或 C）开头。
• 对于非注释行，1) 前五列必须为空或包含数字标签，2) 第六列保留用于连续字符，3) 可执行语句必须从第七列开始，4) 超出的任何代码第 72 列被忽略。

Fortran 90 引入了自由形式布局，消除了列和行长度限制并确定了！作为评论字符。以下代码片段显示了 LAPACK 例程 dlacpy 的参考实现：

pip install numpy

下一个代码片段显示了相同的例程，但使用 Fortran 90 中引入的自由格式布局实现。

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

正如我们所观察到的，通过删除列限制并摆脱使用全部大写字母书写说明符的 F77 约定，现代 Fortran 代码更加明显一致，因此更具可读性。

标记控制结构

LAPACK 例程中的另一个常见做法是使用标记控制结构。例如，考虑以下代码片段，其中标签 10 必须与相应的 CONTINUE 匹配。

npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy

Fortran 90 消除了这种做法的需要，并通过允许使用 end do 来结束 do 循环来提高代码可读性。此更改显示在上面提供的 dlacpy 的自由形式版本中。

假定大小的数组

为了能够灵活地处理不同大小的数组，LAPACK 例程通常对具有假定大小的数组进行操作。在上面的 dlacpy 例程中，输入矩阵 A 被声明为具有根据表达式 A(LDA, *) 的假定大小的二维数组。该表达式声明 A 具有 LDA 行数，并使用 * 作为占位符来指示第二维的大小由调用程序确定。

使用假定大小数组的一个后果是编译器无法对未指定的维度执行边界检查。因此，当前的最佳实践是使用显式接口和假定形状数组（例如，A(LDA,:)）以防止越界内存访问。也就是说，当需要将子矩阵传递给其他函数时，使用假定形状数组可能会出现问题，因为这样做需要切片，这通常会导致编译器创建数组数据的内部副本。

迁移到 Fortran 95

不用说，我花了一段时间才适应 LAPACK 惯例并采用 LAPACK 思维方式。然而，作为一个纯粹主义者，如果我无论如何都要移植例程，我至少想将我设法移植的那些例程带入更现代的时代，以期提高代码的可读性和未来的维护。因此，在与 stdlib 维护者讨论后，我决定将例程迁移到 Fortran 95，虽然它不是最新、最好的 Fortran 版本，但似乎在保持原始实现的外观和感觉之间取得了适当的平衡，确保（足够好）向后兼容性，并利用更新的语法功能。

测试覆盖率

采用自下而上的方法来添加 LAPACK 支持的问题之一是，LAPACK 中通常不存在针对较低级别实用程序例程的显式单元测试。 LAPACK 的测试套件主要采用分层测试原理，其中假设测试较高级别的例程以确保其依赖的较低级别例程作为整个工作流程的一部分正常运行。虽然有人可能会说，针对较低级别的例程，重点关注集成测试而不是单元测试是合理的，因为为每个例程添加测试可能会增加 LAPACK 测试框架的维护负担和复杂性，但这意味着我们不能轻易依赖先前的测试单元测试的艺术，并且必须自己为每个较低级别的例程提供全面的独立单元测试。

文档

与测试覆盖率类似，在 LAPACK 本身之外，找到展示较低级别例程使用的真实记录示例具有挑战性。虽然 LAPACK 例程之前始终有一个文档注释，提供输入参数和可能的返回值的描述，但如果没有代码示例，则可视化和理解预期的输入和输出值可能具有挑战性，特别是在处理专门的矩阵时。虽然缺少单元测试和记录示例都不是世界末日，但这意味着向 stdlib 添加 LAPACK 支持将比我预期的更加困难。编写基准、测试、示例和文档只会需要更多的时间和精力，这可能会限制我在实习期间可以实施的例程数量。

内存布局

在线性存储器中存储矩阵元素时，有两种选择：连续存储列或连续存储行（见图 2）。前一种内存布局称为列优先顺序，后者称为行优先顺序。

图 2：演示以 (a) 列优先（Fortran 样式）或 (b) 行优先（C 样式）顺序将矩阵元素存储在线性存储器中的示意图。选择使用哪种布局很大程度上是一个约定问题。

选择使用哪种布局很大程度上是一个约定问题。例如，Fortran 以列优先顺序存储元素，C 以行优先顺序存储元素。更高级别的库，例如 NumPy 和 stdlib，支持列优先顺序和行优先顺序，允许您在数组创建期间配置多维数组的布局。

pip install numpy

虽然两种内存布局本质上都不比另一种更好，但根据底层存储模型的约定安排数据以确保顺序访问对于确保最佳性能至关重要。现代 CPU 能够比非顺序数据更有效地处理顺序数据，这主要归功于 CPU 缓存，而 CPU 缓存又利用了引用的空间局部性。

为了演示顺序与非顺序元素访问的性能影响，请考虑以下函数，该函数将所有元素从 MxN 矩阵 A 复制到另一个 MxN 矩阵 B，并且假设矩阵元素存储在列优先中订购。

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

设 A 和 B 为以下 3x2 矩阵：

A=⎣⎡​ 135​ 246​ ⎦ ⎤​, B=⎣ ⎡​000​ 000​⎦⎤ ​

当 A 和 B 都按列优先顺序存储时，我们可以按如下方式调用复制例程：

pip install numpy

但是，如果 A 和 B 都按行优先顺序存储，则调用签名将更改为

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

请注意，在后一种情况下，我们无法在最内层循环中按顺序访问元素，因为 da0 是 2，da1 是 -5，对于 db0 和 db1 也是如此。相反，数组索引“指针”在返回到线性内存中较早的元素之前会重复向前跳过，其中 ia = {0, 2, 4, 1, 3, 5} 和 ib 相同。在图 3 中，我们展示了非顺序访问对性能的影响。

图 3：当 copy 假定根据列优先顺序进行顺序元素访问时，向 copy 提供方列优先矩阵与行优先矩阵时的性能比较。 x 轴列举增加的矩阵大小（即元素数量）。所有比率均相对于相应矩阵大小的列主要结果进行归一化。

从图中，我们可以观察到列优先和行优先的性能大致相当，直到我们对具有超过 1e5 个元素 (M = N = ~316) 的方阵进行操作。对于 1e6 个元素 (M = N = ~1000)，提供行优先矩阵进行复制会导致性能下降 25% 以上。对于 1e7 元素 (M = N = ~3160)，我们观察到性能下降超过 85%。显着的性能影响可能归因于在具有大行大小的行主矩阵上操作时引用局部性的降低。

鉴于 LAPACK 是用 Fortran 编写的，它假定列优先访问顺序并相应地实现其算法。这给 stdlib 等库带来了问题，它们不仅支持行主序，而且将其设为默认内存布局。如果我们只是将 LAPACK 的 Fortran 实现移植到 JavaScript，那么提供行主矩阵的用户将会因非顺序访问而遭受不利的性能影响。

为了减轻不利的性能影响，我们借鉴了 BLIS（一个类似 BLAS 的库，支持 BLAS 例程中的行主内存布局和列主内存布局）的想法，并决定在将例程从 Fortran 移植到 JavaScript 时创建修改后的 LAPACK 实现， C 通过每个维度的单独步幅参数显式适应列主内存布局和行主内存布局。对于某些实现，例如与上面定义的复制函数类似的 dlacpy，合并单独且独立的步幅是简单的，通常涉及步幅技巧和循环交换，但是，对于其他实现，由于以下原因，修改结果并不那么简单。专门的矩阵处理、不同的访问模式和组合参数化。

ndarrays

LAPACK 例程主要对存储在线性存储器中的矩阵进行操作，并且根据指定的维度和前导（即第一个）维度的步长来访问其元素。维度分别指定每行和每列中的元素数量。步长指定必须跳过线性内存中的多少个元素才能访问行的下一个元素。 LAPACK 假设属于同一列的元素总是连续的（即在线性存储器中相邻）。图 4 提供了 LAPACK 约定的直观表示（特别是原理图 (a) 和 (b)）。

图 4：示意图说明了 LAPACK 跨步数组约定泛化为非连续跨步数组。 a) 以列优先顺序存储的 5×5 连续矩阵。 b) 以列优先顺序存储的 3×3 非连续子矩阵。通过提供指向第一个索引元素的指针并指定前导（即第一个）维度的步长，可以在 LAPACK 中对子矩阵进行操作。在这种情况下，尽管每列只有三个元素，但由于当存储为较大矩阵的一部分时，线性存储器中的子矩阵元素不连续，因此前导维度的步长为 5。在 LAPACK 中，尾随（即第二个）维度的步长始终假定为统一。 c) 以列主顺序存储的 3×3 非连续子矩阵，具有非单位步长，并将 LAPACK 步长约定推广到前导维度和尾随维度。这种泛化支撑了 stdlib 的多维数组（也称为“ndarrays”）。

NumPy 和 stdlib 等库概括了 LAPACK 的跨步数组约定以支持

1. 最后一个维度中的非单位步长（见图 4 (c)）。 LAPACK 假设矩阵的最后一个维度始终具有单位步长（即列中的元素连续存储在线性内存中）。
2. 任何维度的负步幅。 LAPACK 要求前导矩阵维度的步长为正。
3. 具有两个以上维度的多维数组。 LAPACK 仅明确支持跨步向量和（子）矩阵。

对最后一个维度中非单位步长的支持确保支持 O(1) 创建线性内存的非连续视图，而无需显式数据移动。这些视图通常称为“切片”。作为示例，请考虑以下代码片段，该代码片段使用 stdlib 提供的 API 创建此类视图。

pip install numpy

如果最后一个维度不支持非单位步幅，则不可能从表达式 x['::2,::2'] 返回视图，因为需要将所选元素复制到新的线性内存缓冲区以确保连续性。

图 5：示意图，说明如何使用步幅操作来创建存储在线性存储器中的矩阵元素的翻转和旋转视图。对于所有子原理图，步幅列为 [trailing_dimension,leading_dimension]。每个原理图隐含一个“偏移量”，它指示第一个索引元素的索引，这样，对于矩阵 A，元素 A_ij 是根据 i⋅strides[1] j⋅strides[0] 偏移量解析。 a) 给定一个以列优先顺序存储的 3×3 矩阵，可以操纵前导维度和尾随维度的步幅来创建视图，在该视图中以相反的顺序访问沿一个或多个轴的矩阵元素。 b) 使用类似的步幅操作，可以创建矩阵元素相对于它们在线性内存中的排列的旋转视图。

对负步幅的支持使得元素能够沿一个或多个维度进行 O(1) 反转和旋转（参见图 5）。例如，要从上到下和从左到右翻转矩阵，只需取消步幅即可。在前面的代码片段的基础上，以下代码片段演示了围绕一个或多个轴反转元素。

pip install numpy

负步幅的讨论中隐含的是需要一个“偏移”参数，该参数指示线性内存中第一个索引元素的索引。对于跨步多维数组 A 和跨步列表 s，元素 A_ij⋅⋅⋅n 对应的索引可以根据方程求解

$$\text{idx}=\text{偏移量}i\cdot s_{0}j\cdot s_1\dots n\cdot s_{N-1}$$idx=偏移 i⋅s0​ j⋅s1​ … n⋅sN−1 ​

其中 N 是数组维度数，s_k 对应于第 k 个步幅。

在支持负步幅的 BLAS 和 LAPACK 例程中（仅在对步幅向量进行操作时才支持这种功能（例如，参见上面的 daxpy）），索引偏移量是使用类似于以下代码片段的逻辑来计算的：

pip install numpy

其中 M 是向量元素的数量。这隐含地假设提供的数据指针指向向量的线性存储器的开头。在支持指针的语言中，例如C，为了在线性存储器的不同区域上进行操作，通常在函数调用之前使用指针算术来调整指针，这相对便宜且简单，至少对于一维情况是这样。

例如，返回到上面定义的 c_daxpy，我们可以使用指针算术来限制对线性内存中从输入和输出数组的第十一个和第十六个元素（注意：从零开始索引）开始的五个元素的元素访问，分别如下面的代码片段所示。

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

但是，在 JavaScript 中，不支持二进制缓冲区的显式指针算术，必须显式实例化具有所需字节偏移量的新类型数组对象。在下面的代码片段中，为了获得与上面的 C 示例相同的结果，我们必须解析一个类型化数组构造函数，计算一个新的字节偏移量，计算一个新的类型化数组长度，并创建一个新的类型化数组实例。

npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy

对于大型数组，与访问和操作单个数组元素所花费的时间相比，类型化数组实例化的成本可以忽略不计；然而，对于较小的数组大小，对象实例化会显着影响性能。

因此，为了避免对对象实例化性能产生不利影响，stdlib 将 ndarray 的数据缓冲区与 ndarray 视图开头对应的缓冲区元素的位置解耦。这允许切片表达式 x[2:,3:] 和 x[3:,1:] 返回新的 ndarray 视图而无需需要实例化新的缓冲区实例，如以下代码片段所示。

// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';

// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;

// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');

// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);

由于将数据缓冲区与 ndarray 视图的开头解耦，我们同样试图避免在使用 ndarray 数据调用 LAPACK 例程时实例化新的类型化数组实例。这意味着创建修改后的 LAPACK API 签名，支持所有跨步向量和矩阵的显式偏移参数。

为了简单起见，让我们回到上面定义的 daxpy 的 JavaScript 实现。

pip install numpy

如以下代码片段所示，我们可以修改上述签名和实现，以便将解析第一个索引元素的责任转移给 API 使用者。

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

对于 ndarray，解析发生在 ndarray 实例化期间，使得使用 ndarray 数据调用 daxpy_ndarray 可以直接传递关联的 ndarray 元数据。下面的代码片段对此进行了演示。

npm install @stdlib/ndarray-fancy @stdlib/blas-daxpy

与 BLIS 类似，我们看到了传​​统 LAPACK API 签名（例如，为了向后兼容）和修改后的 API 签名（例如，为了最大限度地减少不利的性能影响）的价值，因此，我们制定了一项计划，提供传统和修改了每个 LAPACK 例程的 API。为了最大限度地减少代码重复，我们的目标是实现一个通用的较低级别的“基本”实现，然后可以由较高级别的 API 进行包装。虽然上面显示的 BLAS 例程 daxpy 的更改可能看起来相对简单，但传统 LAPACK 例程及其预期行为到通用实现的转换通常要复杂得多。

德拉斯普

挑战够了！最终产品是什么样的？！

让我们绕一圈，回到 dlaswp，这是一个 LAPACK 例程，用于根据主元索引列表在输入矩阵上执行一系列行交换。以下代码片段显示了参考 LAPACK Fortran 实现。

// Individually import desired functionality:
import FancyArray from '@stdlib/ndarray-fancy';
import daxpy from '@stdlib/blas-daxpy';

// Define ndarray meta data:
const shape = [4, 4, 4];
const strides = [...];
const offset = 0;

// Define arrays using a "lower-level" fancy array constructor:
const x = new FancyArray('float64', [...], shape, strides, offset, 'row-major');
const y = new FancyArray('float64', [...], shape, strides, offset, 'row-major');

// Perform operation:
daxpy(5.0, x['::2,:,:'], y['::2,:,:']);

为了方便与 C 语言的 Fortran 实现进行交互，LAPACK 提供了一个名为 LAPACKE 的两级 C 接口，它包装了 Fortran 实现，并为行主和列主输入和输出矩阵提供了便利。 dlaswp 的中层接口如以下代码片段所示。

void c_daxpy(const int N, const double alpha, const double *X, const int strideX, double *Y, const int strideY) {
    int ix;
    int iy;
    int i;
    if (N <= 0) {
        return;
    }
    if (alpha == 0.0) {
        return;
    }
    if (strideX < 0) {
        ix = (1-N) * strideX;
    } else {
        ix = 0;
    }
    if (strideY < 0) {
        iy = (1-N) * strideY;
    } else {
        iy = 0;
    }
    for (i = 0; i < N; i++) {
        Y[iy] += alpha * X[ix];
        ix += strideX;
        iy += strideY;
    }
    return;
}

当使用列主矩阵 a 调用时，包装器 LAPACKE_dlaswp_work 只是将提供的参数传递给 Fortran 实现。但是，当使用行主矩阵 a 调用时，包装器必须分配内存、显式转置并将 a 复制到临时矩阵 a_t、重新计算前导维度的步幅、使用 a_t 调用 dlaswp、转置并复制存储在 a_t 中的结果到a，最后释放分配的内存。这是相当大的工作量，并且在大多数 LAPACK 例程中很常见。

以下代码片段显示了移植到 JavaScript 的参考 LAPACK 实现，支持前导和尾随维度步幅、索引偏移以及包含主元索引的步幅向量。

const binary = new UintArray([...]);

const mod = new WebAssembly.Module(binary);

为了提供与传统 LAPACK 具有一致行为的 API，我包装了上述实现并将输入参数调整为“基本”实现，如以下代码片段所示。

pip install numpy

我随后编写了一个单独但类似的包装器，它提供了更直接到 stdlib 多维数组的 API 映射，并在应用枢轴的方向为负时执行一些特殊处理，如以下代码片段所示。

# Import all of NumPy:
import numpy as np

# Define arrays:
x = np.asarray(...)
y = np.asarray(...)

# Perform operation:
y[::2,:,:] += 5.0 * x[::2,:,:]

需要注意的几点：

1. 与传统的 LAPACKE API 相比，在 dlaswp_ndarray 和基础 API 中，matrix_layout（顺序）参数不是必需的，因为可以从提供的步幅推断顺序。
2. 与传统的LAPACKE API相比，当输入矩阵为行主矩阵时，我们不需要将数据复制到临时工作区数组，从而减少了不必要的内存分配。
3. 与直接与 BLAS 和 LAPACK 交互的 NumPy 和 SciPy 等库相比，在 stdlib 中调用 LAPACK 例程时，我们不需要之前将不连续的多维数据复制到临时工作区数组或从临时工作区数组复制非连续的多维数据和调用后，分别。除了与硬件优化的 BLAS 和 LAPACK 接口时之外，所采用的方法有助于最大限度地减少数据移动并确保资源受限的浏览器应用程序的性能。

对于希望利用硬件优化库（例如 OpenBLAS）的服务器端应用程序，我们提供单独的包装器，将通用签名参数适应其优化的 API 等效项。在这种情况下，至少对于足够大的数组来说，创建临时副本是值得的。

当前状态和后续步骤

尽管面临挑战、不可预见的挫折和多次设计迭代，我很高兴地报告，除了上面的 dlaswp 之外，我还能够打开 35 个 PR，添加对各种 LAPACK 例程和相关实用程序的支持。显然不完全是 1,700 个例程，但这是一个好的开始！ :)

尽管如此，未来是光明的，我们对这项工作感到非常兴奋。仍有很大的改进空间和额外的研究和开发。我们尤其渴望

1. 探索工具和自动化。
2. 解决了解析跨多个 stdlib 包的 Fortran 依赖项的源文件时的构建问题。
3. 为 stdlib 的现有 LAPACK 包推出 C 和 Fortran 实现以及本机绑定。
4. 继续发展 stdlib 的模块化 LAPACK 例程库。
5. 确定性能优化的其他领域。

虽然我的 Quansight Labs 实习已经结束，但我的计划是继续添加软件包并推动这项工作。鉴于 LAPACK 的巨大潜力和根本重要性，我们很高兴看到将 LAPACK 引入网络的这一举措继续发展，因此，如果您有兴趣帮助推动这一进程，请随时与我们联系！如果您有兴趣赞助开发，Quansight 的人员将非常乐意与您聊天。

在此，我要感谢 Quansight 提供的这次实习机会。我感到非常幸运能够学到这么多东西。在 Quansight 实习是我长期以来的一个梦想，我很高兴能够实现这个梦想。我要特别感谢 Athan Reines 和 Melissa Mendonça，他们是一位出色的导师，也是一位出色的人！感谢所有 stdlib 核心开发人员以及 Quansight 的其他所有人，感谢你们一路上大大小小的帮助了我。

干杯！

---

stdlib 是一个开源软件项目，致力于提供一整套强大、高性能的库来加速您的项目开发，并让您安心，因为您知道您依赖的是精心制作的高质量软件。

如果您喜欢这篇文章，请给我们一颗星？ GitHub 上并考虑为该项目提供经济支持。您的贡献和持续的支持有助于确保项目的长期成功，我们深表感谢！

以上是网络浏览器中的 LAPACK的详细内容。更多信息请关注PHP中文网其他相关文章！