面向對象在數值計算中如何應用，有哪些具體例子，相較普通函數能夠多大程度提高計算效率？

01-27

面向對象在數值計算中如何應用，有哪些具體例子，相較普通函數能夠多大程度提高計算效率？
另外有沒有像《C++和面向對象數值計算》這本書一樣，詳細介紹 C++ 各種特性在數值計算中應用，並且能夠提高計算效率技巧的書籍？

挺有意思的問題，我也是動手寫了一些大程序才逐漸邊開發邊學習到c++的一些有趣特性。拋磚引玉。不過學術圈絕大多數情況下並沒有IT大廠那樣的老手code review系統，我感覺靠自己折騰進步挺慢的。。。

以老牌線性代數庫Eigen為例，用到的一些特性比如

1。表達式模板與lazy evaluation

相比於手工寫若干BLAS、LAPACK函數，這個特性可以允許做更多的優化

Lazy Evaluation and Aliasing

2。模板類特化

解決不同容器std::allocator分配對象的內存對齊問題，對std::vector 再特化一個。SIMD指令完全依賴這個。

Using STL Containers with Eigen

3。模板與繼承的嵌套（我c++水平比較菜還不能完全理解這麼做的優勢）

比如Matrix類的聲明是這樣的。用戶可以方便的拓展，比如定義自己的permutation matrix、band matrix。

template& class Matrix : public PlainObjectBase& &>

4。模板實現編譯時尺寸已知的矩陣塊操作

Block operations

另外一些庫的例子比如

1。用模板實現的struct of array：

Lunarsong/StructureOfArrays

2。在分散式線性代數容器Tpetra里，定義了一個虛基類TpetraOperator，用戶通過繼承這個虛基類的介面可以非常方便的實現任何matrix free 迭代求解器。並且所有傳遞都可以通過智能指針實現。

https://trilinos.org/docs/r12.12/packages/tpetra/doc/html/classTpetra_1_1Operator.html

3。在FMM庫里實現一個多線程內存池，取代各線程臨時的malloc和free，為SIMD計算部分返回對齊內存：

https://github.com/dmalhotra/pvfmm/blob/develop/include/mem_mgr.hpp

4。用模板實現的SIMD庫：

NumScale/boost.simd

5。比如FFTW庫本身的create_plan並不線程安全，如果用一個Singleton模式的wrapper把所有FFT相關的非線程安全部分都包裹起來加上鎖，各個線程就可以安全的並行執行execute(plan)。如果不用這樣的模式，相關的非線程安全部分散落在各個子函數里，即使要加鎖，在哪裡加鎖合適也不是一個簡單的問題。（我自己的小項目，施工中）。

6。比如涉及MPI傳遞數據結構的部分，可以使用type_traits 里的std::is_trivially_copyable檢查要傳遞的對象類型，結果為真則可以直接memcpy到MPI緩衝區，結果為假再去調用類型的序列化模塊。

簡單來說我的感覺就是，C++相比C和Fortran在多線程、模板、繼承這幾方面的進步允許用戶創建非常靈活的代碼結構和介面，有些用C/Fortran需要一些技巧才能實現的優化在C++上想起來、寫起來都比較自然。

在我的理解中OOP主要是方便程序員的，很多情況下大家討論的都是使用OOP的性能懲罰，不過如果直接說結論的話，不大量使用虛函數的話，理論上這個懲罰多數情況下是可以忽略的。我個人的一點點測試發現同樣的代碼可能g++要比gcc慢10%左右，或許和編譯器的具體實現有關係，但很多時候收益還是更大一點。總體而言C++的特性我可能用得還是比較少，我的代碼風格還是更貼近C，其實純C也可以寫得很像OOP，不過比C++費事很多。前段時間我看了比較多關於性能優化的manual，希望後文說的這點東西能夠有些作用。當然我的計算物理基本上是野生的，所以如果沒有太多參考價值也不用驚訝。

對C++的性能優化我推薦看看這本手冊：

http://www.agner.org/optimize/optimizing_cpp.pdfwww.agner.org

作者講了不少東西。我目前階段的認識是：

1、絕大多數優化編譯器都幫你操心了。包括但不限於常數傳遞、表達式化簡、自動展開小的for循環、自動向量化等等。

2、但是你要寫出能方便編譯器優化的程序。一個最大的阻礙是指針的aliasing。一個例子：

double Foo(double* a, double* b) { return *a+*b; }

這裡a和b兩個指針，在語法上有指向同一塊內存（aliasing）的可能性，但是實際工程中很少有人會這麼使用。然而，編譯器並不知道這一點，所以它會考慮最壞的可能性，不去作更激進的優化來保證語法上的最大兼容性。然而，我們可以向編譯器指明，這兩個指針永遠不會用於指向同一塊內存，來允許編譯器進一步優化。辦法是使用restrict關鍵字。在不同編譯器里具體用法不一樣，在gnu c里可以寫

double Foo(double*__restrict__ a, double*__restrict__ b) { return *a+*b; }

在VS里可以寫

double Foo(double*__restrict a, double*__restrict b) { return *a+*b; }

3、另一個阻礙編譯器優化的因素是浮點數運算。因為浮點數運算的截斷誤差，很多代數上成立的數學恆等式給出的計算結果是不一樣的。（比如，1+1e-16-1和1-1+1e-16是不一樣的），所以很多編譯器傾向於不去優化浮點數運算。然而浮點數計算比整數慢很多（除法又比乘法慢很多），所以這部分優化要程序員自己操心。能用整數運算來代替浮點數運算是最好的情況（比如格點模擬有很多東西都可以rescale成整數），但頻繁的整數/浮點數轉換也是不好的。

4、cache有可能是另一個非常關鍵的性能瓶頸。因為cache line有限，你總是希望一次讀取儘可能多的數據，因此把一起使用的數據組織在一起總是好的。這倒是OOP的優點——一個對象的私有數據總是在一起的。因此一個class str{double a; double b}; str s[1000];要比兩個對象double a[1000], b[1000]更好。另外cache line的設計可能會使array在某個size造成嚴重的cache conflict，這是需要避免的。但相關的計算我還不是很熟悉。

5、現代CPU的一個重要進步就是流水線設計，這意味著前一個指令還沒有完成時，後一個指令就可以開始執行了，甚至CPU可以在一定程度上亂序執行，但前提是你前後的命令沒有相關性。如果後一個指令總是依賴前一個指令的結果，就會浪費一些性能。所以代碼盡量減少前後相關性是要好一點的。比如：

for(int i=0; i&<1000; i++) a[i]=b[i]+c[i];

可能比

for(int i=0; i&<1000; i+=2) { a[i]=b[i]+c[i]; a[i+1]=b[i+1]+c[i+1]; }

要慢一點。

6、於此相應的，在你作條件判斷時，為了讓流水線能夠繼續運轉，CPU會去預測你的判斷結果來加快運行速度。如果條件判斷的正確性是有規律的，CPU根據歷史作出的預測絕大多數都是對的，那麼條件判斷的速度會很快。然而，如果條件判斷完全隨機沒有規律，那麼發現自己預測錯誤的CPU會花非常多時間（12-25個cycle對比0-2個cycle）進行回滾操作，嚴重拖慢運行速度，稱之為branch misprediction。複雜的for循環也會造成這個問題，所以是需要避免的。具體方法五花八門，有時候省掉條件判斷改用look up table可能會更好。

6.5、據說數據內存地址的對齊對性能也有影響，但我現在還不是很確定影響有多大。

7、C++的各種OOP特性通常不會造成性能懲罰，只是過多的繼承可能會讓數據結構非常臃腫，影響cache的效率。然而，虛函數的使用需要CPU花一些cycle去查表尋找正確的實現，如果總是切換實現也可能會造成branch misprediction拖慢性能。

8、然而templete可以作很多黑科技，甚至在一些情況下用來代替虛函數。其原理在於templete是在編譯的時間被解析的，這使得很多原本運行時的branch在編譯期就被消除了。但這類黑魔法我還不會用。

8.5、還有些黑魔法是基於位運算的速度比其他要快。比如

float x; *(int*)x |= 0x80000000; // Set sign bit of x

是一種快速計算-abs(x)的黑魔法。

（注意（引自手冊原文）：

The trick violates the strict aliasing rule of standard C, specifying that two pointers of different types cannot point to the same object (except for char pointers). An optimizing compiler might store the floating point and integer representations in two different registers. You need to check if the compiler does what you want it to. It is safer to use a union, as in example 14.23 page 146.
The trick will fail if the object is treated as bigger than it actually is. This above code will fail if an int uses more bits than a float. (Both use 32 bits in x86 systems).
If you access part of a variable, for example 32 bits of a 64-bit double, then the code will not be portable to platforms that use big endian storage.
If you access a variable in parts, for example if you write a 64-bit double 32 bits at a time, then the code is likely to execute slower than intended because of a store forwarding delay in the CPU (See manual 3: "The microarchitecture of Intel, AMD and VIA CPUs").

）

9、最後也是最重要的是：優化性能之前應當先知道性能的瓶頸在哪裡，不能頭痛醫腳。通常演算法的影響是最大的，只有演算法沒有改進餘地了才會去考慮這些。很多時候多做benchmark總是好的。要profile程序的性能，VS有成熟的工具，Linux下可以使用perf，非常強大。

最後推薦一下同一個作者寫的這個C++的數學庫：

http://www.agner.org/optimize/vectorclass.pdfwww.agner.org

它應該是用到了不少C++的特性。在我的模擬中發現libm的數學函數計算非常非常非常非常耗時間，但這個庫可以利用一些比較先進的指令集（比如SSE和AVX）來優化數學函數的計算。這是因為現代CPU的重要進步除了流水線的設計外，還加入了向量化運算的功能，可以在一個指令內同時做多個數據的運算。比如SSE2以上的指令集支持128bit的向量運算，AVX支持256bit，AVX512甚至支持512bit。也就是說如果你使用支持AVX指令集的CPU，你可以同時進行4個double的運算，這也可以極大加快運行速度。

如果想進一步利用多核處理器的性能，甚至使用顯卡計算，那就是另一個話題了。多CPU並行計算，在科學計算中常用的有多線程共享內存的openMP和多進程不共享內存的MPI。具體教程參見

https://computing.llnl.gov/tutorials/openMP/computing.llnl.gov

https://computing.llnl.gov/tutorials/mpi/computing.llnl.gov

顯卡計算我用了很多CUDA，體驗嘛……………………快的時候是真快，debug的時候是真噁心。

我覺得是沒什麼幫助的……

很多答主都講了很多，但依照題目的限定，面向對象應該不是重點。

lazy evaluation、模板特化，但這不是OOP的內容啊……Java夠OOP吧，Java卻做不到啊。

C++的高性能運算庫依賴於很多C++的特性與編譯期優化。

至於CPU指令集，那更不是C++專利。

通過繼承介面（虛基類）來實現運行時多態，擴展功能，這算是利用了OOP，不過C里傳指針也是慣用方法吧（不過需要對外暴露設計細節），很多C的庫就會提供接入內存分配/IO操作的渠道（設置回調）。

依我看，高性能計算庫應該是面向數據的，對象和對象的交互更弱。

真要說什麼幫助，大概是「繼承」能節省一些代碼量吧，跟模板特化道理差不多。

舉個例子，向量和矩陣相乘，有左乘和右乘，所以這兩方法是屬於誰的？兩者都寫一份相同的方法？還是說藉助向量派生於矩陣來規避這個問題（節省代碼量），然後在運算代碼里加入針對一維向量的特殊優化？

一般來說，設計一個額外的函數來做運算更合理吧？這不就很不OOP嘛……