コレスキー分解のSIMD化

 ボトルネックになっている部分をVisual　Studio　パフォーマンスプロファイラで調べ特定します。その部分のアセンブラコードを調べ、コンパイラが賢くAVX2命令を生成しているかどうかをチェックしました。意図したようにコンパイラが生成してくれればよいのですが、そうでない場合は、IntrinsicsCodeを書くことになります。（コンパイラは、ymmレジスタを使っていませんでした。）以下は、コンパイラの生成が思わしくないので自分で書いた唯一のコードです。

               for (i = 0; i < BLOCK; i += 4) {

               	   __m256d ymm_t0 = _mm256_load_pd(&aa[i + 0 + 0 * BLOCK]);
                   __m256d ymm_t1 = _mm256_load_pd(&aa[i + 0 + 1 * BLOCK]);
                   __m256d ymm_t2 = _mm256_load_pd(&aa[i + 0 + 2 * BLOCK]);
                   __m256d ymm_t3 = _mm256_load_pd(&aa[i + 0 + 3 * BLOCK]);
 
           

                    __m256d ymm_aUnderNow = _mm256_load_pd(aUnderNow);
                    
                    __m256d ymm_aboveNow0 = _mm256_broadcast_sd(  aboveNow);
                    __m256d ymm_aboveNow1 = _mm256_broadcast_sd(aboveNow + 1);
                    __m256d ymm_aboveNow2 = _mm256_broadcast_sd(aboveNow + 2);
                    __m256d ymm_aboveNow3 = _mm256_broadcast_sd(aboveNow + 3);

                    
                    __m256d ymm_tmp = _mm256_broadcast_sd(&(work[k]));
                    __m256d a=_mm256_mul_pd(ymm_aUnderNow, ymm_tmp);

                    ymm_t0 = _mm256_fnmadd_pd(a, ymm_aboveNow0,ymm_t0);//	c=－a＊b＋c; c-=ab;
                    ymm_t1 = _mm256_fnmadd_pd(a, ymm_aboveNow1, ymm_t1);//	c=－a＊b＋c; c-=ab;
                    ymm_t2 = _mm256_fnmadd_pd(a, ymm_aboveNow2, ymm_t2);//	c=－a＊b＋c; c-=ab;
                    ymm_t3 = _mm256_fnmadd_pd(a, ymm_aboveNow3, ymm_t3);//	c=－a＊b＋c; c-=ab;

                    aUnderNow += BLOCK;
                    aboveNow += BLOCK;
                }

                _mm256_store_pd(&aa[i + 0 + 0 * BLOCK], ymm_t0);

BLOCKは定数であり、コンパイラは、勝手にループ展開していましたが全く必要のないことです。（仮想的マシンの話ですが、同時命令実行数が4のとき、上記コードは、ベクタライズなしの同時命令実行数1に対して4x4=16倍の速度になります。） それはともかく、こんな風に書き換えるだけで10%程度速度が向上しました。