[論文レビュー] An efficient way to assemble finite element matrices in vector languages

François Cuvelier, Caroline Japhet+1 more|arXiv (Cornell University)|Jan 14, 2014

Matrix Theory and Algorithms参考文献 14被引用数 4

ひとこと要約

本稿では、MATLAB、Octave、Python、Juliaなどのハイレベルなベクトル言語を用いて、任意の次元における有限要素行列のアセンブリに特化した高度に最適化されたベクトル化アルゴリズムを提示する。ループを排除するための高度なスパース行列演算と対称性の活用により、提案されたOptVSアルゴリズムはスカラ型の場合はC言語の2.9倍以内、ベクトル型の場合は4.1倍以内の性能を達成し、適切に最適化された場合、ベクトル言語がコンパイル済み言語と同等の性能を発揮できることを示している。

ABSTRACT

Efficient Matlab codes in 2D and 3D have been proposed recently to assemble finite element matrices. In this paper we present simple, compact and efficient vectorized algorithms, which are variants of these codes, in arbitrary dimension, without the use of any lower level language. They can be easily implemented in many vector languages (e.g. Matlab, Octave, Python, Scilab, R, Julia, C++ with STL,...). The principle of these techniques is general, we present it for the assembly of several finite element matrices in arbitrary dimension, in the P1 finite element case. We also provide an extension of the algorithms to the case of a system of PDE's. Then we give an extension to piecewise polynomials of higher order. We compare numerically the performance of these algorithms in Matlab, Octave and Python, with that in FreeFEM++ and in a compiled language such as C. Examples show that, unlike what is commonly believed, the performance is not radically worse than that of C : in the best/worst cases, selected vector languages are respectively 2.3/3.5 and 2.9/4.1 times slower than C in the scalar and vector cases. We also present numerical results which illustrate the computational costs of these algorithms compared to standard algorithms and to other recent ones.

研究の動機と目的

• ハイレベルなベクトル言語（MATLAB、Octave、Python、Juliaなど）を用いて、任意の次元における有限要素行列のアセンブリに適した効率的でポータブルかつ読みやすいベクトル化アルゴリズムの開発。
• 科学計算環境におけるループベースの行列アセンブリの性能ボトルネックの克服。
• MATLAB、Octave、Pythonなどの言語におけるベクトル実装が、コンパイル済みCコードと同等の性能を達成できることの実証。
• PDEの連立系および高次有限要素への応用を含め、広範な適用性を確保するためのアプローチの拡張。

提案手法

• 著者らは、ネストされたループを避けるためにスパース行列演算に基づくベクトル化アセンブリ戦略を導入し、要素ごとの配列演算と行列乗算に置き換える。
• コア技術は、グローバルインデックスマッピングとキロンカー型演算を用いて、要素寄与を一括して計算することで、明示的な要素ループを回避する。
• OptVSアルゴリズムは、行列の対称性を活用して重複計算を削減し、効率性とメモリ使用効率を向上させる。
• 本手法は任意の次元に一般化され、係数が変化するPDEの連立系に対しても拡張可能である。
• 本手法は、MATLABやPythonに実装された標準的なスパース行列データ構造（例：CSR/CSC）を用いて実装される。
• アルゴリズムは、2次元および3次元の剛性行列や弾性剛性行列を含む、複数の言語および問題タイプで検証されている。

実験結果

リサーチクエスチョン

• RQ1 MATLAB、Octave、Pythonなどのハイレベル言語におけるベクトル化アルゴリズムは、有限要素行列アセンブリにおいて、コンパイル済みCコードと同等の性能を達成できるか？
• RQ2 問題サイズおよび次元数の増加に伴い、ループベースとベクトル化アセンブリ手法の計算複雑性と性能はどのようにスケーリングするか？
• RQ3 有限要素法における行列アセンブリ最適化のために、対称性とスパース行列演算をどの程度活用できるか？
• RQ4 MATLAB、Octave、Pythonなどの異なるベクトル言語において、同一の最適化アルゴリズムを実装した場合の性能比較はどのようになるか？
• RQ5 提案されたベクトル化戦略は、PDEの連立系および高次有限要素へ一般化可能か？

主な発見

• OptVSアルゴリズムは、O(n_dof)の計算複雑性を達成しており、行列アセンブリにおいて漸近的に最適である。
• 最良のケースでは、スカラ問題ではC言語の2.3倍以内、ベクトル問題では2.9倍以内の速度で実行される。
• 2次元剛性行列において、ベースライン実装はOctaveでOptVSバージョン比最大550倍遅い。
• OptV1バージョンは、異なる言語間でOptVS比40〜550倍遅く、ベクトル化による性能向上の顕著さが示された。
• スパース行列演算を適切に最適化することで、ベクトル言語とC言語の性能差は一般的に予想されているほどではない。
• 本手法は3次元および高次要素に対しても効率的にスケーリングでき、高次P1有限要素への適用にはわずかな修正のみが必要である。