基于超節點LDL分解的大規模結構計算

賴智超　羅曉群　張其林

摘要：采用列壓縮稀疏（Compressed Sparse Column， CSC）矩陣存儲策略對矩陣LDL分解前進行填充元優化排序；基于消去樹進行LDL符號分解，使之獨立于數值分解，避免多余的內存消耗，減少不必要的數值運算.利用矩陣非零元的分布特性分析并實現超節點LDL分解算法，將稀疏矩陣的分解運算變為一系列稠密矩陣運算，并使用優化的BLAS函數庫加速分解.測試表明：算法在成倍地提高計算速度的同時進一步降低內存消耗，適用于大規模的結構計算.

關鍵詞：稀疏矩陣； LDL分解； 消去樹； 符號分解； 超節點

中圖分類號： TP301

文獻標志碼：B

現代結構工程的計算經有限元法處理后往往歸結為一個規模龐大的線性代數方程組的求解.大規模的矩陣存儲和快速的求解是制約大型結構計算的關鍵.

稀疏矩陣是指矩陣中大多數數值為零的矩陣.對于較大規模的稀疏矩陣，非零元占總元素個數比例非常小.結構分析中得到的矩陣就是一個具有很好特性的對稱正定（Symmetric Positive Definite，SPD）稀疏矩陣.本文的關注點是結構工程中的SPD稀疏矩陣的直接法求解計算.目前，已經有不少穩定高效的求解庫，如csparse，taucs，cholmod和superLU等，它們主要是使用cholesky分解和LU分解將稀疏矩陣A分解為LLT或

LU的形式.由于分解后的L中非零元相對于A會增多，所以稱之為“填充元”.[1]因此，一般將稀疏矩陣的分解分為2部分：計算填充元的位置，稱之為“符號分解”；進行實際的數值分解計算，稱之為“數值分解”.

第一個通用快速SPD稀疏矩陣求解器實現庫是SPARSPAK庫[2]，具有里程碑的意義.隨后，稀疏矩陣快速直接解法的算法研究與圖論算法緊密結合，如消去樹[3]性質研究、最小度算法[4]和圖剖分算法[5]等，同時算法的高效性還充分考慮到計算機硬件的內存布局和cache命中率等，因此稀疏矩陣算法得到迅速發展，已經是一個集工程結構、計算數學和計算機應用科學于一體的綜合性學術研究領域.

本文采用列壓縮稀疏（Compressed Sparse Column， CSC）矩陣存儲[6]策略，對矩陣LDL分解前進行優化排序減少填充元，基于消去樹進行LDL符號分解，使之獨立于LDL數值分解.由于非零元相同的稀疏矩陣符號分解結果相同，因此可以重復使用，避免多余的內存消耗，減少不必要的數值運算.進一步充分利用矩陣非零元的分布特性，分析并實現超節點LDL分解算法，使稀疏矩陣的分解運算轉變成一系列稠密矩陣運算，并使用優化的BLAS函數庫加速分解，在成倍地提高計算速度的同時進一步降低內存消耗，適用于大規模結構計算問題.

1 稀疏三角系統Lx=b的求解

本節討論稀疏三角系統Lx=b的求解，其中：

L為稀疏下三角矩陣，b為稀疏向量.該算法是矩陣分解算法內部的基本運算子單元，擁有高效的稀疏右端項求解算法，直接影響矩陣分解算法的效率，因而具有重要意義.

1.1 求解原理

因為CSC格式存儲矩陣的特點是能夠快速地讀取一列中的元素，因此算法以列的方式對矩陣操作，提高求解速度.式（1）～（3）可以滿足要求.

由于右端項b是稀疏的，因此解x大多也是稀疏的.所以，首先要快速地判斷求解結果x中非零元的位置模式，之后采用式（1）～（3）求解x.

1.2 非零元位置模式求解算法

該問題可以抽象為一個圖論數學模型——有向圖的可達性遍歷[7]，其中，

有向圖G構建方式為

摘要：采用列壓縮稀疏（Compressed Sparse Column， CSC）矩陣存儲策略對矩陣LDL分解前進行填充元優化排序；基于消去樹進行LDL符號分解，使之獨立于數值分解，避免多余的內存消耗，減少不必要的數值運算.利用矩陣非零元的分布特性分析并實現超節點LDL分解算法，將稀疏矩陣的分解運算變為一系列稠密矩陣運算，并使用優化的BLAS函數庫加速分解.測試表明：算法在成倍地提高計算速度的同時進一步降低內存消耗，適用于大規模的結構計算.

關鍵詞：稀疏矩陣； LDL分解； 消去樹； 符號分解； 超節點

中圖分類號： TP301

文獻標志碼：B

現代結構工程的計算經有限元法處理后往往歸結為一個規模龐大的線性代數方程組的求解.大規模的矩陣存儲和快速的求解是制約大型結構計算的關鍵.

稀疏矩陣是指矩陣中大多數數值為零的矩陣.對于較大規模的稀疏矩陣，非零元占總元素個數比例非常小.結構分析中得到的矩陣就是一個具有很好特性的對稱正定（Symmetric Positive Definite，SPD）稀疏矩陣.本文的關注點是結構工程中的SPD稀疏矩陣的直接法求解計算.目前，已經有不少穩定高效的求解庫，如csparse，taucs，cholmod和superLU等，它們主要是使用cholesky分解和LU分解將稀疏矩陣A分解為LLT或

LU的形式.由于分解后的L中非零元相對于A會增多，所以稱之為“填充元”.[1]因此，一般將稀疏矩陣的分解分為2部分：計算填充元的位置，稱之為“符號分解”；進行實際的數值分解計算，稱之為“數值分解”.

第一個通用快速SPD稀疏矩陣求解器實現庫是SPARSPAK庫[2]，具有里程碑的意義.隨后，稀疏矩陣快速直接解法的算法研究與圖論算法緊密結合，如消去樹[3]性質研究、最小度算法[4]和圖剖分算法[5]等，同時算法的高效性還充分考慮到計算機硬件的內存布局和cache命中率等，因此稀疏矩陣算法得到迅速發展，已經是一個集工程結構、計算數學和計算機應用科學于一體的綜合性學術研究領域.

本文采用列壓縮稀疏（Compressed Sparse Column， CSC）矩陣存儲[6]策略，對矩陣LDL分解前進行優化排序減少填充元，基于消去樹進行LDL符號分解，使之獨立于LDL數值分解.由于非零元相同的稀疏矩陣符號分解結果相同，因此可以重復使用，避免多余的內存消耗，減少不必要的數值運算.進一步充分利用矩陣非零元的分布特性，分析并實現超節點LDL分解算法，使稀疏矩陣的分解運算轉變成一系列稠密矩陣運算，并使用優化的BLAS函數庫加速分解，在成倍地提高計算速度的同時進一步降低內存消耗，適用于大規模結構計算問題.

1 稀疏三角系統Lx=b的求解

本節討論稀疏三角系統Lx=b的求解，其中：

L為稀疏下三角矩陣，b為稀疏向量.該算法是矩陣分解算法內部的基本運算子單元，擁有高效的稀疏右端項求解算法，直接影響矩陣分解算法的效率，因而具有重要意義.

1.1 求解原理

因為CSC格式存儲矩陣的特點是能夠快速地讀取一列中的元素，因此算法以列的方式對矩陣操作，提高求解速度.式（1）～（3）可以滿足要求.

由于右端項b是稀疏的，因此解x大多也是稀疏的.所以，首先要快速地判斷求解結果x中非零元的位置模式，之后采用式（1）～（3）求解x.

1.2 非零元位置模式求解算法

該問題可以抽象為一個圖論數學模型——有向圖的可達性遍歷[7]，其中，

有向圖G構建方式為

摘要：采用列壓縮稀疏（Compressed Sparse Column， CSC）矩陣存儲策略對矩陣LDL分解前進行填充元優化排序；基于消去樹進行LDL符號分解，使之獨立于數值分解，避免多余的內存消耗，減少不必要的數值運算.利用矩陣非零元的分布特性分析并實現超節點LDL分解算法，將稀疏矩陣的分解運算變為一系列稠密矩陣運算，并使用優化的BLAS函數庫加速分解.測試表明：算法在成倍地提高計算速度的同時進一步降低內存消耗，適用于大規模的結構計算.

關鍵詞：稀疏矩陣； LDL分解； 消去樹； 符號分解； 超節點

中圖分類號： TP301

文獻標志碼：B

現代結構工程的計算經有限元法處理后往往歸結為一個規模龐大的線性代數方程組的求解.大規模的矩陣存儲和快速的求解是制約大型結構計算的關鍵.

稀疏矩陣是指矩陣中大多數數值為零的矩陣.對于較大規模的稀疏矩陣，非零元占總元素個數比例非常小.結構分析中得到的矩陣就是一個具有很好特性的對稱正定（Symmetric Positive Definite，SPD）稀疏矩陣.本文的關注點是結構工程中的SPD稀疏矩陣的直接法求解計算.目前，已經有不少穩定高效的求解庫，如csparse，taucs，cholmod和superLU等，它們主要是使用cholesky分解和LU分解將稀疏矩陣A分解為LLT或

LU的形式.由于分解后的L中非零元相對于A會增多，所以稱之為“填充元”.[1]因此，一般將稀疏矩陣的分解分為2部分：計算填充元的位置，稱之為“符號分解”；進行實際的數值分解計算，稱之為“數值分解”.

第一個通用快速SPD稀疏矩陣求解器實現庫是SPARSPAK庫[2]，具有里程碑的意義.隨后，稀疏矩陣快速直接解法的算法研究與圖論算法緊密結合，如消去樹[3]性質研究、最小度算法[4]和圖剖分算法[5]等，同時算法的高效性還充分考慮到計算機硬件的內存布局和cache命中率等，因此稀疏矩陣算法得到迅速發展，已經是一個集工程結構、計算數學和計算機應用科學于一體的綜合性學術研究領域.

本文采用列壓縮稀疏（Compressed Sparse Column， CSC）矩陣存儲[6]策略，對矩陣LDL分解前進行優化排序減少填充元，基于消去樹進行LDL符號分解，使之獨立于LDL數值分解.由于非零元相同的稀疏矩陣符號分解結果相同，因此可以重復使用，避免多余的內存消耗，減少不必要的數值運算.進一步充分利用矩陣非零元的分布特性，分析并實現超節點LDL分解算法，使稀疏矩陣的分解運算轉變成一系列稠密矩陣運算，并使用優化的BLAS函數庫加速分解，在成倍地提高計算速度的同時進一步降低內存消耗，適用于大規模結構計算問題.

1 稀疏三角系統Lx=b的求解

本節討論稀疏三角系統Lx=b的求解，其中：

L為稀疏下三角矩陣，b為稀疏向量.該算法是矩陣分解算法內部的基本運算子單元，擁有高效的稀疏右端項求解算法，直接影響矩陣分解算法的效率，因而具有重要意義.

1.1 求解原理

因為CSC格式存儲矩陣的特點是能夠快速地讀取一列中的元素，因此算法以列的方式對矩陣操作，提高求解速度.式（1）～（3）可以滿足要求.

由于右端項b是稀疏的，因此解x大多也是稀疏的.所以，首先要快速地判斷求解結果x中非零元的位置模式，之后采用式（1）～（3）求解x.

1.2 非零元位置模式求解算法

該問題可以抽象為一個圖論數學模型——有向圖的可達性遍歷[7]，其中，

有向圖G構建方式為