一條聽起來很怪的規則——而且提出的人明明寫得出來
你已經花了好幾級在學那套機器：帶選軸的 LU 分解、求特徵值的 QR 演算法、O(N log N) 的 FFT、解超大稀疏系統的共軛梯度法。所以當有人告訴你專業科學計算的第一條規則時，可能會覺得像被背叛了：別自己重造輪子。重點不是你做不到——你現在是真的懂這些方法了——而是在正式產品程式裡重新實作它們，幾乎永遠是把力氣用錯了地方，而且暗藏危險。把標準函式庫蓋起來的那群人喊這句話喊得最大聲，因為他們最清楚：要把細節微妙地搞錯，有多少種方式。
把它想成外科手術。懂解剖學是必要的——但懂解剖學不等於決定要替自己開刀。你為了學會而寫的那個課本版高斯消去法，是很好的教學工具，也是很好的雛形。但 LAPACK 裡頭那個版本，已經被四十年來的數值分析家、編譯器專家和硬體工程師把每個角落都拿去極限測試過。這道落差，就是本篇要談的全部。
這個堆疊，一層一層看
你真正站立的東西，是一個數值軟體堆疊：一座由各層疊成的塔，每一層都信任它下面那一層，並把那層對上面的層隱藏起來。最底層坐著 BLAS，也就是 Basic Linear Algebra Subprograms（基本線性代數副程式）——一小套向量與矩陣運算的「詞彙」（內積、矩陣-向量乘積、矩陣-矩陣乘積），其他一切都用這套詞彙來表達。再上一層是 LAPACK，它把那些大演算法——解 A x = b、最小平方、特徵值、SVD——完全寫成一連串的 BLAS 呼叫。再上面才是你實際打字的那些親切環境：Python 裡的 NumPy 和 SciPy、Julia 內建的陣列語言、MATLAB、R。當你在它們任何一個裡頭寫下一行求解時，你常常是直接往下呼叫到底下同一套 LAPACK 和 BLAS。
your code: numpy.linalg.solve(A, b) <- one line you write ------------------------------------------------------------- SciPy / NumPy / Julia / MATLAB <- arrays, slicing, plots LAPACK dgesv, dgeqrf, dsyev, dgesdd <- the algorithms BLAS ddot, dgemv, dgemm <- the vocabulary hardware SIMD registers, cache, cores <- where flops happen
這樣切分有一個很漂亮的工程理由。BLAS 分成三個等級：第 1 級是向量對向量的運算（資料量 O(n)、浮點運算量 O(n)）；第 2 級是矩陣對向量（兩者都是 O(n^2)）；第 3 級是矩陣對矩陣（資料量 O(n^2)，但浮點運算量 O(n^3)）。只有第 3 級每從記憶體抓一個數，就做了夠多算術，足以讓現代處理器忙起來——回想 HPC 那一級講過，一個受限於記憶體的核心，瓶頸是快取未命中而非浮點運算量。所以 LAPACK 被刻意改寫成盡量把時間花在第 3 級的矩陣-矩陣呼叫裡。廠商接著為自家晶片出一套手工調校的 BLAS，同一份 LAPACK 程式碼瞬間就快了十倍。這道乾淨的接縫——上面是演算法、下面是隨硬體而變的速度——就是這個切分能存活四十年的原因。
四十年的抓蟲，到底替你買到了什麼
別自己重造輪子的最深理由不是速度——而是你先前學到的那種「浮點意義下的正確性」。回想整個領域的主方程式：準確度 = 條件 × 穩定度。條件屬於問題本身，你改不了它；穩定度屬於演算法，而且它殘酷地容易搞丟。像 LAPACK 求解器那樣的函式庫常式是後向穩定的：它回傳的答案，恰好是某個與你提出的問題只差一個捨入誤差距離的問題的精確解。這是有限精度演算法能提供的最強誠實保證，而要證明它得靠真正的研究。你那個跳過選軸、或加總順序隨便的手寫迴圈，給不了你這種保證——而且它是默默地失敗的，回傳一個看起來合理卻是錯的答案，而不是當掉。
不過要誠實面對這個界線，因為初學者最容易在這裡過度信任函式庫。後向穩定並不代表答案準確。如果你的問題是病態的——比方說矩陣的條件數接近 10^8——那麼即使是一個完美、後向穩定的求解器，也會損失你約 16 位雙精度數字中的大約 8 位，純粹是因為這個問題會放大那無可避免的輸入捨入。函式庫把它的工作做得無懈可擊；是問題本身很難。好的函式庫會把條件數交到你手上，讓你能預見這件事；手寫的求解器通常讓你對此完全看不見。函式庫替你買到的是一個值得信賴的演算法，而不是一個值得信賴的問題。
函式庫止步、你的判斷接手之處
「別自己重造輪子」是一條關於「實作」的規則，而不是關於「理解」的——而它有實實在在的邊界。函式庫會把你字面上問的那個問題的正確答案交給你，但它無法知道那是不是「對的問題」，也無法知道你餵給它的是不是垃圾。你學過的所有陷阱，在函式庫這條線之上仍然成立。如果你要對等距點做高次多項式擬合，再完美的函式庫也擺脫不了 Runge 現象，擬合會在邊緣發散；解方是改選 Chebyshev 節點或樣條，這是個函式庫沒辦法替你做的建模決定。如果你以為把有限差分的步長 h 無限縮小就能更準，捨入下限終究會讓答案變差，有沒有函式庫都一樣。
也確實有該自己寫的時候。如果你的問題有某種沒有函式庫能呈現的特殊結構——一種不尋常的稀疏樣式、一個自訂的預條件子、一個你只會「作用」卻不知怎麼「儲存」的無矩陣運算子——那麼你大可名正言順地自己組外層迴圈。但請注意這門手藝：你只蓋那個新穎的部分，而裡頭每一個標準步驟仍然呼叫 BLAS、LAPACK 或某個稀疏套件。你做的是外殼，從不重造核心。而當你這麼做時，你會用對待陌生人程式碼一樣的懷疑來對待自己的程式——這正是下一篇的主題「驗證」存在的意義。
如何在不迷失主線的情況下挑工具
你最常碰到的三種環境，骨子裡是同一具引擎，握在手上的感覺卻不同。NumPy 和 SciPy 給你 Python 龐大的生態系，外加包在 C 與 Fortran 函式庫外面那層又薄又快的封裝——你大部分時間都花在編譯過的程式碼裡，所以只要你讓運算對整個陣列向量化、而不是逐元素跑迴圈，那個慢吞吞的 Python 直譯器就很少咬到你。Julia 是為這類工作從頭設計的：它把你自己寫的高階程式碼即時編譯成接近 C 的速度，於是「慢的雛形語言」和「快的正式語言」之間那道老舊而痛苦的分裂就被軟化了。MATLAB 是經典的數值工作台，打磨精良、教學普及，而它坐的也正是同一套 LAPACK。沒有放諸四海皆準的最佳選擇——只有對你的團隊、你的函式庫和你的問題而言對的那一個。
把這些收攏起來，這條規則就不再像背叛了。你學那些演算法，是為了知道該伸手拿哪個函式庫常式、它的結果代表什麼、以及它會在哪裡失敗——而不是為了把它們重打一遍。本級其餘的內容正是站在這個基礎上：知道你腳下的堆疊值得信賴，你就能把稀缺的注意力，花在函式庫永遠不會替你回答的那些問題上——我有沒有把方程式解對？我有沒有解對的方程式？有沒有人能重現這個結果？以及，我誠實的誤差棒是多少？