小さなsparsematricesから固有3 sparsematrixを組み立てる

Question

私は結合されたマルチフィジックスシステムのジャコビアを組み立てています。ヤコビアンは、各システムの対角線上のブロック行列と、結合のための対角ブロックとからなる。 アセンブラを使って分離してブロックし、射影行列で総和して完全なヤコビ行列を得ることが最善の方法です。 擬似コード（J [i]は対角要素、C [ij]は結合）、Pは完全行列への射影です。

// diagonal blocks 
J.setZero(); 
for(int i=0;i<N;++i){ 
    J+=P[i]J[i]P[i].transpose() 
} 
// off diagonal elements 
for(int i=0;i<N;++i){ 
    for(int j=i+1;j<N;++j){ 
     J+=P[i]C[ij]P[j].transpose() 
     J+=P[j]C[ji]P[i].transpose() 
    } 
} 

これは、全体の20％程度のパフォーマンスを必要とします。これは、組み立てには多すぎます。システムが非常に非線形であるため、毎回ジャコビアを再計算する必要があります。 Valgrindは、リソース消費方法がEigen::internal::assign_sparse_to_sparseで、このメソッドではEigen::SparseMatrix<>::InsertBackByOuterInnerを呼び出していることを示します。

このようなマトリックスをより効率的に組み立てる方法はありますか？

（私もプログラムの開発をコンパイルするために代わりにP J * J.transposeのP *（J P.transposeを（））（）を使用しなければならなかった、間違って何かがすでに存在してもよい）

PS：余分なマトリックスにP.transposeを格納することによって、私は少し良いパフォーマンスを得るが、プログラムの開発の15％のため、まだ加算アカウントが

Answer 1

あなたのコードは次のようになります。NDEBUGと最適化が

編集をオンにしていますインプレースで作業することではるかに高速です。まず、（まだ行っていない場合）、最終的なマトリックスと予備スペースに列あたりの非ゼロの数を見積もる：

int nnz_per_col = ...; 
J.reserve(VectorXi::Constant(n_cols, nnz_per_col)); 

列ごとNNZの数は非常に不均一であるならば、あなたも計算することができますそれは列ごとに：fooは、インデックスの適切なマッピングを実装

for each block B[k] 
    for each elements i,j 
    J.coeffRef(foo(i),foo(j)) += B[k](i,j) 

：

VectorXi nnz_per_col(n_cols); 
for each j 
    nnz_per_col(j) = ...; 
J.reserve(nnz_per_col); 

手動要素を挿入します。

そして次の反復のため、予約する必要はありませんが、構造を維持しながら、あなたがゼロに係数値を設定する必要があります。

J.coeffs().setZero();