大規模線性方程組求解的案例
大規模稀疏矩陣線性方程組求解可以有多快！
對于工業軟件研發尤其是CAE軟件研發來說，線性方程組的求解精度和速度較為重要，在方程組規模上來以后，以Krylov子空間為基礎的PCG，Bicgstab，Gmres等方法相對于直接法在求解效率和精度控制上有較大優勢。這也是abaqus等商業軟件在對于規模較大的模型時采用迭代法求解通常會比直接法求解更快的其中一個原因。
以自由度數為41w的如下分布的稀疏矩陣為例，
該稀疏矩陣形成的系數方程組在matlab中采用直接求解為1.5s,采用一定預處理下的PCG求解時間為0.55s。
用筆者自行開發的稀疏矩陣PCG求解器運行此算例，方程組求解時間為0.628s。
此時，整個程序的運行時間瓶頸實際上并不在方程求解，而在于從文件中讀取稀疏矩陣對應的數據。一般情況下采用fscanf讀取數據會快于用fstream讀取。改用fscanf讀取數據后，程序總運行時間從原來的16s變為9s。
【完】
歡迎關注公眾號 有限元術
一個講有限元技術的公眾號
展開 千萬自由度大規模線性方程組并行求解庫PETSc的介紹
同時，當ksptype為preonly時，表明方程組不采用迭代求解，只采用預處理方法直接求解方程組，因此此時實際上就實現了線性方程組的直接法求解。
【數值算法】系數矩陣非對稱時，線性方程組如何求解？-穩定雙共軛梯度法(Bicgstab)求解線性方程組
enddo
write(*,*)"the solution of equation:"
write(*,"(es18.8)")x
end subroutine bicgstab
依據上述過程編寫程序，計算前述非對稱矩陣線性方程組求解結果：
采用matlab求解該方程組的解：
通過對比可知11次迭代已經獲得即為準確的結果。實際上，對于該方法也可以通過一定的預處理方式，使得其所需要的迭代次數更少。
以上，就是穩定雙共軛梯度法求解線性方程組的內容，感謝您的閱讀！
歡迎關注公眾號 有限元術
展開 說說工業軟件研發的“難點”（1）
其中生產控制類和嵌入式工業軟件歸為工業軟件沒問題，但是把運營管理類和協同集成類也歸為工業軟件就有點離譜了，因為這類公司實在太多了，且很多是上市公司，在國家重視工業軟件，且成為風口的大環境下，不知道是否有意為之。
這里介紹的技術“難點”來源于作者長期研發過程中碰到的問題，包含了軟件研發內容，偏向于實際工程應用。
1. 如何定位仿真和試驗之間的誤差：這是個系統問題，一般而言，仿真結果很難和試驗完全匹配，原因是多方面的，很多時候仿真也不是追求完全和試驗精準匹配，而是得出符合規律的結論。
2. 如何提升超大模型的性能：一直以來，超大模型處理是任何一款軟件容易出現瓶頸的內容，也是公眾號介紹的重點，參見超大模型處理系列文章 工業軟件研發中處理超大模型(9)--ChatGPT介紹的方法 和 一篇文章入門仿真軟件性能優化
3. 軟件架構設計：軟件架構設計是軟件產品的框架性內容，需要經驗豐富的軟件研發經驗。參見CAE軟件架構設計
4. 網格自適應加密相關：自適應網格加密和幾何，網格，求解器緊密相關，不同領域的自適應加密重點不盡相同，也是需要長期積累的技術項。參見 深入理解數值計算網格(全篇)
5. 高效生成六面體：目前為止，商業CFD軟件可以全自動生成六面體網格，而結構分析中很難全自動生成高質量的六面體網格。
6. 大規模線性方程組求解效率和準確性：大規模線性方程組的高效穩定求解仍然是世界性難題，也是公眾號介紹的重點，參見 一篇文章入門大規模線性方程組求解
7. 高性能計算，分布式/并行計算：這塊涉及到各種工具和第三方庫的使用，是需要長期試驗積累的技術項。
8. 幾何清理和修復：現在仍然沒有一個工具能夠完全高效，自動化進行幾何清理和修復。未來也是可以借助AI的內容。參見 深入剖析三維幾何內核(全篇）
9.
展開 
轉貼：變參數非線性方程組的求解！
變參數非線性方程組的求解！
對于求解非線性方程組一般用fsolve命令就可以了，但是對于方程組中某一系數是變化的，該怎么求呢？
%定義方程組如下，其中k為變量
function F = myfun(x,k)
H=0.32;
Pc0=0.23;W=0.18;
F=[Pc0+H*(1+1.5*(x(1)/W-1)-0.5*(x(1)/W-1)^3)-x(2);
x(1)-k*sqrt(x(2))];
%求解過程
H=0.32;
Pc0=0.23;W=0.18;
x0 = [2*W; Pc0+2*H]; % 取初值
options = optimset('Display','off');
k=0:0.01:1; % 變量取值范圍[0 1]
for i=1:1:length(k)
kk=k(i);
x = fsolve(@(x) myfun(x,kk), x0, options);%求解非線性方程組
x1(i)=x(1);
x2(i)=x(2);
end
plot(k,x1,'-b',k,x2,'-r');
xlabel('k')
legend('x1','x2')
[ 本帖最后由 studyboy 于 2006-7-30 17:38 編輯 ]
圖片附件: k-x1.x2.bmp (2006-7-5 23:07, 689.12 K)
展開 【數值算法】共軛梯度法求解線性方程組
在有限元程序開發中，線性方程組的求解是一個重要組成部分。在百萬自由度大規模計算的情況下，線性方程組的高效快速求解對整個求解器的計算效率有著至關重要的作用。無論實際上計算的是線性問題，還是各種非線性問題，最終都需要落實到線性方程組的求解上去。非線性方程組的求解實際上往往就是多次求解線性方程組。
目前，線性方程組的求解主要分為直接法和迭代法兩種。
在之前的文章[數值算法與編程]高斯消去法中，我們討論的高斯消去法就是直接法的一種。而本文即將討論的共軛梯度法，是迭代法的一種，并且，其屬于目前求解對稱線性方程組的主要迭代方法。各大商業有限元軟件，在面臨對稱線性方程組的求解時幾乎都會選用各種變化形式的共軛梯度法進行求解。
共軛梯度法的具體原理和算法如下：
假定要求解的對稱線性方程組是：
其中，A是對稱正定的系數矩陣。
則實際上待求的解也是方程
取得最小值的時候的解。
求該方程的最小值的常見方法是最速下降法，該方法算法偽代碼如下：
該方法實際上是沿著負梯度方向進行搜索，直至殘量接近0，較為簡便，但是在條件數很大時，該方法收斂很慢。
展開 大型稀疏線性方程組求解技術——工業仿真的底層核心
背景
在工業仿真領域，對各種現實世界的問題進行數值模擬時，如流體動力學分析、電磁場仿真、結構力學應力應變分析等，其控制方程通常是偏微分方程組，在經過不同方法的隱式離散之后最終都可轉化為大型稀疏線性方程組。隨著人們對計算精度要求的不斷提高，方程組的階數也從上千階、幾十萬階提高到百萬、千萬階甚至更高，所需的計算量以及存儲需求也隨之迅速膨脹。根據一般經驗，方程組求解時間會占總計算時間的70%以上，往往是整個計算過程中的性能瓶頸。如果說求解器是工業CAE軟件的核心模塊，那么大型稀疏線性方程組的求解技術將毫無疑問是底層求解器的核心。
NASA翼型網格經過離散得到的稀疏矩陣（素材來源于網絡）
方法
眾所周知，稀疏線性方程組的求解方法可以分為直接法和迭代法 ，兩類方法各有優劣，特點比較如下：
迭代法[1]：
對于不同類型稀疏矩陣表現差異較大，存在收斂性與收斂速度問題，催生了許多預處理技術（Preconditioners）；
對原矩陣的編輯很少，SpMV(Sparse matrix-vector multiplication)是其核心運算；
內存需求小，求解速度較快，算法復雜度低；
較易實現并行化。
直接法[2]：
通用、穩定；通過前后處理，能夠保證計算的收斂性與精度；
對原矩陣的編輯多（分解、排序、縮放等）；
內存需求大，求解速度慢，算法復雜度更高；
并行度有限。
其中迭代法的種類很多，可以分為定常（Stationary）迭代法與非定常迭代法[3]。
展開 單元細分源代碼和線性方程組求解器
最近寫程序，需要一個高效簡單的線性方程組求解器，就在別人的基礎上改造了一套代碼。
這套有兩個版本:fortran/C++，主要用于求解稀疏對稱矩陣；
C++ version 來源于網絡。 fortran version 由Esimulate改寫。 fortran version 的輸入文件為C++ version\LDL\Matrix 中任意文件，名字固定為A01。
ldlsolver.rar
技術分享︱大型稀疏線性方程組求解技術——工業仿真的底層核心
一、背景
在工業仿真領域，對各種現實世界的問題進行數值模擬時，如流體動力學分析、電磁場仿真、結構力學應力應變分析等，其控制方程通常是偏微分方程組，在經過不同方法的隱式離散之后最終都可轉化為大型稀疏線性方程組。隨著人們對計算精度要求的不斷提高，方程組的階數也從上千階、幾十萬階提高到百萬、千萬階甚至更高，所需的計算量以及存儲需求也隨之迅速膨脹。根據一般經驗，方程組求解時間會占總計算時間的70%以上，往往是整個計算過程中的性能瓶頸。如果說求解器是工業CAE軟件的核心模塊，那么大型稀疏線性方程組的求解技術將毫無疑問是底層求解器的核心。
NASA翼型網格經過離散得到的稀疏矩陣（素材來源于網絡）
二、方法
眾所周知，稀疏線性方程組的求解方法可以分為直接法和迭代法 ，兩類方法各有優劣，特點比較如下：
迭代法[1]：
對于不同類型稀疏矩陣表現差異較大，存在收斂性與收斂速度問題，催生了許多預處理技術（Preconditioners）；
對原矩陣的編輯很少，SpMV(Sparse matrix-vector multiplication)是其核心運算；
內存需求小，求解速度較快，算法復雜度低；
較易實現并行化。
直接法[2]：
通用、穩定；通過前后處理，能夠保證計算的收斂性與精度；
對原矩陣的編輯多（分解、排序、縮放等）；
內存需求大，求解速度慢，算法復雜度更高；
并行度有限。
其中迭代法的種類很多，可以分為定常（Stationary）迭代法與非定常迭代法[3]。
展開 工業軟件研發中處理超大模型(6)--有限元求解器
從工程學角度看，并行計算的本質是提升性價比，GPU的應用對于數值計算而言是異構架構一大成功。NVIDIA提供的CUDA計算語言成為GPU計算的主流工具。
有限元求解器主要操作對象是矩陣，對于大矩陣計算，涉及到的是矩陣的分塊和分解以及并行計算。矩陣的分塊和分解是兩個不同的概念，注意區分。
求解線性方程組，不管是迭代法還是直接法，最多的運算是矩陣和矩陣乘法，矩陣和向量乘法。乘法的并行計算是各個基礎庫的計算核心功能之一。
圖畫分
對于一個超大矩陣，通常需要分塊，分塊的目的是在已有硬件資源基礎上，讓每塊獲得合理的計算資源，達到并行計算的目的。而并行計算一個核心問題就是負載均衡，即整體上每個計算資源要與計算內容匹配，避免部分硬件資源長期閑置。這就要用到圖畫分工具。
圖論：一般將網格單元抽象成一張圖的頂點，頂點的權重代表網格單元的計算量，網格單元之間的數據依賴關系抽象為圖的邊，邊的權重代表網格單元之間需要交換的數據量，這張圖反映了問題的基本計算量和數據依賴關系，可以作為該問題的負載模型。對于稀疏矩陣，該圖可以描述向量乘法的基本計算量和數據依賴關系，稀疏矩陣向量乘法的負載均衡可以轉化為圖的多路劃分問題。在一些線性方程求解庫中，可以看到相應的圖劃分工具。
5. 線性方程組求解
線性方程組求解是隱式有限元方法求解的核心。考慮到大規模稀疏對稱矩陣線性方程組求解方法的復雜性以及求解庫繁多，作為系列文章，后面專門分多文介紹。先簡要介紹兩個解線性方程組的工具庫，MUMPS和OpenBLAS，具體使用方法不展開。
展開 一篇文章入門“求解器”開發(全篇)
考慮到求解器的數據規模，在設計定義數據之初，須使用數據指針，并避免數據之間深拷貝。
矩陣結構定義
Eigen庫提供了很多基本數據結構，可以直接使用，但從性能，效率以及擴展性看并不最優。如果是商業軟件開發，基礎數據結構最好自己定義。
7. 線性方程組求解
對于很多求解器來說，常用的
FEM,FVM,MOM,BEM等需要求解線性方程組。求解線性方程組一般不會自己寫，而是調用一些已有的求解庫。關于線性方程組的求解也一直是公眾號介紹的重點，參考
兩本線性方程組相關書籍pdf下載
一篇文章入門大規模線性方程組求解
大規模線性方程組解法簡介
除了以上專門介紹的線性方程組解法，在求解器相關內容介紹中，也通常會順帶介紹某些特殊的解法。
關于線性方程組求解庫的選擇，之前有過推薦：
入門Eigen
根據實際情況，可用MKL, OpenBlas，MUMPS
其它一些比如PETSc，Hypre，Trillion等稍重量級的也可以試用
某些專門加速的商業庫也可以使用
如果是需要長期從事線性方程組求解工作的，建議可以從最基礎的求解庫
BLAS,LAPACK,ScaLPAPCK用起，很多使用第三方庫碰到的問題和疑惑，根源都能從這些基礎庫里找到答案。另外很多比較老的庫也可以使用，比如SuiteSparse，SuperLu，Pardiso，TAUCS等，了解每種庫的使用特點以及適用范圍。
方程組求解的準確以及性能，和矩陣的數據特點，矩陣規模，硬件環境緊密相關。需要注意的是，在數據結構上做好封裝，方便調用不同庫之間進行切換。
展開 
智能熱流體仿真軟件AICFD 2023R2新版本功能介紹
此外，新版本還優化了部分湍流模型如Realizable K-E模型，提升求解精度和效率，加快收斂速度。
（a）總壓等值面圖（顏色代表速度大小）
（b）聲壓級頻譜
圖5 某乘用車氣動噪聲仿真分析
4）
優化多相流VOF算法
AICFD 提供多種多相流模型，支持水管理、船舶靜水阻力、混合流動等多相流仿真分析。其中，VOF模型是一種應用于固定歐拉網格的界面跟蹤技術，適合模擬分層或自由表面流動。AICFD 2023R2優化了VOF模型算法，并與實驗數據對比驗證了算法的合理性，提升了計算精度。
（a）體積分數云圖（紅色為水，藍色為空氣）
（b）總阻力仿真結果與實驗對比
圖6 船舶靜水阻力計算
5）
新增多相流VOF空化模型（BETA版），可模擬空化現象
在多相流VOF模型的基礎上，AICFD 2023R2新增Zwart空化模型（BETA版），可用于模擬旋轉機械中由于壓力快速改變引起的空泡產生和破裂現象（即空化現象），并評估空化對旋轉機械效率的影響。
圖7 翼型流動仿真與實驗結果對比，得到的空穴長度相對試驗值的誤差小于10%
6）新增代數多重網格（AMG）預處理策略（BETA版），提升收斂速度
大規模稀疏線性方程組的求解是流體仿真的關鍵技術之一，AICFD采用的是共軛梯度(CG)類迭代法求解技術，并采用合適的預處理策略改善稀疏矩陣的條件數，從而提升迭代法的收斂性和可靠性。
展開 CAE仿真軟件發展趨勢
功能和架構：
從大的方面來說，CAE軟件通常是綜合了數學、物理、電化學等多學科理論基礎、計算機科學、工程應用實踐三方面的知識。對于現實客觀世界的建模通常依靠物理、電化學等基礎學科理論，深入理解描述物理、電化學等過程的原理理論。原理理論的開拓通常需要無數數學家、物理學家、化學家以及科研工作者持續的探索和歸納，不斷提出對客觀規律更為精準的描述方法。電化學過程目前普遍采用Newman提出的偽二維模型來描述，遵循電荷守恒和物料守恒。流體仿真中描述流動的底層物理規律是動量、質量守恒，其中N-S方程被用來描述動量守恒過程，連續性方程描述質量守恒。每一種物理場都可以用若干模型和方程來描述，當涉及到多個物理場時，還要建立耦合模型，構成偏微分方程組再進行求解。
這些現實客觀世界的建模所設立的方程通常為偏微分方程，從這一點來說，偏微分方程是人類用來描述客觀世界的工具，而CAE仿真軟件就是要通過計算機科學技術來實現對客觀世界的建模、求解、結果展示以及優化設計。
CAE軟件從產品角度來看有三個重要的組成部分，分別是前處理器、求解器和后處理器。后處理器用于展示CAE求解結果，這里通常運用計算機圖形學技術來開發。CAE軟件的核心在于對現實現象的精準建模以及實現高保真度的數值模擬求解。CAE軟件的仿真質量和數值模擬的質量直接相關。
數值模擬計算一般分為三個環節：網格離散（網格劃分）、邊界條件設定和求解過程本身。數值模擬求解的核心在于網格離散和數值算法，因此前處理器中的網格劃分模塊和求解器是CAE軟件中最為關鍵的兩個部分。
首先，數值計算需要輸入高質量的網格數據，否則就會是“garbage in, garbage out”，網格剖分的好壞直接決定了仿真質量的高低，因此網格生成的算法技術非常重要。而且網格生成的密度還要結合實際工程應用場景，平衡計算精度和計算效率。
展開 技術分享︱突破大規模CFD仿真瓶頸：UNAP代數求解庫性能實測與優化解析
<p><br></p><p class="ql-align-center"><img class="ztext-gif" width="640" role="presentation" src="https://pic1.zhimg.com/v2-4535bc19aaf1c155e5894f226a8af668_b.webp" data-thumbnail="https://pic1.zhimg.com/v2-4535bc19aaf1c155e5894f226a8af668_b.jpg" data-size="normal" alt="動圖" style="display: block;"></p><h2><strong>01 引言：CFD大規模仿真效率瓶頸</strong></h2><h3><strong>1.1 </strong><strong style="background-color: rgba(1, 0, 0, 0);">背景與目標</strong></h3><p> 隨著工業數值模擬在航空航天、能源、船舶及裝備制造等領域的廣泛應用，對計算流體力學（Computational Fluid Dynamics, CFD）軟件的性能要求日益提升。某基于有限體積方法的通用 CFD 仿真軟件，能夠模擬穩態與瞬態下的復雜流動，并支持旋轉機械流動分析與先進的共軛熱傳遞（CHT）模型。結合特定的前處理器，該軟件可以快速生成計算網格，從而為用戶提供從前處理到數值求解的完整解決方案。</p><p> 然而，隨著工程問題規模的不斷擴大，傳統代數求解器在處理大規模稀疏線性方程組時暴露出性能瓶頸。
展開 