稀疏矩陣的案例
基于 MATLAB 的 ANSYS Harwell-Boeing 格式稀疏矩陣提取工具 —— 剛度矩陣與質(zhì)量矩陣 ￥30
4.得到K稀疏矩陣后可通過spy(K)來進(jìn)行稀疏矩陣的結(jié)構(gòu)可視化，以及通過命令K_full=full(K)將稀疏矩陣轉(zhuǎn)換為完整（密集）矩陣。
從單元連接關(guān)系到節(jié)點(diǎn)鄰接點(diǎn)-有限元采用稀疏矩陣求解的前置工作
在有限元求解中，最終通常要求解的是一個(gè)關(guān)于場變量的線性方程組，在常見的位移場有限元中，要求解的是各個(gè)節(jié)點(diǎn)的位移，該線性方程組的系數(shù)矩陣通常稱為剛度矩陣，方程組右邊通常稱為右端項(xiàng)或者荷載向量。一般情況下，由于網(wǎng)格劃分后并不是所有節(jié)點(diǎn)都兩兩連接，因此實(shí)際上最終形成的整體剛度矩陣中大部分元素為0，這種矩陣稱為稀疏矩陣。在有限元求解中，對于這種系數(shù)矩陣為稀疏矩陣的方程組，一種常見的方法是僅保存剛度矩陣的非0元素到內(nèi)存中，0元素不保存，這樣就可以以更小的內(nèi)存保存大型結(jié)構(gòu)的剛度矩陣。
那具體矩陣中有多少元素為0，就可以認(rèn)為其是稀疏矩陣呢？這個(gè)界限實(shí)際上比較模糊，有文獻(xiàn)給出如下定義：如果矩陣的A的非0元素?cái)?shù)量為O(n)，其中n是A的階數(shù)，則矩陣為稀疏矩陣。
稀疏矩陣經(jīng)常通過非0元素分布圖表示其稀疏性質(zhì)，以下是兩個(gè)常見的稀疏矩陣的分布圖：
在有限元分析中，非0元素的分布，實(shí)際上主要取決于單元的節(jié)點(diǎn)連接，以下圖中的單元連接為例：
假設(shè)圖中每個(gè)節(jié)點(diǎn)一個(gè)自由度，則整體剛度矩陣為16x16的矩陣，而具體非0元素的分布，可以通過單元連接得到鄰接點(diǎn)得到，所謂鄰接點(diǎn)，指的是相對于當(dāng)前單元位于同一單元內(nèi)的所有點(diǎn)的集合。以節(jié)點(diǎn)6為例，其鄰接點(diǎn)是1,2,3,5,7,9,10,11。
獲得上述鄰接點(diǎn)后，剛度矩陣中第6行的非0元素的位置實(shí)際上就確定了：k(6,1),k(6,2),k(6,3),k(6,5),k(6,6),k(6,7),k(6,9),k(6,10),k(6,11)。
在實(shí)際采用稀疏矩陣求解有限元問題時(shí)，獲得上述非0元素位置后，就可以對剛度矩陣采用稀疏矩陣存儲(chǔ)，常見的存儲(chǔ)方式有COO，CSR，CSC和DIA等。
展開 從單元連接關(guān)系到節(jié)點(diǎn)鄰接點(diǎn)-有限元中形成稀疏矩陣求解的前置工作
在有限元求解中，最終通常要求解的是一個(gè)關(guān)于場變量的線性方程組，在常見的位移場有限元中，要求解的是各個(gè)節(jié)點(diǎn)的位移，該線性方程組的系數(shù)矩陣通常稱為剛度矩陣，方程組右邊通常稱為右端項(xiàng)或者荷載向量。一般情況下，由于網(wǎng)格劃分后并不是所有節(jié)點(diǎn)都兩兩連接，因此實(shí)際上最終形成的整體剛度矩陣中大部分元素為0，這種矩陣稱為稀疏矩陣。在有限元求解中，對于這種系數(shù)矩陣為稀疏矩陣的方程組，一種常見的方法是僅保存剛度矩陣的非0元素到內(nèi)存中，0元素不保存，這樣就可以以更小的內(nèi)存保存大型結(jié)構(gòu)的剛度矩陣。
那具體矩陣中有多少元素為0，就可以認(rèn)為其是稀疏矩陣呢？這個(gè)界限實(shí)際上比較模糊，有文獻(xiàn)給出如下定義：如果矩陣的A的非0元素?cái)?shù)量為O(n)，其中n是A的階數(shù)，則矩陣為稀疏矩陣。
稀疏矩陣經(jīng)常通過非0元素分布圖表示其稀疏性質(zhì)，以下是兩個(gè)常見的稀疏矩陣的分布圖：
在有限元分析中，非0元素的分布，實(shí)際上主要取決于單元的節(jié)點(diǎn)連接，以下圖中的單元連接為例：
假設(shè)圖中每個(gè)節(jié)點(diǎn)一個(gè)自由度，則整體剛度矩陣為16x16的矩陣，而具體非0元素的分布，可以通過單元連接得到鄰接點(diǎn)得到，所謂鄰接點(diǎn)，指的是相對于當(dāng)前單元位于同一單元內(nèi)的所有點(diǎn)的集合。以節(jié)點(diǎn)6為例，其鄰接點(diǎn)是1,2,3,5,7,9,10,11。
獲得上述鄰接點(diǎn)后，剛度矩陣中第6行的非0元素的位置實(shí)際上就確定了：k(6,1),k(6,2),k(6,3),k(6,5),k(6,6),k(6,7),k(6,9),k(6,10),k(6,11)。
在實(shí)際采用稀疏矩陣求解有限元問題時(shí)，獲得上述非0元素位置后，就可以對剛度矩陣采用稀疏矩陣存儲(chǔ)，常見的存儲(chǔ)方式有COO，CSR，CSC和DIA等。
展開 大規(guī)模稀疏矩陣線性方程組求解可以有多快！
以自由度數(shù)為41w的如下分布的稀疏矩陣為例，
該稀疏矩陣形成的系數(shù)方程組在matlab中采用直接求解為1.5s,采用一定預(yù)處理下的PCG求解時(shí)間為0.55s。
用筆者自行開發(fā)的稀疏矩陣PCG求解器運(yùn)行此算例，方程組求解時(shí)間為0.628s。
此時(shí)，整個(gè)程序的運(yùn)行時(shí)間瓶頸實(shí)際上并不在方程求解，而在于從文件中讀取稀疏矩陣對應(yīng)的數(shù)據(jù)。一般情況下采用fscanf讀取數(shù)據(jù)會(huì)快于用fstream讀取。改用fscanf讀取數(shù)據(jù)后，程序總運(yùn)行時(shí)間從原來的16s變?yōu)?s。
【完】
歡迎關(guān)注公眾號(hào) 有限元術(shù)
一個(gè)講有限元技術(shù)的公眾號(hào)
展開 
稀疏矩陣高性能求解器PARDISO的使用
PARDISO是一個(gè)商用的高速的稀疏矩陣線性方程組直接法求解器，其可以用于大型稀疏對稱或非對稱線性方程組的求解，商用版支持的并行方式包括共享內(nèi)存，分布式內(nèi)存和NVIDA的GPU并行。多次實(shí)踐已經(jīng)表明，PARDISO求解器具有優(yōu)秀的性能，在各種稀疏矩陣直接求解器的比較中通常處于第一梯隊(duì)，其他常見的直接法求解器包括SuiteSparse，MUMPS，SPOOLES等，商用有限元軟件COMSOL的官方文檔甚至直接指出在COMSOL中PARDISO求解速度快于MUMPS。。
對于線性方程組求解來說，一般常見的求解方法可分為直接法和迭代法。直接法常見的方法是高斯消去，LU分解法等，在求解過程中不需要迭代，但是在求解中可能需要較大內(nèi)存，商業(yè)軟件中經(jīng)常使用的一種直接法叫做多波前法。迭代法又可分為兩個(gè)大類，分裂迭代和Krylov迭代。前者通過將系數(shù)矩陣進(jìn)行拆分形成迭代列式，常見的是雅克比迭代，高斯賽德爾迭代等，分裂迭代這一類型的迭代法在商業(yè)軟件中很少直接使用，但是經(jīng)常作為Krylov迭代的預(yù)處理基礎(chǔ)；后者以Krylov子空間為基礎(chǔ)上，常見的包括共軛梯度法，廣義最小殘差法等，Krylov迭代法在商業(yè)軟件中較為常見。
在稀疏矩陣直接法求解中，一個(gè)常見的技術(shù)是矩陣重排，通過矩陣重排，可以使得矩陣的帶寬減小（即非0元素位置相對更趨近于對角線），從而使得直接求解效率更高，一種常見的矩陣重排方法是Cuthill_McKee，常見的可以實(shí)現(xiàn)矩陣重排的庫是METIS。以下是采用Cuthill_McKee重排前和重排后的非0元素分布。
重排前非0元素分布
重排后非0元素分布
商用版本的PARDISO個(gè)人一般較少直接使用（主要是要買）。實(shí)際中較為廣泛使用的其中一個(gè)版本是配置在MKL庫的PARDISO版本。
展開 工業(yè)軟件研發(fā)中處理超大模型(6)--有限元求解器
CSR和CSC數(shù)據(jù)結(jié)構(gòu)可以進(jìn)一步壓縮數(shù)據(jù)，是最常用的稀疏矩陣存儲(chǔ)方式。
以上圖片來源于網(wǎng)絡(luò)
以CSR為例， 在上圖4*4矩陣中：
row offsets表示行偏移值：
第一行從0開始，偏移值為0；
第二行，之前有兩個(gè)值1和7，因此偏移值為2；
同理，第三行，前面有1，7，2，8四個(gè)數(shù)，第三個(gè)偏移為4
column indices和values則分別表示列索引和對應(yīng)值：
比如按順序，第0列的第一個(gè)數(shù)為1， 第1列的第一個(gè)為7，第二個(gè)為2。以此類推。
可以看出CSR是一種全局編碼，相比COO，進(jìn)一步減少了行索引值占用數(shù)據(jù)。CSC數(shù)據(jù)結(jié)構(gòu)類似。
3. DIA(Diagonal Sparse Matrix)
即對角稀疏矩陣格式
對于矩陣元分布在主對角線以及主對角線附近的稀疏矩陣，采取這種存儲(chǔ)方式比較有效，但通用性不強(qiáng)。
還有其它很多稀疏矩陣存儲(chǔ)方式，一般是以上幾種方式的衍生，擴(kuò)充或者組合。
稀疏矩陣的數(shù)據(jù)結(jié)構(gòu)定義
矩陣數(shù)據(jù)結(jié)構(gòu)是最基礎(chǔ)的數(shù)據(jù)結(jié)構(gòu)。從計(jì)算角度看，稀疏矩陣一般單獨(dú)定義，單獨(dú)實(shí)現(xiàn)接口；但是從平臺(tái)和系統(tǒng)架構(gòu)看，最好是將數(shù)據(jù)結(jié)構(gòu)融入到系統(tǒng)設(shè)計(jì)中。
看幾個(gè)關(guān)于矩陣計(jì)算的需求：
普通矩陣和稀疏矩陣進(jìn)行加法計(jì)算
將自定義的稀疏矩陣數(shù)據(jù)傳到不同的線性方程組求解庫求解
稀疏矩陣求逆
4.按照業(yè)務(wù)封裝成模塊供其它模塊調(diào)用
這些需求從計(jì)算角度看是要盡量避免的，但是從業(yè)務(wù)和程序邏輯看，又是合理和很可能出現(xiàn)的。
展開 UNAT加速庫：突破異構(gòu)計(jì)算瓶頸，實(shí)現(xiàn)跨平臺(tái)高效仿真
例如，非結(jié)構(gòu)網(wǎng)格鄰接矩陣采用分圖法，結(jié)構(gòu)網(wǎng)格采用二維或者三維分塊法。
分圖法：是一種用于將圖（Graph）或圖相關(guān)的數(shù)據(jù)結(jié)構(gòu)進(jìn)行劃分的方法，它可以將大規(guī)模的圖或網(wǎng)絡(luò)劃分成若干個(gè)子圖，以便在不同的計(jì)算節(jié)點(diǎn)上并行處理。分圖法的目的是盡可能均衡地劃分圖的節(jié)點(diǎn)和邊，以減少通信開銷并提高計(jì)算效率。
分塊法：用于劃分結(jié)構(gòu)化數(shù)據(jù)（如矩陣），在分塊法中，數(shù)據(jù)被劃分成若干個(gè)塊，每個(gè)塊包含一定數(shù)量的數(shù)據(jù)元素。這種劃分可以幫助提高數(shù)據(jù)訪問的局部性，減少緩存失效和提高計(jì)算效率。
一維兩色染色
二維四色染色（考慮對角也有臨近關(guān)系）
四、試結(jié)果展示
01算子性能測試
在算子性能測試中，我們選取了單向量累加、向量乘加、向量取倒數(shù)、計(jì)算調(diào)和級(jí)數(shù)4個(gè)算子進(jìn)行測試。此次測試是在單核環(huán)境下完成的，以確保測試結(jié)果的準(zhǔn)確性和可靠性。
我們對4個(gè)算子在神威平臺(tái)進(jìn)行了測試，測試結(jié)果如下圖所示：
神威平臺(tái)下算子加速比圖
我們可以看出單核環(huán)境下使用算子和不使用算子之間存在顯著的性能差異：4個(gè)算子取得了較高的加速比，加速比均在15以上，其中，調(diào)和級(jí)數(shù)的優(yōu)化效果最為明顯，達(dá)到了160的加速比。
02稀疏矩陣向量乘
在數(shù)值模擬中，求解大規(guī)模稀疏線性方程組是非常重要的一個(gè)環(huán)節(jié)。在迭代求解過程中，稀疏矩陣向量乘法是耗時(shí)最長的計(jì)算核心之一，存在嚴(yán)重的數(shù)據(jù)局部性差、寫沖突、負(fù)載不均衡等問題。因此，稀疏矩陣向量乘法已經(jīng)成為了當(dāng)前性能優(yōu)化的難點(diǎn)和研究熱點(diǎn)。
我們基于神威平臺(tái)，隨機(jī)選取了來自電磁、流體等領(lǐng)域的典型稀疏矩陣算例，測試了近60種矩陣，運(yùn)行結(jié)果如下圖所示：
根據(jù)上圖，我們可看出大部分算例取得了較好的加速比，然而，也有部分算例的加速比并不理想，這可能與內(nèi)部非零元素的數(shù)量以及矩陣的規(guī)模有關(guān)。
展開 大型稀疏線性方程組求解技術(shù)——工業(yè)仿真的底層核心
背景
在工業(yè)仿真領(lǐng)域，對各種現(xiàn)實(shí)世界的問題進(jìn)行數(shù)值模擬時(shí)，如流體動(dòng)力學(xué)分析、電磁場仿真、結(jié)構(gòu)力學(xué)應(yīng)力應(yīng)變分析等，其控制方程通常是偏微分方程組，在經(jīng)過不同方法的隱式離散之后最終都可轉(zhuǎn)化為大型稀疏線性方程組。隨著人們對計(jì)算精度要求的不斷提高，方程組的階數(shù)也從上千階、幾十萬階提高到百萬、千萬階甚至更高，所需的計(jì)算量以及存儲(chǔ)需求也隨之迅速膨脹。根據(jù)一般經(jīng)驗(yàn)，方程組求解時(shí)間會(huì)占總計(jì)算時(shí)間的70%以上，往往是整個(gè)計(jì)算過程中的性能瓶頸。如果說求解器是工業(yè)CAE軟件的核心模塊，那么大型稀疏線性方程組的求解技術(shù)將毫無疑問是底層求解器的核心。
NASA翼型網(wǎng)格經(jīng)過離散得到的稀疏矩陣（素材來源于網(wǎng)絡(luò)）
方法
眾所周知，稀疏線性方程組的求解方法可以分為直接法和迭代法 ，兩類方法各有優(yōu)劣，特點(diǎn)比較如下：
迭代法[1]：
對于不同類型稀疏矩陣表現(xiàn)差異較大，存在收斂性與收斂速度問題，催生了許多預(yù)處理技術(shù)（Preconditioners）；
對原矩陣的編輯很少，SpMV(Sparse matrix-vector multiplication)是其核心運(yùn)算；
內(nèi)存需求小，求解速度較快，算法復(fù)雜度低；
較易實(shí)現(xiàn)并行化。
直接法[2]：
通用、穩(wěn)定；通過前后處理，能夠保證計(jì)算的收斂性與精度；
對原矩陣的編輯多（分解、排序、縮放等）；
內(nèi)存需求大，求解速度慢，算法復(fù)雜度更高；
并行度有限。
其中迭代法的種類很多，可以分為定常（Stationary）迭代法與非定常迭代法[3]。
展開 技術(shù)分享︱基于非結(jié)構(gòu)網(wǎng)格的仿真——太湖之光上的巨大挑戰(zhàn)
稀疏矩陣向量乘：
用COO格式稀疏矩陣向量乘進(jìn)行測試。隨著矩陣元素維度變化，UNAT最高加速比和單核組Flops分別達(dá)到近27 倍和3.5 GFlops（0.46%理論性能）。與手工優(yōu)化的代碼相比，UNAT API實(shí)現(xiàn)了其大約70%的性能。考慮到兼容性和通用性帶來的額外開銷，其效果是相當(dāng)可觀的。
COO稀疏矩陣向量乘加速比及浮點(diǎn)性能
某CFD求解器加速：
一個(gè)開發(fā)者僅用1周時(shí)間，采用UNAT對其中主要6個(gè)熱點(diǎn)核心進(jìn)行了加速開發(fā)。熱點(diǎn)核心眾核加速比平均在19倍左右，整體加速也達(dá)到了10倍。
某CFD求解器優(yōu)化加速
OpenFOAM風(fēng)資源評估應(yīng)用：
利用UNAT加速工具，我們對OpenFOAM某風(fēng)資源評估應(yīng)用大部分熱點(diǎn)進(jìn)行了加速，眾核加速比在8-15倍之間。由于UNAT卓越貢獻(xiàn)，最終該應(yīng)用在太湖之光上取得了整體4倍加速。
OpenFOAM風(fēng)資源評估應(yīng)用優(yōu)化加速
燃燒仿真應(yīng)用：
利用UNAT加速工具，單個(gè)開發(fā)者 2周內(nèi)完成了15萬行某燃燒仿真程序優(yōu)化，使得該程序最終在太湖之光上取得了整體5.4倍加速。由于性能提升，該軟件可使用10億級(jí)網(wǎng)格，在合理的時(shí)間內(nèi)，進(jìn)行航空發(fā)動(dòng)機(jī)全環(huán)燃燒室高保真數(shù)值模擬。
?
燃燒仿真應(yīng)用優(yōu)化加速
參考文獻(xiàn)：
[1] Vincent, P., Witherden, F., Vermeire, B., Park, J. S., & Iyer, A. (2016). Towards Green Aviation with Python at Petascale.
展開 基于非結(jié)構(gòu)網(wǎng)格的仿真——太湖之光上的巨大挑戰(zhàn)
稀疏矩陣向量乘：用COO格式稀疏矩陣向量乘進(jìn)行測試。隨著矩陣元素維度變化，UNAT最高加速比和單核組Flops分別達(dá)到近27 倍和3.5 GFlops（0.46%理論性能）。與手工優(yōu)化的代碼相比，UNAT API實(shí)現(xiàn)了其大約70%的性能。考慮到兼容性和通用性帶來的額外開銷，其效果是相當(dāng)可觀的。
COO稀疏矩陣向量乘加速比及浮點(diǎn)性能
某CFD求解器加速：一個(gè)開發(fā)者僅用1周時(shí)間，采用UNAT對其中主要6個(gè)熱點(diǎn)核心進(jìn)行了加速開發(fā)。熱點(diǎn)核心眾核加速比平均在19倍左右，整體加速也達(dá)到了10倍。
某CFD求解器優(yōu)化加速
OpenFOAM風(fēng)資源評估應(yīng)用：利用UNAT加速工具，我們對OpenFOAM某風(fēng)資源評估應(yīng)用大部分熱點(diǎn)進(jìn)行了加速，眾核加速比在8-15倍之間。由于UNAT卓越貢獻(xiàn)，最終該應(yīng)用在太湖之光上取得了整體4倍加速。
OpenFOAM風(fēng)資源評估應(yīng)用優(yōu)化加速
燃燒仿真應(yīng)用：利用UNAT加速工具，單個(gè)開發(fā)者 2周內(nèi)完成了15萬行某燃燒仿真程序優(yōu)化，使得該程序最終在太湖之光上取得了整體5.4倍加速。由于性能提升，該軟件可使用10億級(jí)網(wǎng)格，在合理的時(shí)間內(nèi)，進(jìn)行航空發(fā)動(dòng)機(jī)全環(huán)燃燒室高保真數(shù)值模擬。
燃燒仿真應(yīng)用優(yōu)化加速
參考文獻(xiàn)：
[1] Vincent, P., Witherden, F., Vermeire, B., Park, J. S., & Iyer, A. (2016).
展開 技術(shù)分享︱大型稀疏線性方程組求解技術(shù)——工業(yè)仿真的底層核心
一、背景
在工業(yè)仿真領(lǐng)域，對各種現(xiàn)實(shí)世界的問題進(jìn)行數(shù)值模擬時(shí)，如流體動(dòng)力學(xué)分析、電磁場仿真、結(jié)構(gòu)力學(xué)應(yīng)力應(yīng)變分析等，其控制方程通常是偏微分方程組，在經(jīng)過不同方法的隱式離散之后最終都可轉(zhuǎn)化為大型稀疏線性方程組。隨著人們對計(jì)算精度要求的不斷提高，方程組的階數(shù)也從上千階、幾十萬階提高到百萬、千萬階甚至更高，所需的計(jì)算量以及存儲(chǔ)需求也隨之迅速膨脹。根據(jù)一般經(jīng)驗(yàn)，方程組求解時(shí)間會(huì)占總計(jì)算時(shí)間的70%以上，往往是整個(gè)計(jì)算過程中的性能瓶頸。如果說求解器是工業(yè)CAE軟件的核心模塊，那么大型稀疏線性方程組的求解技術(shù)將毫無疑問是底層求解器的核心。
NASA翼型網(wǎng)格經(jīng)過離散得到的稀疏矩陣（素材來源于網(wǎng)絡(luò)）
二、方法
眾所周知，稀疏線性方程組的求解方法可以分為直接法和迭代法 ，兩類方法各有優(yōu)劣，特點(diǎn)比較如下：
迭代法[1]：
對于不同類型稀疏矩陣表現(xiàn)差異較大，存在收斂性與收斂速度問題，催生了許多預(yù)處理技術(shù)（Preconditioners）；
對原矩陣的編輯很少，SpMV(Sparse matrix-vector multiplication)是其核心運(yùn)算；
內(nèi)存需求小，求解速度較快，算法復(fù)雜度低；
較易實(shí)現(xiàn)并行化。
直接法[2]：
通用、穩(wěn)定；通過前后處理，能夠保證計(jì)算的收斂性與精度；
對原矩陣的編輯多（分解、排序、縮放等）；
內(nèi)存需求大，求解速度慢，算法復(fù)雜度更高；
并行度有限。
其中迭代法的種類很多，可以分為定常（Stationary）迭代法與非定常迭代法[3]。
展開 
多體動(dòng)力學(xué)歷史的重要一頁
他也使用了稀疏矩陣方法，保證約束方程在位移、速度、加速度層面都是滿足的。稀疏矩陣配合廣義坐標(biāo)分解，并不需要對雅克比矩陣求逆，奠定了另外一款軟件DADS（Dynamic Analysis and Design System）。Wehage在愛荷華大學(xué)當(dāng)助力教授，后來去了TARDEC，1995年，又加入了卡特彼勒。隨著1979年DADS的開發(fā)，他于1985年指導(dǎo)開發(fā)了位于德國IABG的NUSTAR MBD程序，在1991年指導(dǎo)開發(fā)了TRADEC的SOVAS MBD程序，1997年指導(dǎo)了卡特彼勒的TomSim MBD程序。
Roger Wehage在伊利諾斯大學(xué)的講話2016.5.27
References
[1] Orlandea, N., Chace, M. A., and Calahan, D. A., 1977, “A Sparsity Oriented Approach to Dynamic Analysis and Design of Mechanical Systems,” ASME J. Eng. Ind., 99, pp. 773-784.
[2] Uicker, J. J., Ravani, B., and Sheth, P. N., 2013, Matrix Method in the Design Analysis of Mechanisms and Multibody Systems, Cambridge University Press, New York, NY.
[3] Sheth, P. N., and Uicker, J.
展開 Julia：高效易用的數(shù)值計(jì)算/優(yōu)化編程語言
1//1 : 0//1 for i in 1:L];
P = zeros(t_max,L)
dl = 0.5 * vec(ones(L-1,1))
dl[1] = 0
d = vec(zeros(L,1))
d[1] = 1
d[L] = 0.5
du = 0.5 * vec(ones(L-1,1))
Trans = Tridiagonal(dl,d,du)
P[1,:] = Trans * vec(p)
for i = 2:t_max
P[i,:] = Trans * vec(P[i-1,:])
endend
注意到在稀疏矩陣那里我們利用了三對角矩陣(Tridiagonal)的數(shù)據(jù)結(jié)構(gòu)，而benchmarkTools都做了幾百次的sample然后輸出該代碼塊運(yùn)算時(shí)間的各種統(tǒng)計(jì)量。我們于是知道稀疏矩陣快于array comprehension快于naive的稠密矩陣。在數(shù)值計(jì)算中，使用最有效的數(shù)據(jù)結(jié)構(gòu)和代碼命令永遠(yuǎn)是很重要的。至于如果想驗(yàn)證（sanity check）這三個(gè)代碼塊是否給出了一致的輸出結(jié)果，可以使用@test_approx_eq命令，這里不多贅述。
例子2：稀疏矩陣，自定義數(shù)據(jù)結(jié)構(gòu)，分布式數(shù)據(jù)結(jié)構(gòu)（初步）
我們在上例中看到了一類特殊的稀疏矩陣——三對角矩陣類型的應(yīng)用。這里我們指出更一般的稀疏矩陣類型在Julia中名為sparse()。
展開 【JY】土木工程振型求解之蘭索斯法(Lanczos法)
而且計(jì)算速度非常高效，因此針對稀疏矩陣高效的特征值、特征向量的計(jì)算方法可以采用可存儲(chǔ)量個(gè)大，計(jì)算速度較快的Lanczos迭代法。
提取ANSYS中的剛度矩陣和質(zhì)量矩陣到MATLAB中進(jìn)行二次開發(fā) ￥88
對稱矩陣只存儲(chǔ)下三角元素，如結(jié)構(gòu)剛度矩陣為對稱矩陣，Harwell-Boeing格式則僅記錄下三角元素。
通過上述解釋我們可以知道該剛度矩陣有1027個(gè)列指針，有7978個(gè)行索引，有1026個(gè)右邊項(xiàng)，矩陣的維度為1026×1026。
要將上述文件讀取到MATLAB中并轉(zhuǎn)化為常規(guī)的稀疏矩陣，需要了解Harwell-Boeing壓縮儲(chǔ)存格式及重構(gòu)方法，這里面采用的是CSC方式，以列指針、行索引的方式去重構(gòu)稀疏矩陣。
