FEMの疎なマトリックスを効率的に組み立てる方法

Question

すべてを扱う、 私の質問を読んでいただきありがとうございます。 私はEigen3.3.4（http://eigen.tuxfamily.org/index.php?title=Main_Page）を使っていくつかのFEMコードを書いています。

私はEigen3.3.4のドキュメントを読んでおり、このウェブサイト（http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html）には と書かれていますが、追加コピーを避けて高性能を得るにはRef<MatrixBase>を使用してください。疎行列のための私のFEMコードでそう

は、部品を組み立て、の関数であるとしましょう：

Uは変位を表す

FormFE(const Ref<VectorXd> &U,const Ref<VectorXd> &V, 
Ref<SparseMatrix<double> > AMATRIX,Ref<VectorXd> RHS) 

、Vは速度項を表します。 AMATRIXは私の疎な行列、RHSは残差項です。

それから私は、だから私は試した（私は（私は）初期化するための値をゼロに設定し、すべての非ゼロ要素とその値が含まれているtripletListを持っている） 組み立て最初の前に私のマトリクス状を初期化しよう：

をAMATRIX.setFromTriplets(ZeroTripList.begin(),ZeroTripList.end());

しかし、私はエラーを持っている：

class Eigen::Ref<Eigen::SparseMatrix<double, 0, int> >’ has no member named ‘setFromTriplets 

だから、どのように私はこの問題を解決するだろうか？

私の解決策の一つは、使用している：

FormFE(const Ref<VectorXd> &U,const Ref<VectorXd> &V, 
SparseMatrix<double> &AMATRIX,Ref<VectorXd> RHS) 

これは非常に正常に動作しているが、私はそれが効率的であるかどうかわかりません。 私はCPPで非常によくないよ：P実は

、私の質問は以下のとおりです。

1. が効率的に（特にFEM計算のための）固有の使用方法、私はそれぞれにほとんどどこでも固有のVectorXdとMatrixXd使用します私のFEM関連の機能。
2. SparseMatrixを効率的に組み立てるには？
3. FEMアセンブル用にOpenMP並列化を行うことはできますか？
4. C++ベースのFEMコーディング（ライブラリの推奨または有用なアイデア）の参考にしてください。

ありがとうございます。 よろしくお願いいたします。

Answer 1

はい、ここをクリックしてSparseMatrix<double> &を渡してください。 Ref<SparseMatrix>の目的は、Map<SparseMatrix> ...

使用setFromTripletsもOKパフォーマンスを得るようにしてくださいするために行うには正しいことで、サブスパース行列のように、SparseMatrixにressemble組み立てオブジェクトを渡すです。 mat.insert(i,j) = val;を使用して要素を直接挿入すると、適切に実行すると（つまり、予約と適切な挿入オーダーの適切な呼び出し）、最大2倍高速になる可能性があります。しかし、それが間違っている場合、それはx100倍遅くなる可能性もあります。文書を見てください。SparseMatrix::insertで、それは、OpenMPを使用して行列を充填することも可能であるが、これはもっと慎重かつ厳格であることを要求し、ここで、このための典型的なパターンです：もちろん

int n_cols = ??, n_rows = ??; 
std::vector<int> nnz_per_col(n_cols); 
// set each nnz_per_col[j] to the exact number 
// of non-zero entries in the j-th column (or more, but NOT less) 
SparseMatrix<double> mat(n_rows, n_cols); 
#pragma omp parallel for 
for(int j=0; j<cols; ++j) { 
    for each non zero entry i in the j-th column { 
    // preferably with increasing i 
    double val_i_j = ...; 
    mat.insert(i,j) = val_i_j; 
    } 
} 

は、あなたもあれば、行単位で作業することができますそれはあなたのために簡単です。この場合は、SparseMatrix<double,RowMajor>を使用してください。もちろん、このパターンを調整して列や行などのブロックを処理することもできます。

アセンブリの場合、高密度マトリックス/ベクターで作業する必要がある場合は、固定サイズ。 MatrixXd/VectorXdを使用する代わりに、静的に割り当てられるタイプのMatrix<double,N,M>とMatrix<double,N,1>を使用することをお勧めします。これにより、多数のメモリ割り当て/割り当て解除が防止されます。

最後に、最も重要な勧告：あなたは、パフォーマンスを心配している場合、あなたのコードを最適化した場合の時間と労力を調査する前に、あなたのコードをプロファイリングすることを忘れないでください。また、コンパイラの最適化をONにして常にベンチ/プロファイル。