這篇文章最初發(fā)佈在 Quansight Labs 部落格上，經(jīng) Quansight 許可修改並重新發(fā)佈於此處。

Web 應(yīng)用程式正在迅速崛起，成為高效能科學(xué)運(yùn)算和支援 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)路獨(dú)特限制量身定制的全面、高品質(zhì)的 LAPACK 實(shí)作尚未實(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)頁瀏覽器成為數(shù)值運(yùn)算和機(jī)器學(xué)習(xí)的一流環(huán)境提供關(guān)鍵的建構(gòu)模組，並預(yù)示著更強(qiáng)大的人工智慧網(wǎng)路應(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ù)組資料結(jié)構(gòu)以及數(shù)學(xué)、統(tǒng)計(jì)和線性代數(shù)的相關(guān)例程，但使用 JavaScript 而不是 Python 作為其主要腳本語言。因此，stdlib 專注於 Web 生態(tài)系統(tǒng)及其應(yīng)用程式開發(fā)範(fàn)例。這種關(guān)注需要一些有趣的設(shè)計(jì)和專案架構(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)行陣列操作。如果您正在部署一個(gè)僅需要 NumPy 的 ndarray 物件及其一些操作例程的應(yīng)用程序，那麼安裝和捆綁所有 NumPy 意味著包含大量「死代碼」。用 Web 開發(fā)的術(shù)語來說，我們會(huì)說 NumPy 不是「tree shakeable」。對(duì)於正常的 NumPy 安裝，這表示至少需要 30MB 的磁碟空間，對(duì)於排除所有偵錯(cuò)語句的自訂建置至少需要 15MB 的磁碟空間。對(duì)於 SciPy，這些數(shù)字可以分別增加到 130MB 和 50MB。不用說，在 Web 應(yīng)用程式中僅提供幾個(gè)功能的 15MB 庫是不可能的，特別是對(duì)於需要將 Web 應(yīng)用程式部署到網(wǎng)路連接較差或記憶體受限的裝置的開發(fā)人員而言。

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

舉個(gè)例子，假設(shè)您正在使用兩個(gè)矩陣堆疊（即三維立方體的二維切片），並且您想要選擇每隔一個(gè)切片並執(zhí)行常見的BLAS 運(yùn)算y = a * x，其中x 和y是ndarray，a 是標(biāo)量常數(shù)。要使用 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ū)０赴惭b為整體庫之外，您還可以將各種功能單元安裝為單獨(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)頁瀏覽器中的 LAPACK 有什麼關(guān)係呢？嗯，去年夏天我們的目標(biāo)之一是應(yīng)用 stdlib 精神——小型、範(fàn)圍狹窄的軟體包，只做一件事，並做好一件事——將 LAPACK 引入網(wǎng)路。

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

不幸的是，這種方法存在幾個(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，並且產(chǎn)生的程式碼需要大量修補(bǔ)。開發(fā)基於 LLVM 的 Fortran 編譯器的工作正在進(jìn)行中，但仍存在差距和複雜的工具鏈。
2. 如同上面關(guān)於 Web 應(yīng)用程式中的整體函式庫的討論中所提到的，LAPACK 的龐大是問題的一部分。即使編譯問題已解決，在僅需要使用一兩個(gè) LAPACK 例程的 Web 應(yīng)用程式中包含包含所有 LAPACK 的單一 WebAssembly 二進(jìn)位檔案也意味著大量死程式碼，導(dǎo)致載入時(shí)間變慢並增加記憶體消耗。
3. 雖然人們可以嘗試將單一 LAPACK 例程編譯為獨(dú)立的 WebAssembly 二進(jìn)位文件，但這樣做可能會(huì)導(dǎo)致二進(jìn)位膨脹，因?yàn)槎鄠€(gè)獨(dú)立二進(jìn)位檔案可能包含來自公共依賴項(xiàng)的重複程式碼。為了減輕二進(jìn)位膨脹，可以嘗試執(zhí)行模組拆分。在這種情況下，首先將公共依賴項(xiàng)分解為包含共享程式碼的獨(dú)立二進(jìn)位文件，然後為各個(gè) API 產(chǎn)生單獨(dú)的二進(jìn)位。雖然在某些情況下適用，但這種方法很快就會(huì)變得笨拙，因?yàn)檫@種方法需要在加載時(shí)通過將一個(gè)或多個(gè)模組的導(dǎo)出與一個(gè)或多個(gè)其他模組的導(dǎo)入縫合在一起來連結(jié)各個(gè)WebAssembly 模組。這不僅很乏味，而且這種方法還會(huì)帶來效能損失，因?yàn)楫?dāng) WebAssembly 例程呼叫匯入的匯出時(shí)，它們現(xiàn)在必須跨入 JavaScript，而不是保留在 WebAssembly 中。聽起來很複雜？是的！
4. 除了專門對(duì)標(biāo)量輸入?yún)?shù)進(jìn)行操作的WebAssembly 模組（例如，計(jì)算單一數(shù)字的正弦值）之外，每個(gè)WebAssembly 模組實(shí)例都必須與WebAssembly 記憶體關(guān)聯(lián)，該記憶體以64KiB 的固定增量分配（即「頁面」） 」）。重要的是，截至這篇博文，WebAssembly 內(nèi)存只能增長而永遠(yuǎn)不會(huì)收縮。由於當(dāng)前沒有向主機(jī)釋放內(nèi)存的機(jī)制，因此WebAssembly 應(yīng)用程式的記憶體佔(zhàn)用只會(huì)增加。
最後，雖然 WebAssembly 功能強(qiáng)大，但它需要更陡峭的學(xué)習(xí)曲線和一組更複雜且經(jīng)常快速發(fā)展的工具鏈。在最終用戶應(yīng)用程式中，JavaScript（一種 Web 原生動(dòng)態(tài)編譯程式語言）和 WebAssembly 之間的介面進(jìn)一步增加了複雜性，尤其是在必須執(zhí)行手動(dòng)記憶體管理時(shí)。

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

pip install numpy

編譯到 WebAssembly 並將 WebAssembly 二進(jìn)位檔案載入到我們的 Web 應(yīng)用程式後，我們需要執(zhí)行一系列步驟，然後才能從 JavaScript 呼叫 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,:,:]

接下來，我們需要定義模組記憶體並建立一個(gè)新的 WebAssembly 模組實(shí)例。

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

建立模組實(shí)例後，我們現(xiàn)在可以呼叫導(dǎo)出的 BLAS 例程。但是，如果資料是在模組記憶體之外定義的，我們首先需要將該資料複製到記憶體實(shí)例，並且始終以小端位元組順序執(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)在資料已寫入模組內(nèi)存，我們可以呼叫 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 記憶體「指標(biāo)」（即位元組偏移）的下游函式庫（例如D3）以進(jìn)行視覺化或進(jìn)一步分析，我們需要從模組複製資料記憶體回原始輸出數(shù)組。

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

const mod = new WebAssembly.Module(binary);

光計(jì)算 y = a*x y 就需要大量工作。相反，與純 JavaScript 實(shí)作相比，它可能類似於以下程式碼片段。

// 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)，我們可以使用外部定義的資料直接呼叫 daxpy，而無需上面 WebAssembly 範(fàn)例中所需的資料移動(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ù)期的，考慮到所需的資料移動(dòng)，還會(huì)對(duì)效能產(chǎn)生負(fù)面影響，如下圖所示。

圖 1：對(duì)於增加數(shù)組長度（x 軸）的 BLAS 例程 daxpy，stdlib 的 C、JavaScript 和 WebAssembly (Wasm) 實(shí)現(xiàn)的效能比較。在Wasm（複製）基準(zhǔn)測(cè)試中，輸入和輸出資料被複製到Wasm記憶體中或從Wasm記憶體複製，導(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è)試中，輸入和輸出資料直接在WebAssembly 模組記憶體中分配和操作，並且在Wasm（複製）基準(zhǔn)測(cè)試中，輸入和輸出資料複製到WebAssembly 模組記憶體中或從WebAssembly 模組記憶體中複製，如上所述。從圖表中，我們可以觀察到以下幾點(diǎn)：

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

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

徹底的模組化

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

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

我的任務(wù)是為這樣的努力奠定基礎(chǔ)，解決問題並找到差距，並希望讓我們更接近網(wǎng)路上的高效能線性代數(shù)。但是徹底的模組化是什麼樣的呢？這一切都始於基本的功能單元，套件。

stdlib 中的每個(gè)套件都是獨(dú)立的，包含共同本地化的測(cè)試、基準(zhǔn)測(cè)試、範(fàn)例、文件、建立文件和關(guān)聯(lián)的元資料（包括任何依賴項(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í)作（即參考 LAPACK）的效能。
• docs：包含輔助文件的資料夾，其中包括 REPL 幫助文字和定義類型化 API 簽署的 TypeScript 聲明。
• examples：包含可執(zhí)行示範(fàn)程式碼的資料夾，除了充當(dāng)文件之外，還可以幫助開發(fā)人員健全性檢查實(shí)作行為。
• include：包含 C 頭檔案的資料夾。
• lib：包含 JavaScript 來源實(shí)作的資料夾，其中 index.js 作為套件入口點(diǎn)，其他 *.js 檔案定義內(nèi)部實(shí)作模組。
• src：包含 C 和 Fortran 來源實(shí)作的資料夾。每個(gè)模組化 LAPACK 套件都應(yīng)包含稍作修改的 Fortran 參考實(shí)作（F77 到自由格式的 Fortran）。 C 檔案包括遵循Fortran 參考實(shí)作的純C 實(shí)作、用於呼叫Fortran 參考實(shí)作的包裝器、用於在伺服器端應(yīng)用程式中呼叫硬體最佳化函式庫（例如OpenBLAS）的包裝器，以及用於呼叫已編譯的本機(jī)器綁定來自Node.js 中的JavaScript 或相容的伺服器端JavaScript 執(zhí)行時(shí)期的C。
• test：包含單元測(cè)試的資料夾，用於測(cè)試 JavaScript 和本機(jī)實(shí)作中的預(yù)期行為。本機(jī)實(shí)作的測(cè)試是用 JavaScript 編寫的，並利用本機(jī)綁定實(shí)作 JavaScript 和 C/Fortran 之間的互通性。
• binding.gyp/include.gypi：用於編譯 Node.js 原生插件的建置文件，它提供了 JavaScript 和原生程式碼之間的橋樑。
• manifest.json：stdlib內(nèi)部C和Fortran編譯原始檔包管理的設(shè)定檔。
• package.json：包含包元資料的文件，包括外部包依賴項(xiàng)的枚舉以及用於基於瀏覽器的 Web 應(yīng)用程式的純 JavaScript 實(shí)現(xiàn)的路徑。
• README.md：包含套件的主要文件的文件，其中包括 API 簽章以及 JavaScript 和 C 介面的範(fàn)例。

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

多階段方法

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

為了確定首先定位哪些 LAPACK 例程，我解析了 LAPACK 的 Fortran 原始碼以產(chǎn)生呼叫圖。這使我能夠推斷出每個(gè) LAPACK 例程的依賴關(guān)係樹。有了這張圖，我就進(jìn)行了拓?fù)渑判?，從而幫助我識(shí)別沒有依賴性的例程，並且這些例程將不可避免地成為其他例程的構(gòu)建塊。雖然我選擇一個(gè)特定的高級(jí)例程並向後工作的深度優(yōu)先方法將使我能夠?qū)崿F(xiàn)特定的功能，但這種方法可能會(huì)導(dǎo)致我在嘗試實(shí)現(xiàn)日益複雜的例程時(shí)陷入困境。透過專注於圖表的“葉子”，我可以優(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)代互動(dòng)式 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í)作中。這些約定中最明顯的一種是格式。

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

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

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

pip install numpy

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

# 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 必須與對(duì)應(yīng)的 CONTINUE 相符。

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

Fortran 90 消除了這種做法的需要，並透過允許使用 end do 來結(jié)束 do 迴圈來提高程式碼可讀性。此變更顯示在上面提供的 dlacpy 的自由形式版本中。

假定大小的數(shù)組

為了能夠靈活地處理不同大小的數(shù)組，LAPACK 例程通常會(huì)對(duì)具有假定大小的數(shù)組進(jìn)行操作。在上面的 dlacpy 例程中，輸入矩陣 A 被宣告為具有根據(jù)表達(dá)式 A(LDA, *) 的假定大小的二維數(shù)組。此表達(dá)式宣告 A 具有 LDA 行數(shù)，並使用 * 作為佔(zhàn)位符來指示第二維的大小由呼叫程式?jīng)Q定。

使用假定大小數(shù)組的一個(gè)後果是編譯器無法對(duì)未指定的維度執(zhí)行邊界檢查。因此，目前的最佳實(shí)踐是使用明確介面和假定形狀數(shù)組（例如，A(LDA,:)）以防止越界記憶體存取。也就是說，當(dāng)需要將子矩陣傳遞給其他函數(shù)時(shí)，使用假定形狀數(shù)組可能會(huì)出現(xiàn)問題，因?yàn)檫@樣做需要切片，這通常會(huì)導(dǎo)致編譯器建立數(shù)組資料的內(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í)記錄範(fàn)例具有挑戰(zhàn)性。雖然LAPACK 例程之前始終有一個(gè)文件註釋，提供輸入?yún)?shù)和可能的返回值的描述，但如果沒有程式碼範(fàn)例，則可視化和理解預(yù)期的輸入和輸出值可能具有挑戰(zhàn)性，特別是在處理專門的矩陣時(shí)。雖然缺少單元測(cè)試和記錄範(fàn)例都不是世界末日，但這意味著向 stdlib 添加 LAPACK 支援將比我預(yù)期的更加困難。編寫基準(zhǔn)、測(cè)試、範(fàn)例和文件只會(huì)需要更多的時(shí)間和精力，這可能會(huì)限制我在實(shí)習(xí)期間可以實(shí)施的例程數(shù)量。

記憶體佈局

在線性記憶體中儲(chǔ)存矩陣元素時(shí)，有兩種選擇：連續(xù)儲(chǔ)存列或連續(xù)儲(chǔ)存行（請(qǐng)參閱圖 2）。前一種記憶體佈局稱為列優(yōu)先順序，後者稱為行優(yōu)先順序。

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

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

pip install numpy

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

為了示範(fàn)順序與非順序元素存取的效能影響，請(qǐng)考慮以下函數(shù)，該函數(shù)將所有元素從 MxN 矩陣 A 複製到另一個(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 矩陣：

0 00 >0000> >???

當(dāng) A 和 B 都按列優(yōu)先順序儲(chǔ)存時(shí)，我們可以如下呼叫複製例程：

pip install numpy

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

# 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ù)組索引「指標(biāo)」在返回線性記憶體中較早的元素之前會(huì)重複向前跳過，其中 ia = {0, 2, 4, 1, 3, 5} 和 ib 相同。在圖 3 中，我們展示了非順序存取對(duì)效能的影響。

圖 3：當(dāng) copy 假設(shè)依照列優(yōu)先順序進(jìn)行順序元素存取時(shí)，向 copy 提供方列優(yōu)先矩陣與行優(yōu)先矩陣時(shí)的效能比較。 x 軸列舉增加的矩陣大?。丛?cái)?shù)）。所有比率均相對(duì)於對(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)行複製會(huì)導(dǎo)致效能下降 25% 以上。對(duì)於 1e7 元素 (M = N = ~3160)，我們觀察到效能下降超過 85%。顯著的效能影響可能歸因於在具有大行大小的行主矩陣上操作時(shí)引用局部性的降低。

鑑於 LAPACK 是用 Fortran 編寫的，它假定列優(yōu)先存取順序並相應(yīng)地實(shí)現(xiàn)其演算法。這給 stdlib 等函式庫帶來了問題，它們不僅支援行主序，而且將其設(shè)為預(yù)設(shè)記憶體佈局。如果我們只是將 LAPACK 的 Fortran 實(shí)作移植到 JavaScript，那麼提供行主矩陣的使用者將會(huì)因非順序存取而遭受不利的效能影響。

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

ndarrays

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

圖 4：示意圖說明了 LAPACK 跨步數(shù)組約定泛化為非連續(xù)跨步數(shù)組。 a) 以列優(yōu)先順序儲(chǔ)存的 5×5 連續(xù)矩陣。 b) 以列優(yōu)先順序儲(chǔ)存的 3×3 非連續(xù)子矩陣。透過提供指向第一個(gè)索引元素的指標(biāo)並指定前導(dǎo)（即第一個(gè)）維度的步長，可以在 LAPACK 中對(duì)子矩陣進(jìn)行操作。在這種情況下，儘管每列只有三個(gè)元素，但由於當(dāng)儲(chǔ)存為較大矩陣的一部分時(shí)，線性記憶體中的子矩陣元素不連續(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ǔ)存在線性記憶體中）。
2. 任何維度的負(fù)步幅。 LAPACK 要求前導(dǎo)矩陣維度的步長為正。
3. 具有兩個(gè)以上維度的多維數(shù)組。 LAPACK 僅明確支援跨步向量和（子）矩陣。

對(duì)最後一個(gè)維度中非單位步長的支援確保支援 O(1) 建立線性記憶體的非連續(xù)視圖，而無需明確資料移動(dòng)。這些視圖通常稱為“切片”。作為範(fàn)例，請(qǐng)考慮以下程式碼片段，該程式碼片段使用 stdlib 提供的 API 建立此類視圖。

pip install numpy

如果最後一個(gè)維度不支援非單位步幅，則不可能從表達(dá)式x['::2,::2'] 返回視圖，因?yàn)樾枰獙⑺x元素複製到新的線性記憶體緩衝區(qū)以確保連續(xù)性。

圖 5：示意圖，說明如何使用步幅操作來建立儲(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)維度和尾隨維度的步幅來建立視圖，在該視圖中以相反的順序存取沿一個(gè)或多個(gè)軸的矩陣元素。 b) 使用類似的步幅操作，可以建立矩陣元素相對(duì)於它們?cè)诰€性記憶體中的排列的旋轉(zhuǎn)視圖。

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

pip install numpy

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

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

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

pip install numpy

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

例如，回到上面定義的c_daxpy，我們可以使用指標(biāo)算術(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)算術(shù)，必須明確實(shí)例化具有所需位元組偏移量的新類型數(shù)組物件。在下面的程式碼片段中，為了獲得與上面的C 範(fàn)例相同的結(jié)果，我們必須解析一個(gè)類型化數(shù)組構(gòu)造函數(shù)，計(jì)算一個(gè)新的位元組偏移量，計(jì)算一個(gè)新的類型化數(shù)組長度，並建立一個(gè)新的類型化數(shù)組實(shí)例。

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

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

因此，為了避免對(duì)物件實(shí)例化效能產(chǎn)生不利影響，stdlib 將 ndarray 的資料緩衝區(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,:,:']);

由於將資料緩衝區(qū)與 ndarray 視圖的開頭解耦，我們同樣試圖避免在使用 ndarray 資料呼叫 LAPACK 程式時(shí)實(shí)例化新的類型化數(shù)組實(shí)例。這意味著建立修改後的 LAPACK API 簽名，支援所有跨步向量和矩陣的明確偏移參數(shù)。

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

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 資料呼叫 daxpy_ndarray 可以直接傳遞關(guān)聯(lián)的 ndarray 元資料。下面的程式碼片段對(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。為了最大限度地減少程式碼重複，我們的目標(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)換通常要複雜得多。

德拉斯普

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

讓我們繞一圈，回到 dlaswp，這是一個(gè) LAPACK 例程，用於根據(jù)主元索引列表在輸入矩陣上執(zhí)行一系列行交換。以下程式碼片段顯示了參考 LAPACK Fortran 實(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,:,:']);

為了方便與 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 呼叫時(shí)，包裝器 LAPACKE_dlaswp_work 只是將提供的參數(shù)傳遞給 Fortran 實(shí)作。但是，當(dāng)使用行主矩陣a 呼叫時(shí)，包裝器必須分配記憶體、明確轉(zhuǎn)置並將a 複製到臨時(shí)矩陣a_t、重新計(jì)算前導(dǎo)維度的步幅、使用a_t 呼叫dlaswp、轉(zhuǎn)置並複製儲(chǔ)存在a_t中的結(jié)果到a，最後釋放分配的記憶體。這是相當(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í)作並將輸入?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)榭梢詮奶峁┑牟椒茢囗樞颉?
2. 與傳統(tǒng)的LAPACKE API相比，當(dāng)輸入矩陣為行主矩陣時(shí)，我們不需要將資料複製到臨時(shí)工作區(qū)數(shù)組，從而減少了不必要的記憶體分配。
3. 與直接與BLAS 和LAPACK 互動(dòng)的NumPy 和SciPy 等函式庫相比，在stdlib 中呼叫LAPACK 程式時(shí)，我們不需要之前將不連續(xù)的多維資料複製到臨時(shí)工作區(qū)數(shù)組或從臨時(shí)工作區(qū)數(shù)組複製非連續(xù)的多維資料和呼叫後，分別。除了與硬體優(yōu)化的 BLAS 和 LAPACK 介面時(shí)之外，所採用的方法有助於最大限度地減少資料移動(dòng)並確保資源受限的瀏覽器應(yīng)用程式的效能。

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

目前狀態(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í)的建置問題。
3. 為 stdlib 的現(xiàn)有 LAPACK 套件推出 C 和 Fortran 實(shí)作以及本機(jī)綁定。
4. 繼續(xù)發(fā)展 stdlib 的模組化 LAPACK 程式庫。
5. 確定效能最佳化的其他領(lǐng)域。

雖然我的 Quansight Labs 實(shí)習(xí)已經(jīng)結(jié)束，但我的計(jì)劃是繼續(xù)添加軟體包並推動(dòng)這項(xiàng)工作。鑑於 LAPACK 的巨大潛力和根本重要性，我們很高興看到將 LAPACK 引入網(wǎng)路的這一舉措繼續(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)大、高效能的函式庫來加速您的專案開發(fā)，並讓您安心，因?yàn)槟滥蕾嚨氖蔷难u作的高品質(zhì)軟體。

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

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