[論文レビュー] An efficient way to perform the assembly of finite element matrices in Matlab and Octave

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

Manufacturing Process and Optimization参考文献 8被引用数 17

ひとこと要約

本稿では、MATLABおよびOctave向けに、完全にベクトル化されループを含まないアルゴリズムを提示する。この手法は、2次元のP1ラグランジュ有限要素および線形弾性問題における質量行列と剛性行列の効率的なアセンブリを実現し、FreeFEM++などの専用ソフトウェアと同等の性能を達成している。低レベルのコーディングを回避し、最適化されたスパース行列処理とアレイベースのインデキシングを活用することで、標準的および先行のベクトル化手法を凌駕する。

ABSTRACT

We describe different optimization techniques to perform the assembly of finite element matrices in Matlab and Octave, from the standard approach to recent vectorized ones, without any low level language used. We finally obtain a simple and efficient vectorized algorithm able to compete in performance with dedicated software such as FreeFEM++. The principle of this assembly algorithm is general, we present it for different matrices in the P1 finite elements case and in linear elasticity. We present numerical results which illustrate the computational costs of the different approaches

研究の動機と目的

• MATLABおよびOctaveにおける有限要素行列アセンブリの非効率性を是正する。これは、ループベースの実装に起因する。
• MATLAB R2012a/bで発生するスパース行列の再割り当てに起因するパフォーマンス劣化を克服する。
• ループを排除し、アレイ操作を用いることで、専用FEMソフトウェアと同等またはそれ以上の性能を達成するベクトル化されたアセンブリアルゴリズムを開発する。
• 本手法の一般性と拡張可能性を示し、Pk要素および3次元問題への応用を含め、線形弾性問題への適用を検証する。
• MATLABおよびOctaveと完全に互換性があり、効率的かつ容易にカスタマイズ可能なオープンソース実装を提供する。

提案手法

• 標準のループベースのアセンブリを、グローバル行列要素のための事前計算済みインデックス配列（Ig, Jg, Kg）を用いた完全なベクトル化アプローチに置き換える。
• 明示的なループを回避し、すべての要素に対して同時に要素行列（例：質量行列や剛性行列）を行列アレイ演算で計算する。
• MATLAB R2012a/bにおけるパフォーマンス劣化を防ぐために、sparse(nq,nq)の代わりにspalloc(nq,nq,9*nme)を用いてスパース行列の初期化を最適化する。
• アフィン変換と数値積分を、要素固有の関数（例：ElemStiffMatP1）を用いて実行し、局所行列を計算する。
• すべての要素寄与を一度の関数呼び出しで効率的にアセンブルするために、sparse(Ig, Jg, Kg, nq, nq)を用いてグローバル行列を構築する。
• 言語固有の構文や低レベルな構文を避けることで、MATLABおよびGNU Octaveの両方と互換性を保つ。

実験結果

リサーチクエスチョン

• RQ1 MATLAB/Octaveにおける完全なベクトル化・ループフリーなアルゴリズムは、FreeFEM++などの専用FEMソフトウェアと同等のパフォーマンスを達成できるか？
• RQ2 なぜMATLAB R2012a/bにおける標準的なスパース行列初期化は、以前のバージョンと比較して顕著なパフォーマンス劣化を引き起こすのか？
• RQ3 数値積分の公式や低レベル言語を用いずに、どのようにして有限要素行列アセンブリにベクトル化を適用できるか？
• RQ4 提案手法は、他の有限要素タイプ（例：Pk、3次元）および物理的問題（例：線形弾性）へどの程度一般化可能か？
• RQ5 効率的なスパース行列初期化は、全体のアセンブリパフォーマンスにどのような役割を果たし、どのように最適化できるか？

主な発見

• 提案されたベクトル化アルゴリズム（OptV2）は、2次元P1有限要素においてFreeFEM++と同等のパフォーマンスを達成し、同じハードウェア上での計算時間はFreeFEM++の1.5倍以内に収まる。
• MATLAB R2012a/bにおけるスパース行列初期化のパフォーマンス劣化は、著者らによって確認され、MathWorksに報告されており、spallocがその解決策であると特定された。
• sparseの代わりにspallocを使用することで、1600自由度では最大44％の計算時間短縮が達成され、40000自由度では最大2.95倍の高速化が実現された（R2011aと比較）。
• ベクトル化アプローチ（OptV2）は、古典的なループベース手法（OptV0）および先行のベクトル化手法（例：[2, 3, 12, 17]）を上回り、速度とシンプルさの両面で優れている。
• 本手法は線形弾性問題への応用を通じて、一般性と再利用可能性が明確に示された。
• 実装全体がMATLABおよびGNU Octaveと完全に互換性があり、広範な利用可能性と再現可能性を実現している。