Eigenスパース行列、既に最適化されていた

C++で置き換えてみたのですが、圧倒的に元のEigenの方が速かったです。それで、Eigen核心の部分のコードをDisassembleして見たのですが、次のように僅か5命令で、疎行列とベクタをロードし積和演算を行うことが達成されていました。見事なコードです。 ソースはC++であり、どのように書いたらこのようなコードが生成されるのか、謎ですが、真似て書いても、中々、下のようなコードは生成されませんでした。(少なくとも二重ベクタループでは駄目でした。)

C++のvectorを使った二重ループの速度比較 - Qiita

00007FF6340034A0  movsxd      rax,dword ptr [r11+rcx*4]  
00007FF6340034A4  movsd       xmm0,mmword ptr [r9+rax*8]  
00007FF6340034AA  mulsd       xmm0,mmword ptr [r10+rcx*8]  
00007FF6340034B0  inc         rcx  
00007FF6340034B3  addsd       xmm1,xmm0  
00007FF6340034B7  cmp         rcx,rdx  
00007FF634003ま

また、マルチスレッド化は、OpenMPを用いた記述になっており、推定要素数に応じてマルチスレッドが駆動されます。実際、今回のインスタンスについては、ほぼ、スレッド数倍速になっており、既に大変よく効いていることが分かりました。

static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
  {
    LhsEval lhsEval(lhs);
    
    Index n = lhs.outerSize();
#ifdef EIGEN_HAS_OPENMP
    Eigen::initParallel();
    Index threads = Eigen::nbThreads();
#endif
    
    for(Index c=0; c1 && lhsEval.nonZerosEstimate() > 20000)
      {
        #pragma omp parallel for schedule(dynamic,(n+threads*4-1)/(threads*4)) num_threads(threads)
        for(Index i=0; i
  
  
  
 疎行列とベクタの積は、SpMVというみたいです。いくつか論文がヒットし、実装を見てみました。 
  
[PDF] Autotuning Runtime Specialization For Sparse Matrix-Vector Multiplication by Buse Yilmaz, Baris Aktemur, M. Garzarán, Samuel N. Kamin, Mustafa Furkan Kıraç · 10.1145/2851500 (citationsy.com)
A sparse matrix‐vector multiplication method with low preprocessing cost - Aktemur - 2018 - Concurrency and Computation: Practice and Experience - Wiley Online Library
動的な並列実行機構を用いたSpMV実装の性能評価 | 文献情報 | J-GLOBAL 科学技術総合リンクセンター (jst.go.jp)
疎行列の問題は、疎であるが故に、メモリからのロードがネックになり最近のCPUの演算能力を生かせないことにあります。上記のような試みはなされていますが、本質的には、解決が難しいということです。Eigenのコードにしても、残る改善の余地は、ループアンローリングしかありません。しかし、これをC++で記述しても、Eigen程の精緻なコードを生成しないので、アセンブラで記述するしかないでしょう。
アンローリング｜行列積高速化#13｜やまもと｜note
アンローリングは、行うことにしても必要な性能(IPMの3倍速)は、得られそうにありません。そこで、さらなる改善としてNSPに特化した行列演算を試みることにしました。SIMDの命令セットを見たのですが、残念ながら効力があるだろうと思われるのは、AVXという命令セット以降です。最近のCPUならば、備えているとは思いますが、以降の検討は、AVXにすることにしました。（巨大インスタンスについては、それ以外に立ち向かう方法を思いつきません。）
記述は、アセンブラで行いますが、VisualStudioのアセンブラの使い方が良く分からず、JITアセンブラXbyakを使います。試してみたのですが、良い感じで使えそうでした。
また、現在の開発環境は古いCPUでAVXの評価が出来ないので、中古サーバを手配し開発環境を更新します。
にても、大規模NSP問題を解決するために、FirstOrderSolverを導入、FirstOrderSolverの問題を解決するために、疎行列演算問題を解く羽目に。疎行列問題を解決するために、SIMD
に手を出すことになるとは思ってもみませんでした。