這篇文章最初發(fā)布在 Quansight Labs 博客上，經(jīng) Quansight 許可修改并重新發(fā)布于此處。

Web 應(yīng)用程序正在迅速崛起，成為高性能科學(xué)計(jì)算和支持 AI 的最終用戶體驗(yàn)的新領(lǐng)域。支撐機(jī)器學(xué)習(xí)/人工智能革命的是??線性代數(shù)，它是關(guān)于線性方程及其在向量空間和矩陣中的表示的數(shù)學(xué)分支。 LAPACK（“Linear Algebra Package”）是數(shù)值線性代數(shù)的基礎(chǔ)軟件庫，提供強(qiáng)大的、經(jīng)過實(shí)戰(zhàn)檢驗(yàn)的常見矩陣運(yùn)算實(shí)現(xiàn)。盡管 LAPACK 是大多數(shù)數(shù)值計(jì)算編程語言和庫的基礎(chǔ)組件，但針對(duì)網(wǎng)絡(luò)獨(dú)特限制量身定制的全面、高質(zhì)量的 LAPACK 實(shí)現(xiàn)尚未實(shí)現(xiàn)。那是……到現(xiàn)在為止。

今年早些時(shí)候，我有幸成為 Quansight Labs 的暑期實(shí)習(xí)生，該實(shí)驗(yàn)室是 Quansight 的公益部門，也是科學(xué) Python 生態(tài)系統(tǒng)的領(lǐng)導(dǎo)者。在實(shí)習(xí)期間，我致力于為 stdlib 添加初始 LAPACK 支持，stdlib 是一個(gè)用 C 和 JavaScript 編寫的科學(xué)計(jì)算基礎(chǔ)庫，并針對(duì)在 Web 瀏覽器和其他 Web 原生環(huán)境（例如 Node.js 和 Deno）中使用進(jìn)行了優(yōu)化。在這篇博文中，我將討論我的旅程、一些預(yù)期和意外（！）的挑戰(zhàn)以及未來的道路。我希望，如果運(yùn)氣好的話，這項(xiàng)工作能夠?yàn)槭咕W(wǎng)絡(luò)瀏覽器成為數(shù)值計(jì)算和機(jī)器學(xué)習(xí)的一流環(huán)境提供關(guān)鍵的構(gòu)建模塊，并預(yù)示著更強(qiáng)大的人工智能網(wǎng)絡(luò)應(yīng)用程序的未來。

聽起來很有趣嗎？我們走吧！

什么是 stdlib？

熟悉 LAPACK 的本博客讀者可能不太熟悉 Web 技術(shù)的狂野世界。對(duì)于那些來自數(shù)值和科學(xué)計(jì)算領(lǐng)域并熟悉科學(xué) Python 生態(tài)系統(tǒng)的人來說，最簡(jiǎn)單的方式是將 stdlib 視為 NumPy 和 SciPy 模型中的開源科學(xué)計(jì)算庫。它提供多維數(shù)組數(shù)據(jù)結(jié)構(gòu)以及數(shù)學(xué)、統(tǒng)計(jì)和線性代數(shù)的相關(guān)例程，但使用 JavaScript 而不是 Python 作為其主要腳本語言。因此，stdlib 專注于 Web 生態(tài)系統(tǒng)及其應(yīng)用程序開發(fā)范例。這種關(guān)注需要一些有趣的設(shè)計(jì)和項(xiàng)目架構(gòu)決策，這使得 stdlib 與為數(shù)值計(jì)算設(shè)計(jì)的更傳統(tǒng)的庫相比相當(dāng)獨(dú)特。

以 NumPy 為例，NumPy 是一個(gè)單一的整體庫，其中所有組件（除了可選的第三方依賴項(xiàng)（例如 OpenBLAS）之外）形成一個(gè)不可分割的單元。如果不安裝所有 NumPy，就無法簡(jiǎn)單地安裝 NumPy 例程來進(jìn)行數(shù)組操作。如果您正在部署一個(gè)僅需要 NumPy 的 ndarray 對(duì)象及其一些操作例程的應(yīng)用程序，那么安裝和捆綁所有 NumPy 意味著包含大量“死代碼”。用 Web 開發(fā)的術(shù)語來說，我們會(huì)說 NumPy 不是“tree shakeable”。對(duì)于正常的 NumPy 安裝，這意味著至少需要 30MB 的磁盤空間，對(duì)于排除所有調(diào)試語句的自定義構(gòu)建至少需要 15MB 的磁盤空間。對(duì)于 SciPy，這些數(shù)字可以分別增加到 130MB 和 50MB。不用說，在 Web 應(yīng)用程序中僅提供幾個(gè)功能的 15MB 庫是不可能的，特別是對(duì)于需要將 Web 應(yīng)用程序部署到網(wǎng)絡(luò)連接較差或內(nèi)存受限的設(shè)備的開發(fā)人員而言。

考慮到 Web 應(yīng)用程序開發(fā)的獨(dú)特限制，stdlib 采用自下而上的設(shè)計(jì)方法，其中每個(gè)功能單元都可以獨(dú)立于代碼庫中不相關(guān)和未使用的部分進(jìn)行安裝和使用。通過采用可分解的軟件架構(gòu)和徹底的模塊化，stdlib 使用戶能夠安裝和使用他們所需要的內(nèi)容，除了所需的 API 集及其顯式依賴項(xiàng)之外幾乎沒有多余的代碼，從而確保更小的內(nèi)存占用、捆綁尺寸和更快的部署。

舉個(gè)例子，假設(shè)您正在使用兩個(gè)矩陣堆棧（即三維立方體的二維切片），并且您想要選擇每隔一個(gè)切片并執(zhí)行常見的 BLAS 運(yùn)算 y = a * x，其中 x 和 y 是 ndarray，a 是標(biāo)量常量。要使用 NumPy 執(zhí)行此操作，您首先需要安裝所有 NumPy

pip install numpy

然后執(zhí)行各種操作

# 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，除了能夠?qū)㈨?xiàng)目安裝為整體庫之外，您還可以將各種功能單元安裝為單獨(dú)的包

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

然后執(zhí)行各種操作

// 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,:,:']);

重要的是，您不僅可以獨(dú)立安裝 stdlib 超過 4,000 個(gè)軟件包中的任何一個(gè)，還可以通過分叉關(guān)聯(lián)的 GitHub 存儲(chǔ)庫來修復(fù)、改進(jìn)和重新混合其中任何一個(gè)軟件包（例如，請(qǐng)參閱 @stdlib/ndarray-fancy） ）。通過定義顯式的抽象層和依賴樹，stdlib 使您可以自由地為您的應(yīng)用程序選擇正確的抽象層。在某些方面，這是一個(gè)簡(jiǎn)單的想法，如果您習(xí)慣了傳統(tǒng)的科學(xué)軟件庫設(shè)計(jì)，也許是非正統(tǒng)的想法，但是，當(dāng)與 Web 平臺(tái)緊密集成時(shí)，它會(huì)產(chǎn)生強(qiáng)大的后果并創(chuàng)造令人興奮的新可能性！

WebAssembly 怎么樣？

好吧，也許你的興趣已經(jīng)激起； stdlib 似乎很有趣。但這與網(wǎng)絡(luò)瀏覽器中的 LAPACK 有什么關(guān)系呢？嗯，去年夏天我們的目標(biāo)之一是應(yīng)用 stdlib 精神——小型、范圍狹窄的軟件包，只做一件事，并做好一件事——將 LAPACK 引入網(wǎng)絡(luò)。

但是等等，你說！這是一項(xiàng)極端的任務(wù)。 LAPACK 規(guī)模龐大，約有 1,700 個(gè)例程，在合理的時(shí)間范圍內(nèi)實(shí)施其中的 10% 都是一項(xiàng)重大挑戰(zhàn)。將 LAPACK 編譯為 WebAssembly 不是更好嗎？WebAssembly 是 C、Go 和 Rust 等編程語言的可移植編譯目標(biāo)，可以在網(wǎng)絡(luò)上部署，然后就到此為止了？

不幸的是，這種方法存在幾個(gè)問題。

1. 將 Fortran 編譯為 WebAssembly 仍然是一個(gè)正在積極開發(fā)的領(lǐng)域（請(qǐng)參閱 1、2、3、4 和 5）。在撰寫本文時(shí)，常見的方法是使用 f2c 將 Fortran 編譯為 C，然后執(zhí)行單獨(dú)的編譯步驟將 C 轉(zhuǎn)換為 WebAssembly。然而，這種方法是有問題的，因?yàn)?f2c 僅完全支持 Fortran 77，并且生成的代碼需要大量修補(bǔ)。開發(fā)基于 LLVM 的 Fortran 編譯器的工作正在進(jìn)行中，但仍然存在差距和復(fù)雜的工具鏈。
2. 正如上面關(guān)于 Web 應(yīng)用程序中的整體庫的討論中提到的，LAPACK 的龐大是問題的一部分。即使編譯問題得到解決，在僅需要使用一兩個(gè) LAPACK 例程的 Web 應(yīng)用程序中包含包含所有 LAPACK 的單個(gè) WebAssembly 二進(jìn)制文件也意味著大量死代碼，導(dǎo)致加載時(shí)間變慢并增加內(nèi)存消耗。
3. 雖然人們可以嘗試將單個(gè) LAPACK 例程編譯為獨(dú)立的 WebAssembly 二進(jìn)制文件，但這樣做可能會(huì)導(dǎo)致二進(jìn)制膨脹，因?yàn)槎鄠€(gè)獨(dú)立二進(jìn)制文件可能包含來自公共依賴項(xiàng)的重復(fù)代碼。為了減輕二進(jìn)制膨脹，可以嘗試執(zhí)行模塊拆分。在這種情況下，首先將公共依賴項(xiàng)分解為包含共享代碼的獨(dú)立二進(jìn)制文件，然后為各個(gè) API 生成單獨(dú)的二進(jìn)制文件。雖然在某些情況下適用，但這種方法很快就會(huì)變得笨拙，因?yàn)檫@種方法需要在加載時(shí)通??過將一個(gè)或多個(gè)模塊的導(dǎo)出與一個(gè)或多個(gè)其他模塊的導(dǎo)入縫合在一起來鏈接各個(gè) WebAssembly 模塊。這不僅很乏味，而且這種方法還會(huì)帶來性能損失，因?yàn)楫?dāng) WebAssembly 例程調(diào)用導(dǎo)入的導(dǎo)出時(shí)，它們現(xiàn)在必須跨入 JavaScript，而不是保留在 WebAssembly 中。聽起來很復(fù)雜？是的！
4. 除了專門對(duì)標(biāo)量輸入?yún)?shù)進(jìn)行操作的 WebAssembly 模塊（例如，計(jì)算單個(gè)數(shù)字的正弦值）之外，每個(gè) WebAssembly 模塊實(shí)例都必須與 WebAssembly 內(nèi)存關(guān)聯(lián)，該內(nèi)存以 64KiB 的固定增量分配（即“頁面”） ”）。重要的是，截至這篇博文，WebAssembly 內(nèi)存只能增長而永遠(yuǎn)不會(huì)收縮。由于當(dāng)前沒有向主機(jī)釋放內(nèi)存的機(jī)制，因此 WebAssembly 應(yīng)用程序的內(nèi)存占用只會(huì)增加。這兩個(gè)方面的結(jié)合增加了分配從未使用過的內(nèi)存的可能性以及內(nèi)存泄漏的普遍性。
5. 最后，雖然 WebAssembly 功能強(qiáng)大，但它需要更陡峭的學(xué)習(xí)曲線和一組更復(fù)雜且經(jīng)?？焖侔l(fā)展的工具鏈。在最終用戶應(yīng)用程序中，JavaScript（一種 Web 原生動(dòng)態(tài)編譯編程語言）和 WebAssembly 之間的接口進(jìn)一步增加了復(fù)雜性，尤其是在必須執(zhí)行手動(dòng)內(nèi)存管理時(shí)。

為了幫助說明最后一點(diǎn)，讓我們回到 BLAS 例程 daxpy，它執(zhí)行運(yùn)算 y = a*x y，其中 x 和 y 是跨步向量，a 是標(biāo)量常量。如果用 C 實(shí)現(xiàn)，基本實(shí)現(xiàn)可能類似于以下代碼片段。

pip install numpy

編譯到 WebAssembly 并將 WebAssembly 二進(jìn)制文件加載到我們的 Web 應(yīng)用程序中后，我們需要執(zhí)行一系列步驟，然后才能從 JavaScript 調(diào)用 c_daxpy 例程。首先，我們需要實(shí)例化一個(gè)新的 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,:,:]

接下來，我們需要定義模塊內(nèi)存并創(chuàng)建一個(gè)新的 WebAssembly 模塊實(shí)例。

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

創(chuàng)建模塊實(shí)例后，我們現(xiàn)在可以調(diào)用導(dǎo)出的 BLAS 例程。但是，如果數(shù)據(jù)是在模塊內(nèi)存之外定義的，我們首先需要將該數(shù)據(jù)復(fù)制到內(nèi)存實(shí)例，并且始終以小端字節(jié)順序執(zhí)行此操作。

// 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,:,:']);

現(xiàn)在數(shù)據(jù)已寫入模塊內(nèi)存，我們可以調(diào)用 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;
}

最后，如果我們需要將結(jié)果傳遞到不支持 WebAssembly 內(nèi)存“指針”（即字節(jié)偏移）的下游庫（例如 D3）以進(jìn)行可視化或進(jìn)一步分析，我們需要從模塊復(fù)制數(shù)據(jù)內(nèi)存回原始輸出數(shù)組。

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

const mod = new WebAssembly.Module(binary);

僅僅計(jì)算 y = a*x y 就需要大量工作。相反，與純 JavaScript 實(shí)現(xiàn)相比，它可能類似于以下代碼片段。

// 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 實(shí)現(xiàn)，我們可以使用外部定義的數(shù)據(jù)直接調(diào)用 daxpy，而無需上面 WebAssembly 示例中所需的數(shù)據(jù)移動(dòng)。

// 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 方法不僅不太符合人體工程學(xué)，而且正如所預(yù)期的，考慮到所需的數(shù)據(jù)移動(dòng)，還會(huì)對(duì)性能產(chǎn)生負(fù)面影響，如下圖所示。

圖 1：對(duì)于增加數(shù)組長度（x 軸）的 BLAS 例程 daxpy，stdlib 的 C、JavaScript 和 WebAssembly (Wasm) 實(shí)現(xiàn)的性能比較。在Wasm（復(fù)制）基準(zhǔn)測(cè)試中，輸入和輸出數(shù)據(jù)被復(fù)制到Wasm內(nèi)存中或從Wasm內(nèi)存中復(fù)制，導(dǎo)致性能較差。

在上圖中，我顯示了 BLAS 例程 daxpy 的 stdlib 的 C、JavaScript 和 WebAssembly (Wasm) 實(shí)現(xiàn)的性能比較，以增加數(shù)組長度（沿 x 軸枚舉）。 y 軸顯示相對(duì)于基線 C 實(shí)現(xiàn)的標(biāo)準(zhǔn)化速率。在 Wasm 基準(zhǔn)測(cè)試中，輸入和輸出數(shù)據(jù)直接在 WebAssembly 模塊內(nèi)存中分配和操作，并且在 Wasm（復(fù)制）基準(zhǔn)測(cè)試中，輸入和輸出數(shù)據(jù)復(fù)制到 WebAssembly 模塊內(nèi)存中或從 WebAssembly 模塊內(nèi)存中復(fù)制，如上所述。從圖表中，我們可以觀察到以下幾點(diǎn)：

1. 一般來說，得益于高度優(yōu)化的即時(shí) (JIT) 編譯器，JavaScript 代碼在精心編寫的情況下，執(zhí)行速度僅比本機(jī)代碼慢 2 到 3 倍。對(duì)于松散類型、動(dòng)態(tài)編譯的編程語言來說，這個(gè)結(jié)果令人印象深刻，并且至少對(duì)于 daxpy 來說，在不同的數(shù)組長度上保持一致。
2. 隨著數(shù)據(jù)大小以及 WebAssembly 模塊中花費(fèi)的時(shí)間量的增加，WebAssembly 可以接近原生（~1.5 倍）速度。此結(jié)果更符合預(yù)期的 WebAssembly 性能。
3. 雖然 WebAssembly 可以實(shí)現(xiàn)接近本機(jī)的速度，但數(shù)據(jù)移動(dòng)要求可能會(huì)對(duì)性能產(chǎn)生不利影響，正如在 daxpy 中觀察到的那樣。在這種情況下，精心設(shè)計(jì)的 JavaScript 實(shí)現(xiàn)可以避免此類要求，即使不是更好，也可以實(shí)現(xiàn)相同的性能，就像 daxpy 的情況一樣。

總體而言，WebAssembly 可以提供性能改進(jìn)；然而，該技術(shù)并不是靈丹妙藥，需要謹(jǐn)慎使用才能實(shí)現(xiàn)預(yù)期收益。即使提供卓越的性能，這種收益也必須與復(fù)雜性增加、潛在的更大的捆綁包大小和更復(fù)雜的工具鏈的成本相平衡。對(duì)于許多應(yīng)用程序來說，簡(jiǎn)單的 JavaScript 實(shí)現(xiàn)就可以了。

徹底的模塊化

現(xiàn)在我已經(jīng)起訴了反對(duì)將整個(gè) LAPACK 編譯為 WebAssembly 的案件并就此結(jié)束，那我們?cè)撛趺崔k？好吧，如果我們要擁抱 stdlib 精神，那么我們就需要徹底的模塊化。

擁抱徹底的模塊化意味著要認(rèn)識(shí)到最好的東西是高度上下文相關(guān)的，并且根據(jù)用戶應(yīng)用程序的需求和限制，開發(fā)人員需要靈活地選擇正確的抽象。如果開發(fā)人員正在編寫 Node.js 應(yīng)用程序，這可能意味著綁定到硬件優(yōu)化庫，例如 OpenBLAS、Intel MKL 或 Apple Accelerate，以獲得卓越的性能。如果開發(fā)人員正在部署需要一小組數(shù)字例程的 Web 應(yīng)用程序，那么 JavaScript 可能是適合該工作的工具。如果開發(fā)人員正在開發(fā)大型的資源密集型 WebAssembly 應(yīng)用程序（例如，用于圖像編輯或游戲引擎），那么能夠輕松編譯單個(gè)例程作為大型應(yīng)用程序的一部分將至關(guān)重要。簡(jiǎn)而言之，我們想要一個(gè)完全模塊化的 LAPACK。

我的任務(wù)是為這樣的努力奠定基礎(chǔ)，解決問題并找到差距，并希望讓我們更接近網(wǎng)絡(luò)上的高性能線性代數(shù)。但是徹底的模塊化是什么樣的呢？這一切都始于基本的功能單元，包??。

stdlib 中的每個(gè)包都是獨(dú)立的，包含共同本地化的測(cè)試、基準(zhǔn)測(cè)試、示例、文檔、構(gòu)建文件和關(guān)聯(lián)的元數(shù)據(jù)（包括任何依賴項(xiàng)的枚舉），并與外界定義清晰的 API 界面。為了向 stdlib 添加 LAPACK 支持，這意味著為每個(gè) LAPACK 例程創(chuàng)建一個(gè)單獨(dú)的獨(dú)立包，其結(jié)構(gòu)如下：

pip install numpy

簡(jiǎn)單來說，

• benchmark：包含微基準(zhǔn)的文件夾，用于評(píng)估相對(duì)于參考實(shí)現(xiàn)（即參考 LAPACK）的性能。
• docs：包含輔助文檔的文件夾，其中包括 REPL 幫助文本和定義類型化 API 簽名的 TypeScript 聲明。
• examples：包含可執(zhí)行演示代碼的文件夾，除了充當(dāng)文檔之外，還可以幫助開發(fā)人員健全性檢查實(shí)現(xiàn)行為。
• include：包含 C 頭文件的文件夾。
• lib：包含 JavaScript 源實(shí)現(xiàn)的文件夾，其中 index.js 作為包入口點(diǎn)，其他 *.js 文件定義內(nèi)部實(shí)現(xiàn)模塊。
• src：包含 C 和 Fortran 源實(shí)現(xiàn)的文件夾。每個(gè)模塊化 LAPACK 包都應(yīng)包含稍作修改的 Fortran 參考實(shí)現(xiàn)（F77 到自由格式的 Fortran）。 C 文件包括遵循 Fortran 參考實(shí)現(xiàn)的純 C 實(shí)現(xiàn)、用于調(diào)用 Fortran 參考實(shí)現(xiàn)的包裝器、用于在服務(wù)器端應(yīng)用程序中調(diào)用硬件優(yōu)化庫（例如 OpenBLAS）的包裝器，以及用于調(diào)用已編譯的本機(jī)綁定來自 Node.js 中的 JavaScript 或兼容的服務(wù)器端 JavaScript 運(yùn)行時(shí)的 C。
• test：包含單元測(cè)試的文件夾，用于測(cè)試 JavaScript 和本機(jī)實(shí)現(xiàn)中的預(yù)期行為。本機(jī)實(shí)現(xiàn)的測(cè)試是用 JavaScript 編寫的，并利用本機(jī)綁定實(shí)現(xiàn) JavaScript 和 C/Fortran 之間的互操作。
• binding.gyp/include.gypi：用于編譯 Node.js 原生插件的構(gòu)建文件，它提供了 JavaScript 和原生代碼之間的橋梁。
• manifest.json：stdlib內(nèi)部C和Fortran編譯源文件包管理的配置文件。
• package.json：包含包元數(shù)據(jù)的文件，包括外部包依賴項(xiàng)的枚舉以及用于基于瀏覽器的 Web 應(yīng)用程序的純 JavaScript 實(shí)現(xiàn)的路徑。
• README.md：包含包的主要文檔的文件，其中包括 API 簽名以及 JavaScript 和 C 接口的示例。

考慮到 stdlib 苛刻的文檔和測(cè)試要求，為每個(gè)例程添加支持是一項(xiàng)相當(dāng)大的工作量，但最終結(jié)果是健壯的、高質(zhì)量的，最重要的是，適合作為科學(xué)計(jì)算基礎(chǔ)的模塊化代碼在現(xiàn)代網(wǎng)絡(luò)上。但預(yù)賽就夠了！讓我們言歸正傳吧！

多階段方法

基于之前為 stdlib 添加 BLAS 支持的努力，我們決定在添加 LAPACK 支持時(shí)遵循類似的多階段方法，其中我們首先優(yōu)先考慮 JavaScript 實(shí)現(xiàn)及其相關(guān)的測(cè)試和文檔，然后，一旦存在測(cè)試和文檔、回填 C 和 Fortran 實(shí)現(xiàn)以及任何與硬件優(yōu)化庫相關(guān)的本機(jī)綁定。這種方法使我們能夠在董事會(huì)上提出一些早期觀點(diǎn)，可以這么說，快速將 API 提供給用戶，建立強(qiáng)大的測(cè)試程序和基準(zhǔn)，并在深入研究構(gòu)建工具鏈和性能之前調(diào)查工具和自動(dòng)化的潛在途徑優(yōu)化。但從哪里開始呢？

為了確定首先定位哪些 LAPACK 例程，我解析了 LAPACK 的 Fortran 源代碼以生成調(diào)用圖。這使我能夠推斷出每個(gè) LAPACK 例程的依賴關(guān)系樹。有了這張圖，我就進(jìn)行了拓?fù)渑判?，從而幫助我識(shí)別沒有依賴性的例程，并且這些例程將不可避免地成為其他例程的構(gòu)建塊。雖然我選擇一個(gè)特定的高級(jí)例程并向后工作的深度優(yōu)先方法將使我能夠?qū)崿F(xiàn)特定的功能，但這種方法可能會(huì)導(dǎo)致我在嘗試實(shí)現(xiàn)日益復(fù)雜的例程時(shí)陷入困境。通過關(guān)注圖表的“葉子”，我可以優(yōu)先考慮常用例程（即具有高入度的例程），從而通過解鎖稍后提供多個(gè)更高級(jí)別例程的能力來最大化我的影響我的努力或其他貢獻(xiàn)者的努力。

有了我的計(jì)劃，我很高興能夠開始工作。對(duì)于我的第一個(gè)例程，我選擇了 dlaswp，它根據(jù)提供的主元索引列表在一般矩形矩陣上執(zhí)行一系列行交換，并且它是 LAPACK 的 LU 分解例程的關(guān)鍵構(gòu)建塊。那就是我的挑戰(zhàn)開始的時(shí)候......

挑戰(zhàn)

舊版 Fortran

在 Quansight Labs 實(shí)習(xí)之前，我是（現(xiàn)在仍然是?。㎜Fortran 的定期貢獻(xiàn)者，LFortran 是一個(gè)基于 LLVM 構(gòu)建的現(xiàn)代交互式 Fortran 編譯器，我對(duì)自己的 Fortran 技能相當(dāng)有信心。然而，我面臨的首要挑戰(zhàn)之一就是理解現(xiàn)在被認(rèn)為是“遺留”的 Fortran 代碼。我在下面強(qiáng)調(diào)了三個(gè)最初的障礙。

格式化

LAPACK 最初是用 FORTRAN 77 (F77) 編寫的。雖然該庫在版本 3.2 (2008) 中轉(zhuǎn)移到 Fortran 90，但舊約定仍然保留在參考實(shí)現(xiàn)中。這些約定中最明顯的一種是格式。

編寫 F77 程序的開發(fā)人員使用從穿孔卡繼承的固定表單布局來實(shí)現(xiàn)此目的。這種布局對(duì)字符列的使用有嚴(yán)格的要求：

• 占據(jù)整行的注釋必須以第一列中的特殊字符（例如 *、! 或 C）開頭。
• 對(duì)于非注釋行，1) 前五列必須為空或包含數(shù)字標(biāo)簽，2) 第六列保留用于連續(xù)字符，3) 可執(zhí)行語句必須從第七列開始，4) 超出的任何代碼第 72 列被忽略。

Fortran 90 引入了自由形式布局，消除了列和行長度限制并確定了！作為評(píng)論字符。以下代碼片段顯示了 LAPACK 例程 dlacpy 的參考實(shí)現(xiàn)：

pip install numpy

下一個(gè)代碼片段顯示了相同的例程，但使用 Fortran 90 中引入的自由格式布局實(shí)現(xiàn)。

# 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 約定，現(xiàn)代 Fortran 代碼更加明顯一致，因此更具可讀性。

標(biāo)記控制結(jié)構(gòu)

LAPACK 例程中的另一個(gè)常見做法是使用標(biāo)記控制結(jié)構(gòu)。例如，考慮以下代碼片段，其中標(biāo)簽 10 必須與相應(yīng)的 CONTINUE 匹配。

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

Fortran 90 消除了這種做法的需要，并通過允許使用 end do 來結(jié)束 do 循環(huán)來提高代碼可讀性。此更改顯示在上面提供的 dlacpy 的自由形式版本中。

假定大小的數(shù)組

為了能夠靈活地處理不同大小的數(shù)組，LAPACK 例程通常對(duì)具有假定大小的數(shù)組進(jìn)行操作。在上面的 dlacpy 例程中，輸入矩陣 A 被聲明為具有根據(jù)表達(dá)式 A(LDA, *) 的假定大小的二維數(shù)組。該表達(dá)式聲明 A 具有 LDA 行數(shù)，并使用 * 作為占位符來指示第二維的大小由調(diào)用程序確定。

使用假定大小數(shù)組的一個(gè)后果是編譯器無法對(duì)未指定的維度執(zhí)行邊界檢查。因此，當(dāng)前的最佳實(shí)踐是使用顯式接口和假定形狀數(shù)組（例如，A(LDA,:)）以防止越界內(nèi)存訪問。也就是說，當(dāng)需要將子矩陣傳遞給其他函數(shù)時(shí)，使用假定形狀數(shù)組可能會(huì)出現(xiàn)問題，因?yàn)檫@樣做需要切片，這通常會(huì)導(dǎo)致編譯器創(chuàng)建數(shù)組數(shù)據(jù)的內(nèi)部副本。

遷移到 Fortran 95

不用說，我花了一段時(shí)間才適應(yīng) LAPACK 慣例并采用 LAPACK 思維方式。然而，作為一個(gè)純粹主義者，如果我無論如何都要移植例程，我至少想將我設(shè)法移植的那些例程帶入更現(xiàn)代的時(shí)代，以期提高代碼的可讀性和未來的維護(hù)。因此，在與 stdlib 維護(hù)者討論后，我決定將例程遷移到 Fortran 95，雖然它不是最新、最好的 Fortran 版本，但似乎在保持原始實(shí)現(xiàn)的外觀和感覺之間取得了適當(dāng)?shù)钠胶猓_保（足夠好）向后兼容性，并利用更新的語法功能。

測(cè)試覆蓋率

采用自下而上的方法來添加 LAPACK 支持的問題之一是，LAPACK 中通常不存在針對(duì)較低級(jí)別實(shí)用程序例程的顯式單元測(cè)試。 LAPACK 的測(cè)試套件主要采用分層測(cè)試原理，其中假設(shè)測(cè)試較高級(jí)別的例程以確保其依賴的較低級(jí)別例程作為整個(gè)工作流程的一部分正常運(yùn)行。雖然有人可能會(huì)說，針對(duì)較低級(jí)別的例程，重點(diǎn)關(guān)注集成測(cè)試而不是單元測(cè)試是合理的，因?yàn)闉槊總€(gè)例程添加測(cè)試可能會(huì)增加 LAPACK 測(cè)試框架的維護(hù)負(fù)擔(dān)和復(fù)雜性，但這意味著我們不能輕易依賴先前的測(cè)試單元測(cè)試的藝術(shù)，并且必須自己為每個(gè)較低級(jí)別的例程提供全面的獨(dú)立單元測(cè)試。

文檔

與測(cè)試覆蓋率類似，在 LAPACK 本身之外，找到展示較低級(jí)別例程使用的真實(shí)記錄示例具有挑戰(zhàn)性。雖然 LAPACK 例程之前始終有一個(gè)文檔注釋，提供輸入?yún)?shù)和可能的返回值的描述，但如果沒有代碼示例，則可視化和理解預(yù)期的輸入和輸出值可能具有挑戰(zhàn)性，特別是在處理專門的矩陣時(shí)。雖然缺少單元測(cè)試和記錄示例都不是世界末日，但這意味著向 stdlib 添加 LAPACK 支持將比我預(yù)期的更加困難。編寫基準(zhǔn)、測(cè)試、示例和文檔只會(huì)需要更多的時(shí)間和精力，這可能會(huì)限制我在實(shí)習(xí)期間可以實(shí)施的例程數(shù)量。

內(nèi)存布局

在線性存儲(chǔ)器中存儲(chǔ)矩陣元素時(shí)，有兩種選擇：連續(xù)存儲(chǔ)列或連續(xù)存儲(chǔ)行（見圖 2）。前一種內(nèi)存布局稱為列優(yōu)先順序，后者稱為行優(yōu)先順序。

圖 2：演示以 (a) 列優(yōu)先（Fortran 樣式）或 (b) 行優(yōu)先（C 樣式）順序?qū)⒕仃囋卮鎯?chǔ)在線性存儲(chǔ)器中的示意圖。選擇使用哪種布局很大程度上是一個(gè)約定問題。

選擇使用哪種布局很大程度上是一個(gè)約定問題。例如，F(xiàn)ortran 以列優(yōu)先順序存儲(chǔ)元素，C 以行優(yōu)先順序存儲(chǔ)元素。更高級(jí)別的庫，例如 NumPy 和 stdlib，支持列優(yōu)先順序和行優(yōu)先順序，允許您在數(shù)組創(chuàng)建期間配置多維數(shù)組的布局。

pip install numpy

雖然兩種內(nèi)存布局本質(zhì)上都不比另一種更好，但根據(jù)底層存儲(chǔ)模型的約定安排數(shù)據(jù)以確保順序訪問對(duì)于確保最佳性能至關(guān)重要。現(xiàn)代 CPU 能夠比非順序數(shù)據(jù)更有效地處理順序數(shù)據(jù)，這主要?dú)w功于 CPU 緩存，而 CPU 緩存又利用了引用的空間局部性。

為了演示順序與非順序元素訪問的性能影響，請(qǐng)考慮以下函數(shù)，該函數(shù)將所有元素從 MxN 矩陣 A 復(fù)制到另一個(gè) MxN 矩陣 B，并且假設(shè)矩陣元素存儲(chǔ)在列優(yōu)先中訂購。

# Import all of NumPy:
import numpy as np

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

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

設(shè) A 和 B 為以下 3x2 矩陣：

當(dāng) A 和 B 都按列優(yōu)先順序存儲(chǔ)時(shí)，我們可以按如下方式調(diào)用復(fù)制例程：

pip install numpy

但是，如果 A 和 B 都按行優(yōu)先順序存儲(chǔ)，則調(diào)用簽名將更改為

# Import all of NumPy:
import numpy as np

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

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

請(qǐng)注意，在后一種情況下，我們無法在最內(nèi)層循環(huán)中按順序訪問元素，因?yàn)?da0 是 2，da1 是 -5，對(duì)于 db0 和 db1 也是如此。相反，數(shù)組索引“指針”在返回到線性內(nèi)存中較早的元素之前會(huì)重復(fù)向前跳過，其中 ia = {0, 2, 4, 1, 3, 5} 和 ib 相同。在圖 3 中，我們展示了非順序訪問對(duì)性能的影響。

圖 3：當(dāng) copy 假定根據(jù)列優(yōu)先順序進(jìn)行順序元素訪問時(shí)，向 copy 提供方列優(yōu)先矩陣與行優(yōu)先矩陣時(shí)的性能比較。 x 軸列舉增加的矩陣大?。丛?cái)?shù)量）。所有比率均相對(duì)于相應(yīng)矩陣大小的列主要結(jié)果進(jìn)行歸一化。

從圖中，我們可以觀察到列優(yōu)先和行優(yōu)先的性能大致相當(dāng)，直到我們對(duì)具有超過 1e5 個(gè)元素 (M = N = ~316) 的方陣進(jìn)行操作。對(duì)于 1e6 個(gè)元素 (M = N = ~1000)，提供行優(yōu)先矩陣進(jìn)行復(fù)制會(huì)導(dǎo)致性能下降 25% 以上。對(duì)于 1e7 元素 (M = N = ~3160)，我們觀察到性能下降超過 85%。顯著的性能影響可能歸因于在具有大行大小的行主矩陣上操作時(shí)引用局部性的降低。

鑒于 LAPACK 是用 Fortran 編寫的，它假定列優(yōu)先訪問順序并相應(yīng)地實(shí)現(xiàn)其算法。這給 stdlib 等庫帶來了問題，它們不僅支持行主序，而且將其設(shè)為默認(rèn)內(nèi)存布局。如果我們只是將 LAPACK 的 Fortran 實(shí)現(xiàn)移植到 JavaScript，那么提供行主矩陣的用戶將會(huì)因非順序訪問而遭受不利的性能影響。

為了減輕不利的性能影響，我們借鑒了 BLIS（一個(gè)類似 BLAS 的庫，支持 BLAS 例程中的行主內(nèi)存布局和列主內(nèi)存布局）的想法，并決定在將例程從 Fortran 移植到 JavaScript 時(shí)創(chuàng)建修改后的 LAPACK 實(shí)現(xiàn)， C 通過每個(gè)維度的單獨(dú)步幅參數(shù)顯式適應(yīng)列主內(nèi)存布局和行主內(nèi)存布局。對(duì)于某些實(shí)現(xiàn)，例如與上面定義的復(fù)制函數(shù)類似的 dlacpy，合并單獨(dú)且獨(dú)立的步幅是簡(jiǎn)單的，通常涉及步幅技巧和循環(huán)交換，但是，對(duì)于其他實(shí)現(xiàn)，由于以下原因，修改結(jié)果并不那么簡(jiǎn)單。專門的矩陣處理、不同的訪問模式和組合參數(shù)化。

ndarrays

LAPACK 例程主要對(duì)存儲(chǔ)在線性存儲(chǔ)器中的矩陣進(jìn)行操作，并且根據(jù)指定的維度和前導(dǎo)（即第一個(gè)）維度的步長來訪問其元素。維度分別指定每行和每列中的元素?cái)?shù)量。步長指定必須跳過線性內(nèi)存中的多少個(gè)元素才能訪問行的下一個(gè)元素。 LAPACK 假設(shè)屬于同一列的元素總是連續(xù)的（即在線性存儲(chǔ)器中相鄰）。圖 4 提供了 LAPACK 約定的直觀表示（特別是原理圖 (a) 和 (b)）。

圖 4：示意圖說明了 LAPACK 跨步數(shù)組約定泛化為非連續(xù)跨步數(shù)組。 a) 以列優(yōu)先順序存儲(chǔ)的 5×5 連續(xù)矩陣。 b) 以列優(yōu)先順序存儲(chǔ)的 3×3 非連續(xù)子矩陣。通過提供指向第一個(gè)索引元素的指針并指定前導(dǎo)（即第一個(gè)）維度的步長，可以在 LAPACK 中對(duì)子矩陣進(jìn)行操作。在這種情況下，盡管每列只有三個(gè)元素，但由于當(dāng)存儲(chǔ)為較大矩陣的一部分時(shí)，線性存儲(chǔ)器中的子矩陣元素不連續(xù)，因此前導(dǎo)維度的步長為 5。在 LAPACK 中，尾隨（即第二個(gè)）維度的步長始終假定為統(tǒng)一。 c) 以列主順序存儲(chǔ)的 3×3 非連續(xù)子矩陣，具有非單位步長，并將 LAPACK 步長約定推廣到前導(dǎo)維度和尾隨維度。這種泛化支撐了 stdlib 的多維數(shù)組（也稱為“ndarrays”）。

NumPy 和 stdlib 等庫概括了 LAPACK 的跨步數(shù)組約定以支持

1. 最后一個(gè)維度中的非單位步長（見圖 4 (c)）。 LAPACK 假設(shè)矩陣的最后一個(gè)維度始終具有單位步長（即列中的元素連續(xù)存儲(chǔ)在線性內(nèi)存中）。
2. 任何維度的負(fù)步幅。 LAPACK 要求前導(dǎo)矩陣維度的步長為正。
3. 具有兩個(gè)以上維度的多維數(shù)組。 LAPACK 僅明確支持跨步向量和（子）矩陣。

對(duì)最后一個(gè)維度中非單位步長的支持確保支持 O(1) 創(chuàng)建線性內(nèi)存的非連續(xù)視圖，而無需顯式數(shù)據(jù)移動(dòng)。這些視圖通常稱為“切片”。作為示例，請(qǐng)考慮以下代碼片段，該代碼片段使用 stdlib 提供的 API 創(chuàng)建此類視圖。

pip install numpy

如果最后一個(gè)維度不支持非單位步幅，則不可能從表達(dá)式 x['::2,::2'] 返回視圖，因?yàn)樾枰獙⑺x元素復(fù)制到新的線性內(nèi)存緩沖區(qū)以確保連續(xù)性。

圖 5：示意圖，說明如何使用步幅操作來創(chuàng)建存儲(chǔ)在線性存儲(chǔ)器中的矩陣元素的翻轉(zhuǎn)和旋轉(zhuǎn)視圖。對(duì)于所有子原理圖，步幅列為 [trailing_dimension,leading_dimension]。每個(gè)原理圖隱含一個(gè)“偏移量”，它指示第一個(gè)索引元素的索引，這樣，對(duì)于矩陣 A，元素 A_ij 是根據(jù) i?strides[1] j?strides[0] 偏移量解析。 a) 給定一個(gè)以列優(yōu)先順序存儲(chǔ)的 3×3 矩陣，可以操縱前導(dǎo)維度和尾隨維度的步幅來創(chuàng)建視圖，在該視圖中以相反的順序訪問沿一個(gè)或多個(gè)軸的矩陣元素。 b) 使用類似的步幅操作，可以創(chuàng)建矩陣元素相對(duì)于它們?cè)诰€性內(nèi)存中的排列的旋轉(zhuǎn)視圖。

對(duì)負(fù)步幅的支持使得元素能夠沿一個(gè)或多個(gè)維度進(jìn)行 O(1) 反轉(zhuǎn)和旋轉(zhuǎn)（參見圖 5）。例如，要從上到下和從左到右翻轉(zhuǎn)矩陣，只需取消步幅即可。在前面的代碼片段的基礎(chǔ)上，以下代碼片段演示了圍繞一個(gè)或多個(gè)軸反轉(zhuǎn)元素。

pip install numpy

負(fù)步幅的討論中隱含的是需要一個(gè)“偏移”參數(shù)，該參數(shù)指示線性內(nèi)存中第一個(gè)索引元素的索引。對(duì)于跨步多維數(shù)組 A 和跨步列表 s，元素 A_ij???n 對(duì)應(yīng)的索引可以根據(jù)方程求解

其中 N 是數(shù)組維度數(shù)，s_k 對(duì)應(yīng)于第 k 個(gè)步幅。

在支持負(fù)步幅的 BLAS 和 LAPACK 例程中（僅在對(duì)步幅向量進(jìn)行操作時(shí)才支持這種功能（例如，參見上面的 daxpy）），索引偏移量是使用類似于以下代碼片段的邏輯來計(jì)算的：

pip install numpy

其中 M 是向量元素的數(shù)量。這隱含地假設(shè)提供的數(shù)據(jù)指針指向向量的線性存儲(chǔ)器的開頭。在支持指針的語言中，例如C，為了在線性存儲(chǔ)器的不同區(qū)域上進(jìn)行操作，通常在函數(shù)調(diào)用之前使用指針?biāo)阈g(shù)來調(diào)整指針，這相對(duì)便宜且簡(jiǎn)單，至少對(duì)于一維情況是這樣。

例如，返回到上面定義的 c_daxpy，我們可以使用指針?biāo)阈g(shù)來限制對(duì)線性內(nèi)存中從輸入和輸出數(shù)組的第十一個(gè)和第十六個(gè)元素（注意：從零開始索引）開始的五個(gè)元素的元素訪問，分別如下面的代碼片段所示。

# 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 中，不支持二進(jìn)制緩沖區(qū)的顯式指針?biāo)阈g(shù)，必須顯式實(shí)例化具有所需字節(jié)偏移量的新類型數(shù)組對(duì)象。在下面的代碼片段中，為了獲得與上面的 C 示例相同的結(jié)果，我們必須解析一個(gè)類型化數(shù)組構(gòu)造函數(shù)，計(jì)算一個(gè)新的字節(jié)偏移量，計(jì)算一個(gè)新的類型化數(shù)組長度，并創(chuàng)建一個(gè)新的類型化數(shù)組實(shí)例。

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

對(duì)于大型數(shù)組，與訪問和操作單個(gè)數(shù)組元素所花費(fèi)的時(shí)間相比，類型化數(shù)組實(shí)例化的成本可以忽略不計(jì)；然而，對(duì)于較小的數(shù)組大小，對(duì)象實(shí)例化會(huì)顯著影響性能。

因此，為了避免對(duì)對(duì)象實(shí)例化性能產(chǎn)生不利影響，stdlib 將 ndarray 的數(shù)據(jù)緩沖區(qū)與 ndarray 視圖開頭對(duì)應(yīng)的緩沖區(qū)元素的位置解耦。這允許切片表達(dá)式 x[2:,3:] 和 x[3:,1:] 返回新的 ndarray 視圖而無需需要實(shí)例化新的緩沖區(qū)實(shí)例，如以下代碼片段所示。

// 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,:,:']);

由于將數(shù)據(jù)緩沖區(qū)與 ndarray 視圖的開頭解耦，我們同樣試圖避免在使用 ndarray 數(shù)據(jù)調(diào)用 LAPACK 例程時(shí)實(shí)例化新的類型化數(shù)組實(shí)例。這意味著創(chuàng)建修改后的 LAPACK API 簽名，支持所有跨步向量和矩陣的顯式偏移參數(shù)。

為了簡(jiǎn)單起見，讓我們回到上面定義的 daxpy 的 JavaScript 實(shí)現(xiàn)。

pip install numpy

如以下代碼片段所示，我們可以修改上述簽名和實(shí)現(xiàn)，以便將解析第一個(gè)索引元素的責(zé)任轉(zhuǎn)移給 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,:,:]

對(duì)于 ndarray，解析發(fā)生在 ndarray 實(shí)例化期間，使得使用 ndarray 數(shù)據(jù)調(diào)用 daxpy_ndarray 可以直接傳遞關(guān)聯(lián)的 ndarray 元數(shù)據(jù)。下面的代碼片段對(duì)此進(jìn)行了演示。

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

與 BLIS 類似，我們看到了傳??統(tǒng) LAPACK API 簽名（例如，為了向后兼容）和修改后的 API 簽名（例如，為了最大限度地減少不利的性能影響）的價(jià)值，因此，我們制定了一項(xiàng)計(jì)劃，提供傳統(tǒng)和修改了每個(gè) LAPACK 例程的 API。為了最大限度地減少代碼重復(fù)，我們的目標(biāo)是實(shí)現(xiàn)一個(gè)通用的較低級(jí)別的“基本”實(shí)現(xiàn)，然后可以由較高級(jí)別的 API 進(jìn)行包裝。雖然上面顯示的 BLAS 例程 daxpy 的更改可能看起來相對(duì)簡(jiǎn)單，但傳統(tǒng) LAPACK 例程及其預(yù)期行為到通用實(shí)現(xiàn)的轉(zhuǎn)換通常要復(fù)雜得多。

德拉斯普

挑戰(zhàn)夠了！最終產(chǎn)品是什么樣的？！

讓我們繞一圈，回到 dlaswp，這是一個(gè) LAPACK 例程，用于根據(jù)主元索引列表在輸入矩陣上執(zhí)行一系列行交換。以下代碼片段顯示了參考 LAPACK Fortran 實(shí)現(xiàn)。

// 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 實(shí)現(xiàn)進(jìn)行交互，LAPACK 提供了一個(gè)名為 LAPACKE 的兩級(jí) C 接口，它包裝了 Fortran 實(shí)現(xiàn)，并為行主和列主輸入和輸出矩陣提供了便利。 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;
}

當(dāng)使用列主矩陣 a 調(diào)用時(shí)，包裝器 LAPACKE_dlaswp_work 只是將提供的參數(shù)傳遞給 Fortran 實(shí)現(xiàn)。但是，當(dāng)使用行主矩陣 a 調(diào)用時(shí)，包裝器必須分配內(nèi)存、顯式轉(zhuǎn)置并將 a 復(fù)制到臨時(shí)矩陣 a_t、重新計(jì)算前導(dǎo)維度的步幅、使用 a_t 調(diào)用 dlaswp、轉(zhuǎn)置并復(fù)制存儲(chǔ)在 a_t 中的結(jié)果到a，最后釋放分配的內(nèi)存。這是相當(dāng)大的工作量，并且在大多數(shù) LAPACK 例程中很常見。

以下代碼片段顯示了移植到 JavaScript 的參考 LAPACK 實(shí)現(xiàn)，支持前導(dǎo)和尾隨維度步幅、索引偏移以及包含主元索引的步幅向量。

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

const mod = new WebAssembly.Module(binary);

為了提供與傳統(tǒng) LAPACK 具有一致行為的 API，我包裝了上述實(shí)現(xiàn)并將輸入?yún)?shù)調(diào)整為“基本”實(shí)現(xiàn)，如以下代碼片段所示。

pip install numpy

我隨后編寫了一個(gè)單獨(dú)但類似的包裝器，它提供了更直接到 stdlib 多維數(shù)組的 API 映射，并在應(yīng)用樞軸的方向?yàn)樨?fù)時(shí)執(zhí)行一些特殊處理，如以下代碼片段所示。

# Import all of NumPy:
import numpy as np

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

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

需要注意的幾點(diǎn)：

1. 與傳統(tǒng)的 LAPACKE API 相比，在 dlaswp_ndarray 和基礎(chǔ) API 中，matrix_layout（順序）參數(shù)不是必需的，因?yàn)榭梢詮奶峁┑牟椒茢囗樞颉?/li>
2. 與傳統(tǒng)的LAPACKE API相比，當(dāng)輸入矩陣為行主矩陣時(shí)，我們不需要將數(shù)據(jù)復(fù)制到臨時(shí)工作區(qū)數(shù)組，從而減少了不必要的內(nèi)存分配。
3. 與直接與 BLAS 和 LAPACK 交互的 NumPy 和 SciPy 等庫相比，在 stdlib 中調(diào)用 LAPACK 例程時(shí)，我們不需要之前將不連續(xù)的多維數(shù)據(jù)復(fù)制到臨時(shí)工作區(qū)數(shù)組或從臨時(shí)工作區(qū)數(shù)組復(fù)制非連續(xù)的多維數(shù)據(jù)和調(diào)用后，分別。除了與硬件優(yōu)化的 BLAS 和 LAPACK 接口時(shí)之外，所采用的方法有助于最大限度地減少數(shù)據(jù)移動(dòng)并確保資源受限的瀏覽器應(yīng)用程序的性能。

對(duì)于希望利用硬件優(yōu)化庫（例如 OpenBLAS）的服務(wù)器端應(yīng)用程序，我們提供單獨(dú)的包裝器，將通用簽名參數(shù)適應(yīng)其優(yōu)化的 API 等效項(xiàng)。在這種情況下，至少對(duì)于足夠大的數(shù)組來說，創(chuàng)建臨時(shí)副本是值得的。

當(dāng)前狀態(tài)和后續(xù)步驟

盡管面臨挑戰(zhàn)、不可預(yù)見的挫折和多次設(shè)計(jì)迭代，我很高興地報(bào)告，除了上面的 dlaswp 之外，我還能夠打開 35 個(gè) PR，添加對(duì)各種 LAPACK 例程和相關(guān)實(shí)用程序的支持。顯然不完全是 1,700 個(gè)例程，但這是一個(gè)好的開始！ :)

盡管如此，未來是光明的，我們對(duì)這項(xiàng)工作感到非常興奮。仍有很大的改進(jìn)空間和額外的研究和開發(fā)。我們尤其渴望

1. 探索工具和自動(dòng)化。
2. 解決了解析跨多個(gè) stdlib 包的 Fortran 依賴項(xiàng)的源文件時(shí)的構(gòu)建問題。
3. 為 stdlib 的現(xiàn)有 LAPACK 包推出 C 和 Fortran 實(shí)現(xiàn)以及本機(jī)綁定。
4. 繼續(xù)發(fā)展 stdlib 的模塊化 LAPACK 例程庫。
5. 確定性能優(yōu)化的其他領(lǐng)域。

雖然我的 Quansight Labs 實(shí)習(xí)已經(jīng)結(jié)束，但我的計(jì)劃是繼續(xù)添加軟件包并推動(dòng)這項(xiàng)工作。鑒于 LAPACK 的巨大潛力和根本重要性，我們很高興看到將 LAPACK 引入網(wǎng)絡(luò)的這一舉措繼續(xù)發(fā)展，因此，如果您有興趣幫助推動(dòng)這一進(jìn)程，請(qǐng)隨時(shí)與我們聯(lián)系！如果您有興趣贊助開發(fā)，Quansight 的人員將非常樂意與您聊天。

在此，我要感謝 Quansight 提供的這次實(shí)習(xí)機(jī)會(huì)。我感到非常幸運(yùn)能夠?qū)W到這么多東西。在 Quansight 實(shí)習(xí)是我長期以來的一個(gè)夢(mèng)想，我很高興能夠?qū)崿F(xiàn)這個(gè)夢(mèng)想。我要特別感謝 Athan Reines 和 Melissa Mendon?a，他們是一位出色的導(dǎo)師，也是一位出色的人！感謝所有 stdlib 核心開發(fā)人員以及 Quansight 的其他所有人，感謝你們一路上大大小小的幫助了我。

干杯！

---

stdlib 是一個(gè)開源軟件項(xiàng)目，致力于提供一整套強(qiáng)大、高性能的庫來加速您的項(xiàng)目開發(fā)，并讓您安心，因?yàn)槟滥蕾嚨氖蔷闹谱鞯母哔|(zhì)量軟件。

如果您喜歡這篇文章，請(qǐng)給我們一顆星？ GitHub 上并考慮為該項(xiàng)目提供經(jīng)濟(jì)支持。您的貢獻(xiàn)和持續(xù)的支持有助于確保項(xiàng)目的長期成功，我們深表感謝！

以上是網(wǎng)絡(luò)瀏覽器中的 LAPACK的詳細(xì)內(nèi)容。更多信息請(qǐng)關(guān)注PHP中文網(wǎng)其他相關(guān)文章！