網路瀏覽器中的 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 應用程式的記憶體佔用只會增加。
最後，雖然 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=⎡⎡11111 352466666666 ⎤⎤，B=⎡000

0 00 >0000> >⎦⎤​

當 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}1{\dots }_{n\cdot \mathrm{}s}-$$1textrm{idx} = textrm{偏移} i cdot s_0 j cdot s_1 ldots 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中文網其他相關文章！