Help us understand the problem. What is going on with this article?

100万回の小規模行列の掛け算におけるメモリ配置と計算速度について

More than 1 year has passed since last update.

C++を使っておりますが、テーマが大分それてしまいました。
カレンダーで、このタイトルだけ趣旨がずれて浮いてる気がしてます…失礼いたします…

本稿でやったこと

 大量の小規模行列の掛け算に対し、2種類のメモリ配置を試し、どちらが計算速度が速いかを確認しました。
 いつも通り設計していたのですが、速度向上を図って試行錯誤していくうちに、C言語のようなソースコードになってしまいました。いずれ整えます。
 私はHPCの分野は素人ですので、ご意見頂けると有り難いです。

背景

 ※これはテーマからかなり逸れるので、飛ばしても構いません。

 多くの工学分野の解析実務で、FEM(有限要素法)という数値解析手法が広く用いられています。FEMは振動解析や郷土解析などの構造計算の分野でも使われます。私は構造計算でFEMを使うことが多いため、今回それを前提とします。
 FEMによる数値解析においては、限られたリソースで高速に解く必要があり、効率の良いメモリ配置が重要と思います。FEMの書籍には、解説とともにサンプルプログラムが公開されているものがあります(例えば[1])。しかしその解説において、計算速度に対するメモリ配置の説明が無いと私は思っております。そこで、現状、私が想定している2種類のメモリ配置の方法について計算速度を確認し、今後の機会に活かそうかと思った次第です。
 ただし、HPCの分野でよく見られる大規模疎行列の計算の高速化ではなく、大量の小規模行列の演算に着目します。FEMの用途によっては、大量の小規模行列の演算が必要になるからです。

FEMによる構造計算の手順の概説

 ※これもテーマからかなり逸れるので、飛ばしても構いません。

 例えば周波数応答の振動解析の場合、一般的な計算手順はおおよそ下記ようになります。

 1. 解析対象の物体を$n$個の要素に分割する。
 2. 着目したい周波数$ \omega_{j} $について次の計算を行う。
  2-1. 要素の係数行列$ \left[K_e^{\,i}\,\right](i=1,...,n)$を計算する。
  2-2. n個の$\left[K_e^{\,i}\,\right]$を全体係数行列$ \left[K_{global}\,\right]$にAssembleする。
  2-3. $ \left[K_{global}\,\right]$を左辺行列、入力を右辺ベクトルとした連立一次方程式を解き、出力を求める。
  2-4. 2-3の出力と、2-1の計算過程で得られた行列を用いて要素の特性値(応力、歪みなど)を求める。

 2-3の処理は、大規模データの場合、多くの解析プログラムにおいてクリティカルパスになることが多く、高速化の研究対象になっています。ここでは、余り議論の対象とならない2-1、2-4などの小規模行列の処理に着目します。現状この処理がプログラム全体の計算時間に占める割合は大きくないと思いますが、2-3が高速化されれば必然的にこの処理の割合が増えるので、検討の意義はあるかと思います。2-1の$ \left[K_e^{\,i}\,\right]$は下式で計算されます。

$$
\left[K_e^{\,i}\,\right]=\int_{\Omega_e^{\,i}}\left[B_e^{\,i}\,\right]^T\left[D_e^{\,i}\,\right]\left[B_e^{\,i}\,\right]d\Omega \tag{1}
$$

3つの行列の掛け算を積分した形になります。各行列が何を表すかの説明は省略します。
今回、100万要素の計算を前提とし、(1)式をi=1~100万まで100万回計算します。以前なら大規模、今なら中規模という要素数です。

メモリ配置2案

 (1)式は3つの行列の掛け算と積分ですが、行列掛け算は単純化のため3つでなく2つとし、積分は数値積分で行うものとし、検討対象は次式とします。
$$\left[A^{\,i}\,\right]\left[B^{\,i}\,\right]=\left[X^{\,i}\,\right] \tag{2} $$
行列サイズはそれぞれ6x6、6x24、6x24とします。これはFEMのHEXA要素の計算の一部を想定しています。
$\left[A^{\,i}\right]$、$\left[B^{\,i}\right]$、$\left[X^{\,i}\right]$を下記2通りの方法(a),(b)で配置します。

 (a) ABXを1セットとして、n=100万回連続して格納する。
  $[A^{\,1}]\,[B^{\,1}]\,[X^{\,1}]\,\,,[A^{\,2}]\,[B^{\,2}]\,[X^{\,2}]\,\,,...,\,\,[A^{\,n}]\,[B^{\,n}]\,[X^{\,n}] \tag{3} $

 (b) A,B,Xをそれぞれn=100万回連続して格納する。
     $[A^{\,1}]\,[A^{\,2}]\,\,,...,\,\,[A^{\,n}]$、 $[B^{\,1}]\,[B^{\,2}]\,\,,...,\,\,[B^{\,n}]$、 $[X^{\,1}]\,[X^{\,2}]\,\,,...,\,\,[X^{\,n}] \tag{4} $

 $[A^{\,i}]$、$[B^{\,i}]$、$[X^{\,i}]$合わせて倍精度浮動小数点が324個なので、100万個分のメモリ使用量は約2.4GBとなります。キャッシュには乗らないので、メモリ配置の良し悪しでキャッシュヒット率が決まり、速度が決まる、という目論見です。(a)は一回の演算に使う行列ABXが連続して格納されているので、キャッシュヒット率が高くなり高速であると期待しておりました。

実装

 まず、100万要素分のdouble型配列を準備します。

    constexpr size_t    nElm    = 1000000;  //  100万要素

    constexpr int       nAelm   = 36;       //  行列A(6x6)の成分数
    constexpr int       nBelm   = 144;      //  行列B(6x24)の成分数
    constexpr int       nXelm   = 144;      //  行列X(6x24)の成分数
    constexpr int       nElmAll = nBelm + nAelm + nXelm;    //  ABXの成分数の和324個

    constexpr size_t    nSize = nElm * nElmAll;         //  100万要素の行列成分数
    std::vector<double> buf( nSize );                   //  100万要素分のdouble型配列。メモリ使用量は nSize * 8byte = 約2.4GB

値は乱数で初期化しておきます。

    //  -1.0~1.0の乱数設定
    {
        std::random_device  seed_gen;
        std::mt19937        engine( seed_gen() );

        for( size_t i = 0; i < nSize; ++i ){
            buf[ i ] = ( 2.0 * double( engine() ) / double( UINT32_MAX ) ) - 1.0;
        }
    }

固定長の2つの行列掛け算なので、ループは全て展開します。(a)(b)それぞれのAB=Xの掛け算関数の実装を次に示します。(a)では、A行列の先頭ポインタが分かればB、Xの先頭ポインタも分かるので、引数はAの先頭ポインタのみです。一方、(b)については、A,B,Xそれぞれの先頭ポインタが必要です。
(a)(b)の計算の実装はそれぞれ次のようになりました。ソースが長いので閉じております。

(a)の行列掛け算の実装
inline bool mmultABX_a( double* pA )
{
    {
        pA[ 180 ] = pA[ 0 ] * pA[ 36 ] + pA[ 6 ] * pA[ 37 ] + pA[ 12 ] * pA[ 38 ] + pA[ 18 ] * pA[ 39 ] + pA[ 24 ] * pA[ 40 ] + pA[ 30 ] * pA[ 41 ];
        pA[ 181 ] = pA[ 1 ] * pA[ 36 ] + pA[ 7 ] * pA[ 37 ] + pA[ 13 ] * pA[ 38 ] + pA[ 19 ] * pA[ 39 ] + pA[ 25 ] * pA[ 40 ] + pA[ 31 ] * pA[ 41 ];
        pA[ 182 ] = pA[ 2 ] * pA[ 36 ] + pA[ 8 ] * pA[ 37 ] + pA[ 14 ] * pA[ 38 ] + pA[ 20 ] * pA[ 39 ] + pA[ 26 ] * pA[ 40 ] + pA[ 32 ] * pA[ 41 ];
        pA[ 183 ] = pA[ 3 ] * pA[ 36 ] + pA[ 9 ] * pA[ 37 ] + pA[ 15 ] * pA[ 38 ] + pA[ 21 ] * pA[ 39 ] + pA[ 27 ] * pA[ 40 ] + pA[ 33 ] * pA[ 41 ];
        pA[ 184 ] = pA[ 4 ] * pA[ 36 ] + pA[ 10 ] * pA[ 37 ] + pA[ 16 ] * pA[ 38 ] + pA[ 22 ] * pA[ 39 ] + pA[ 28 ] * pA[ 40 ] + pA[ 34 ] * pA[ 41 ];
        pA[ 185 ] = pA[ 5 ] * pA[ 36 ] + pA[ 11 ] * pA[ 37 ] + pA[ 17 ] * pA[ 38 ] + pA[ 23 ] * pA[ 39 ] + pA[ 29 ] * pA[ 40 ] + pA[ 35 ] * pA[ 41 ];
    }
    {
        pA[ 186 ] = pA[ 0 ] * pA[ 42 ] + pA[ 6 ] * pA[ 43 ] + pA[ 12 ] * pA[ 44 ] + pA[ 18 ] * pA[ 45 ] + pA[ 24 ] * pA[ 46 ] + pA[ 30 ] * pA[ 47 ];
        pA[ 187 ] = pA[ 1 ] * pA[ 42 ] + pA[ 7 ] * pA[ 43 ] + pA[ 13 ] * pA[ 44 ] + pA[ 19 ] * pA[ 45 ] + pA[ 25 ] * pA[ 46 ] + pA[ 31 ] * pA[ 47 ];
        pA[ 188 ] = pA[ 2 ] * pA[ 42 ] + pA[ 8 ] * pA[ 43 ] + pA[ 14 ] * pA[ 44 ] + pA[ 20 ] * pA[ 45 ] + pA[ 26 ] * pA[ 46 ] + pA[ 32 ] * pA[ 47 ];
        pA[ 189 ] = pA[ 3 ] * pA[ 42 ] + pA[ 9 ] * pA[ 43 ] + pA[ 15 ] * pA[ 44 ] + pA[ 21 ] * pA[ 45 ] + pA[ 27 ] * pA[ 46 ] + pA[ 33 ] * pA[ 47 ];
        pA[ 190 ] = pA[ 4 ] * pA[ 42 ] + pA[ 10 ] * pA[ 43 ] + pA[ 16 ] * pA[ 44 ] + pA[ 22 ] * pA[ 45 ] + pA[ 28 ] * pA[ 46 ] + pA[ 34 ] * pA[ 47 ];
        pA[ 191 ] = pA[ 5 ] * pA[ 42 ] + pA[ 11 ] * pA[ 43 ] + pA[ 17 ] * pA[ 44 ] + pA[ 23 ] * pA[ 45 ] + pA[ 29 ] * pA[ 46 ] + pA[ 35 ] * pA[ 47 ];
    }
    {
        pA[ 192 ] = pA[ 0 ] * pA[ 48 ] + pA[ 6 ] * pA[ 49 ] + pA[ 12 ] * pA[ 50 ] + pA[ 18 ] * pA[ 51 ] + pA[ 24 ] * pA[ 52 ] + pA[ 30 ] * pA[ 53 ];
        pA[ 193 ] = pA[ 1 ] * pA[ 48 ] + pA[ 7 ] * pA[ 49 ] + pA[ 13 ] * pA[ 50 ] + pA[ 19 ] * pA[ 51 ] + pA[ 25 ] * pA[ 52 ] + pA[ 31 ] * pA[ 53 ];
        pA[ 194 ] = pA[ 2 ] * pA[ 48 ] + pA[ 8 ] * pA[ 49 ] + pA[ 14 ] * pA[ 50 ] + pA[ 20 ] * pA[ 51 ] + pA[ 26 ] * pA[ 52 ] + pA[ 32 ] * pA[ 53 ];
        pA[ 195 ] = pA[ 3 ] * pA[ 48 ] + pA[ 9 ] * pA[ 49 ] + pA[ 15 ] * pA[ 50 ] + pA[ 21 ] * pA[ 51 ] + pA[ 27 ] * pA[ 52 ] + pA[ 33 ] * pA[ 53 ];
        pA[ 196 ] = pA[ 4 ] * pA[ 48 ] + pA[ 10 ] * pA[ 49 ] + pA[ 16 ] * pA[ 50 ] + pA[ 22 ] * pA[ 51 ] + pA[ 28 ] * pA[ 52 ] + pA[ 34 ] * pA[ 53 ];
        pA[ 197 ] = pA[ 5 ] * pA[ 48 ] + pA[ 11 ] * pA[ 49 ] + pA[ 17 ] * pA[ 50 ] + pA[ 23 ] * pA[ 51 ] + pA[ 29 ] * pA[ 52 ] + pA[ 35 ] * pA[ 53 ];
    }
    {
        pA[ 198 ] = pA[ 0 ] * pA[ 54 ] + pA[ 6 ] * pA[ 55 ] + pA[ 12 ] * pA[ 56 ] + pA[ 18 ] * pA[ 57 ] + pA[ 24 ] * pA[ 58 ] + pA[ 30 ] * pA[ 59 ];
        pA[ 199 ] = pA[ 1 ] * pA[ 54 ] + pA[ 7 ] * pA[ 55 ] + pA[ 13 ] * pA[ 56 ] + pA[ 19 ] * pA[ 57 ] + pA[ 25 ] * pA[ 58 ] + pA[ 31 ] * pA[ 59 ];
        pA[ 200 ] = pA[ 2 ] * pA[ 54 ] + pA[ 8 ] * pA[ 55 ] + pA[ 14 ] * pA[ 56 ] + pA[ 20 ] * pA[ 57 ] + pA[ 26 ] * pA[ 58 ] + pA[ 32 ] * pA[ 59 ];
        pA[ 201 ] = pA[ 3 ] * pA[ 54 ] + pA[ 9 ] * pA[ 55 ] + pA[ 15 ] * pA[ 56 ] + pA[ 21 ] * pA[ 57 ] + pA[ 27 ] * pA[ 58 ] + pA[ 33 ] * pA[ 59 ];
        pA[ 202 ] = pA[ 4 ] * pA[ 54 ] + pA[ 10 ] * pA[ 55 ] + pA[ 16 ] * pA[ 56 ] + pA[ 22 ] * pA[ 57 ] + pA[ 28 ] * pA[ 58 ] + pA[ 34 ] * pA[ 59 ];
        pA[ 203 ] = pA[ 5 ] * pA[ 54 ] + pA[ 11 ] * pA[ 55 ] + pA[ 17 ] * pA[ 56 ] + pA[ 23 ] * pA[ 57 ] + pA[ 29 ] * pA[ 58 ] + pA[ 35 ] * pA[ 59 ];
    }
    {
        pA[ 204 ] = pA[ 0 ] * pA[ 60 ] + pA[ 6 ] * pA[ 61 ] + pA[ 12 ] * pA[ 62 ] + pA[ 18 ] * pA[ 63 ] + pA[ 24 ] * pA[ 64 ] + pA[ 30 ] * pA[ 65 ];
        pA[ 205 ] = pA[ 1 ] * pA[ 60 ] + pA[ 7 ] * pA[ 61 ] + pA[ 13 ] * pA[ 62 ] + pA[ 19 ] * pA[ 63 ] + pA[ 25 ] * pA[ 64 ] + pA[ 31 ] * pA[ 65 ];
        pA[ 206 ] = pA[ 2 ] * pA[ 60 ] + pA[ 8 ] * pA[ 61 ] + pA[ 14 ] * pA[ 62 ] + pA[ 20 ] * pA[ 63 ] + pA[ 26 ] * pA[ 64 ] + pA[ 32 ] * pA[ 65 ];
        pA[ 207 ] = pA[ 3 ] * pA[ 60 ] + pA[ 9 ] * pA[ 61 ] + pA[ 15 ] * pA[ 62 ] + pA[ 21 ] * pA[ 63 ] + pA[ 27 ] * pA[ 64 ] + pA[ 33 ] * pA[ 65 ];
        pA[ 208 ] = pA[ 4 ] * pA[ 60 ] + pA[ 10 ] * pA[ 61 ] + pA[ 16 ] * pA[ 62 ] + pA[ 22 ] * pA[ 63 ] + pA[ 28 ] * pA[ 64 ] + pA[ 34 ] * pA[ 65 ];
        pA[ 209 ] = pA[ 5 ] * pA[ 60 ] + pA[ 11 ] * pA[ 61 ] + pA[ 17 ] * pA[ 62 ] + pA[ 23 ] * pA[ 63 ] + pA[ 29 ] * pA[ 64 ] + pA[ 35 ] * pA[ 65 ];
    }
    {
        pA[ 210 ] = pA[ 0 ] * pA[ 66 ] + pA[ 6 ] * pA[ 67 ] + pA[ 12 ] * pA[ 68 ] + pA[ 18 ] * pA[ 69 ] + pA[ 24 ] * pA[ 70 ] + pA[ 30 ] * pA[ 71 ];
        pA[ 211 ] = pA[ 1 ] * pA[ 66 ] + pA[ 7 ] * pA[ 67 ] + pA[ 13 ] * pA[ 68 ] + pA[ 19 ] * pA[ 69 ] + pA[ 25 ] * pA[ 70 ] + pA[ 31 ] * pA[ 71 ];
        pA[ 212 ] = pA[ 2 ] * pA[ 66 ] + pA[ 8 ] * pA[ 67 ] + pA[ 14 ] * pA[ 68 ] + pA[ 20 ] * pA[ 69 ] + pA[ 26 ] * pA[ 70 ] + pA[ 32 ] * pA[ 71 ];
        pA[ 213 ] = pA[ 3 ] * pA[ 66 ] + pA[ 9 ] * pA[ 67 ] + pA[ 15 ] * pA[ 68 ] + pA[ 21 ] * pA[ 69 ] + pA[ 27 ] * pA[ 70 ] + pA[ 33 ] * pA[ 71 ];
        pA[ 214 ] = pA[ 4 ] * pA[ 66 ] + pA[ 10 ] * pA[ 67 ] + pA[ 16 ] * pA[ 68 ] + pA[ 22 ] * pA[ 69 ] + pA[ 28 ] * pA[ 70 ] + pA[ 34 ] * pA[ 71 ];
        pA[ 215 ] = pA[ 5 ] * pA[ 66 ] + pA[ 11 ] * pA[ 67 ] + pA[ 17 ] * pA[ 68 ] + pA[ 23 ] * pA[ 69 ] + pA[ 29 ] * pA[ 70 ] + pA[ 35 ] * pA[ 71 ];
    }
    {
        pA[ 216 ] = pA[ 0 ] * pA[ 72 ] + pA[ 6 ] * pA[ 73 ] + pA[ 12 ] * pA[ 74 ] + pA[ 18 ] * pA[ 75 ] + pA[ 24 ] * pA[ 76 ] + pA[ 30 ] * pA[ 77 ];
        pA[ 217 ] = pA[ 1 ] * pA[ 72 ] + pA[ 7 ] * pA[ 73 ] + pA[ 13 ] * pA[ 74 ] + pA[ 19 ] * pA[ 75 ] + pA[ 25 ] * pA[ 76 ] + pA[ 31 ] * pA[ 77 ];
        pA[ 218 ] = pA[ 2 ] * pA[ 72 ] + pA[ 8 ] * pA[ 73 ] + pA[ 14 ] * pA[ 74 ] + pA[ 20 ] * pA[ 75 ] + pA[ 26 ] * pA[ 76 ] + pA[ 32 ] * pA[ 77 ];
        pA[ 219 ] = pA[ 3 ] * pA[ 72 ] + pA[ 9 ] * pA[ 73 ] + pA[ 15 ] * pA[ 74 ] + pA[ 21 ] * pA[ 75 ] + pA[ 27 ] * pA[ 76 ] + pA[ 33 ] * pA[ 77 ];
        pA[ 220 ] = pA[ 4 ] * pA[ 72 ] + pA[ 10 ] * pA[ 73 ] + pA[ 16 ] * pA[ 74 ] + pA[ 22 ] * pA[ 75 ] + pA[ 28 ] * pA[ 76 ] + pA[ 34 ] * pA[ 77 ];
        pA[ 221 ] = pA[ 5 ] * pA[ 72 ] + pA[ 11 ] * pA[ 73 ] + pA[ 17 ] * pA[ 74 ] + pA[ 23 ] * pA[ 75 ] + pA[ 29 ] * pA[ 76 ] + pA[ 35 ] * pA[ 77 ];
    }
    {
        pA[ 222 ] = pA[ 0 ] * pA[ 78 ] + pA[ 6 ] * pA[ 79 ] + pA[ 12 ] * pA[ 80 ] + pA[ 18 ] * pA[ 81 ] + pA[ 24 ] * pA[ 82 ] + pA[ 30 ] * pA[ 83 ];
        pA[ 223 ] = pA[ 1 ] * pA[ 78 ] + pA[ 7 ] * pA[ 79 ] + pA[ 13 ] * pA[ 80 ] + pA[ 19 ] * pA[ 81 ] + pA[ 25 ] * pA[ 82 ] + pA[ 31 ] * pA[ 83 ];
        pA[ 224 ] = pA[ 2 ] * pA[ 78 ] + pA[ 8 ] * pA[ 79 ] + pA[ 14 ] * pA[ 80 ] + pA[ 20 ] * pA[ 81 ] + pA[ 26 ] * pA[ 82 ] + pA[ 32 ] * pA[ 83 ];
        pA[ 225 ] = pA[ 3 ] * pA[ 78 ] + pA[ 9 ] * pA[ 79 ] + pA[ 15 ] * pA[ 80 ] + pA[ 21 ] * pA[ 81 ] + pA[ 27 ] * pA[ 82 ] + pA[ 33 ] * pA[ 83 ];
        pA[ 226 ] = pA[ 4 ] * pA[ 78 ] + pA[ 10 ] * pA[ 79 ] + pA[ 16 ] * pA[ 80 ] + pA[ 22 ] * pA[ 81 ] + pA[ 28 ] * pA[ 82 ] + pA[ 34 ] * pA[ 83 ];
        pA[ 227 ] = pA[ 5 ] * pA[ 78 ] + pA[ 11 ] * pA[ 79 ] + pA[ 17 ] * pA[ 80 ] + pA[ 23 ] * pA[ 81 ] + pA[ 29 ] * pA[ 82 ] + pA[ 35 ] * pA[ 83 ];
    }
    {
        pA[ 228 ] = pA[ 0 ] * pA[ 84 ] + pA[ 6 ] * pA[ 85 ] + pA[ 12 ] * pA[ 86 ] + pA[ 18 ] * pA[ 87 ] + pA[ 24 ] * pA[ 88 ] + pA[ 30 ] * pA[ 89 ];
        pA[ 229 ] = pA[ 1 ] * pA[ 84 ] + pA[ 7 ] * pA[ 85 ] + pA[ 13 ] * pA[ 86 ] + pA[ 19 ] * pA[ 87 ] + pA[ 25 ] * pA[ 88 ] + pA[ 31 ] * pA[ 89 ];
        pA[ 230 ] = pA[ 2 ] * pA[ 84 ] + pA[ 8 ] * pA[ 85 ] + pA[ 14 ] * pA[ 86 ] + pA[ 20 ] * pA[ 87 ] + pA[ 26 ] * pA[ 88 ] + pA[ 32 ] * pA[ 89 ];
        pA[ 231 ] = pA[ 3 ] * pA[ 84 ] + pA[ 9 ] * pA[ 85 ] + pA[ 15 ] * pA[ 86 ] + pA[ 21 ] * pA[ 87 ] + pA[ 27 ] * pA[ 88 ] + pA[ 33 ] * pA[ 89 ];
        pA[ 232 ] = pA[ 4 ] * pA[ 84 ] + pA[ 10 ] * pA[ 85 ] + pA[ 16 ] * pA[ 86 ] + pA[ 22 ] * pA[ 87 ] + pA[ 28 ] * pA[ 88 ] + pA[ 34 ] * pA[ 89 ];
        pA[ 233 ] = pA[ 5 ] * pA[ 84 ] + pA[ 11 ] * pA[ 85 ] + pA[ 17 ] * pA[ 86 ] + pA[ 23 ] * pA[ 87 ] + pA[ 29 ] * pA[ 88 ] + pA[ 35 ] * pA[ 89 ];
    }
    {
        pA[ 234 ] = pA[ 0 ] * pA[ 90 ] + pA[ 6 ] * pA[ 91 ] + pA[ 12 ] * pA[ 92 ] + pA[ 18 ] * pA[ 93 ] + pA[ 24 ] * pA[ 94 ] + pA[ 30 ] * pA[ 95 ];
        pA[ 235 ] = pA[ 1 ] * pA[ 90 ] + pA[ 7 ] * pA[ 91 ] + pA[ 13 ] * pA[ 92 ] + pA[ 19 ] * pA[ 93 ] + pA[ 25 ] * pA[ 94 ] + pA[ 31 ] * pA[ 95 ];
        pA[ 236 ] = pA[ 2 ] * pA[ 90 ] + pA[ 8 ] * pA[ 91 ] + pA[ 14 ] * pA[ 92 ] + pA[ 20 ] * pA[ 93 ] + pA[ 26 ] * pA[ 94 ] + pA[ 32 ] * pA[ 95 ];
        pA[ 237 ] = pA[ 3 ] * pA[ 90 ] + pA[ 9 ] * pA[ 91 ] + pA[ 15 ] * pA[ 92 ] + pA[ 21 ] * pA[ 93 ] + pA[ 27 ] * pA[ 94 ] + pA[ 33 ] * pA[ 95 ];
        pA[ 238 ] = pA[ 4 ] * pA[ 90 ] + pA[ 10 ] * pA[ 91 ] + pA[ 16 ] * pA[ 92 ] + pA[ 22 ] * pA[ 93 ] + pA[ 28 ] * pA[ 94 ] + pA[ 34 ] * pA[ 95 ];
        pA[ 239 ] = pA[ 5 ] * pA[ 90 ] + pA[ 11 ] * pA[ 91 ] + pA[ 17 ] * pA[ 92 ] + pA[ 23 ] * pA[ 93 ] + pA[ 29 ] * pA[ 94 ] + pA[ 35 ] * pA[ 95 ];
    }
    {
        pA[ 240 ] = pA[ 0 ] * pA[ 96 ] + pA[ 6 ] * pA[ 97 ] + pA[ 12 ] * pA[ 98 ] + pA[ 18 ] * pA[ 99 ] + pA[ 24 ] * pA[ 100 ] + pA[ 30 ] * pA[ 101 ];
        pA[ 241 ] = pA[ 1 ] * pA[ 96 ] + pA[ 7 ] * pA[ 97 ] + pA[ 13 ] * pA[ 98 ] + pA[ 19 ] * pA[ 99 ] + pA[ 25 ] * pA[ 100 ] + pA[ 31 ] * pA[ 101 ];
        pA[ 242 ] = pA[ 2 ] * pA[ 96 ] + pA[ 8 ] * pA[ 97 ] + pA[ 14 ] * pA[ 98 ] + pA[ 20 ] * pA[ 99 ] + pA[ 26 ] * pA[ 100 ] + pA[ 32 ] * pA[ 101 ];
        pA[ 243 ] = pA[ 3 ] * pA[ 96 ] + pA[ 9 ] * pA[ 97 ] + pA[ 15 ] * pA[ 98 ] + pA[ 21 ] * pA[ 99 ] + pA[ 27 ] * pA[ 100 ] + pA[ 33 ] * pA[ 101 ];
        pA[ 244 ] = pA[ 4 ] * pA[ 96 ] + pA[ 10 ] * pA[ 97 ] + pA[ 16 ] * pA[ 98 ] + pA[ 22 ] * pA[ 99 ] + pA[ 28 ] * pA[ 100 ] + pA[ 34 ] * pA[ 101 ];
        pA[ 245 ] = pA[ 5 ] * pA[ 96 ] + pA[ 11 ] * pA[ 97 ] + pA[ 17 ] * pA[ 98 ] + pA[ 23 ] * pA[ 99 ] + pA[ 29 ] * pA[ 100 ] + pA[ 35 ] * pA[ 101 ];
    }
    {
        pA[ 246 ] = pA[ 0 ] * pA[ 102 ] + pA[ 6 ] * pA[ 103 ] + pA[ 12 ] * pA[ 104 ] + pA[ 18 ] * pA[ 105 ] + pA[ 24 ] * pA[ 106 ] + pA[ 30 ] * pA[ 107 ];
        pA[ 247 ] = pA[ 1 ] * pA[ 102 ] + pA[ 7 ] * pA[ 103 ] + pA[ 13 ] * pA[ 104 ] + pA[ 19 ] * pA[ 105 ] + pA[ 25 ] * pA[ 106 ] + pA[ 31 ] * pA[ 107 ];
        pA[ 248 ] = pA[ 2 ] * pA[ 102 ] + pA[ 8 ] * pA[ 103 ] + pA[ 14 ] * pA[ 104 ] + pA[ 20 ] * pA[ 105 ] + pA[ 26 ] * pA[ 106 ] + pA[ 32 ] * pA[ 107 ];
        pA[ 249 ] = pA[ 3 ] * pA[ 102 ] + pA[ 9 ] * pA[ 103 ] + pA[ 15 ] * pA[ 104 ] + pA[ 21 ] * pA[ 105 ] + pA[ 27 ] * pA[ 106 ] + pA[ 33 ] * pA[ 107 ];
        pA[ 250 ] = pA[ 4 ] * pA[ 102 ] + pA[ 10 ] * pA[ 103 ] + pA[ 16 ] * pA[ 104 ] + pA[ 22 ] * pA[ 105 ] + pA[ 28 ] * pA[ 106 ] + pA[ 34 ] * pA[ 107 ];
        pA[ 251 ] = pA[ 5 ] * pA[ 102 ] + pA[ 11 ] * pA[ 103 ] + pA[ 17 ] * pA[ 104 ] + pA[ 23 ] * pA[ 105 ] + pA[ 29 ] * pA[ 106 ] + pA[ 35 ] * pA[ 107 ];
    }
    {
        pA[ 252 ] = pA[ 0 ] * pA[ 108 ] + pA[ 6 ] * pA[ 109 ] + pA[ 12 ] * pA[ 110 ] + pA[ 18 ] * pA[ 111 ] + pA[ 24 ] * pA[ 112 ] + pA[ 30 ] * pA[ 113 ];
        pA[ 253 ] = pA[ 1 ] * pA[ 108 ] + pA[ 7 ] * pA[ 109 ] + pA[ 13 ] * pA[ 110 ] + pA[ 19 ] * pA[ 111 ] + pA[ 25 ] * pA[ 112 ] + pA[ 31 ] * pA[ 113 ];
        pA[ 254 ] = pA[ 2 ] * pA[ 108 ] + pA[ 8 ] * pA[ 109 ] + pA[ 14 ] * pA[ 110 ] + pA[ 20 ] * pA[ 111 ] + pA[ 26 ] * pA[ 112 ] + pA[ 32 ] * pA[ 113 ];
        pA[ 255 ] = pA[ 3 ] * pA[ 108 ] + pA[ 9 ] * pA[ 109 ] + pA[ 15 ] * pA[ 110 ] + pA[ 21 ] * pA[ 111 ] + pA[ 27 ] * pA[ 112 ] + pA[ 33 ] * pA[ 113 ];
        pA[ 256 ] = pA[ 4 ] * pA[ 108 ] + pA[ 10 ] * pA[ 109 ] + pA[ 16 ] * pA[ 110 ] + pA[ 22 ] * pA[ 111 ] + pA[ 28 ] * pA[ 112 ] + pA[ 34 ] * pA[ 113 ];
        pA[ 257 ] = pA[ 5 ] * pA[ 108 ] + pA[ 11 ] * pA[ 109 ] + pA[ 17 ] * pA[ 110 ] + pA[ 23 ] * pA[ 111 ] + pA[ 29 ] * pA[ 112 ] + pA[ 35 ] * pA[ 113 ];
    }
    {
        pA[ 258 ] = pA[ 0 ] * pA[ 114 ] + pA[ 6 ] * pA[ 115 ] + pA[ 12 ] * pA[ 116 ] + pA[ 18 ] * pA[ 117 ] + pA[ 24 ] * pA[ 118 ] + pA[ 30 ] * pA[ 119 ];
        pA[ 259 ] = pA[ 1 ] * pA[ 114 ] + pA[ 7 ] * pA[ 115 ] + pA[ 13 ] * pA[ 116 ] + pA[ 19 ] * pA[ 117 ] + pA[ 25 ] * pA[ 118 ] + pA[ 31 ] * pA[ 119 ];
        pA[ 260 ] = pA[ 2 ] * pA[ 114 ] + pA[ 8 ] * pA[ 115 ] + pA[ 14 ] * pA[ 116 ] + pA[ 20 ] * pA[ 117 ] + pA[ 26 ] * pA[ 118 ] + pA[ 32 ] * pA[ 119 ];
        pA[ 261 ] = pA[ 3 ] * pA[ 114 ] + pA[ 9 ] * pA[ 115 ] + pA[ 15 ] * pA[ 116 ] + pA[ 21 ] * pA[ 117 ] + pA[ 27 ] * pA[ 118 ] + pA[ 33 ] * pA[ 119 ];
        pA[ 262 ] = pA[ 4 ] * pA[ 114 ] + pA[ 10 ] * pA[ 115 ] + pA[ 16 ] * pA[ 116 ] + pA[ 22 ] * pA[ 117 ] + pA[ 28 ] * pA[ 118 ] + pA[ 34 ] * pA[ 119 ];
        pA[ 263 ] = pA[ 5 ] * pA[ 114 ] + pA[ 11 ] * pA[ 115 ] + pA[ 17 ] * pA[ 116 ] + pA[ 23 ] * pA[ 117 ] + pA[ 29 ] * pA[ 118 ] + pA[ 35 ] * pA[ 119 ];
    }
    {
        pA[ 264 ] = pA[ 0 ] * pA[ 120 ] + pA[ 6 ] * pA[ 121 ] + pA[ 12 ] * pA[ 122 ] + pA[ 18 ] * pA[ 123 ] + pA[ 24 ] * pA[ 124 ] + pA[ 30 ] * pA[ 125 ];
        pA[ 265 ] = pA[ 1 ] * pA[ 120 ] + pA[ 7 ] * pA[ 121 ] + pA[ 13 ] * pA[ 122 ] + pA[ 19 ] * pA[ 123 ] + pA[ 25 ] * pA[ 124 ] + pA[ 31 ] * pA[ 125 ];
        pA[ 266 ] = pA[ 2 ] * pA[ 120 ] + pA[ 8 ] * pA[ 121 ] + pA[ 14 ] * pA[ 122 ] + pA[ 20 ] * pA[ 123 ] + pA[ 26 ] * pA[ 124 ] + pA[ 32 ] * pA[ 125 ];
        pA[ 267 ] = pA[ 3 ] * pA[ 120 ] + pA[ 9 ] * pA[ 121 ] + pA[ 15 ] * pA[ 122 ] + pA[ 21 ] * pA[ 123 ] + pA[ 27 ] * pA[ 124 ] + pA[ 33 ] * pA[ 125 ];
        pA[ 268 ] = pA[ 4 ] * pA[ 120 ] + pA[ 10 ] * pA[ 121 ] + pA[ 16 ] * pA[ 122 ] + pA[ 22 ] * pA[ 123 ] + pA[ 28 ] * pA[ 124 ] + pA[ 34 ] * pA[ 125 ];
        pA[ 269 ] = pA[ 5 ] * pA[ 120 ] + pA[ 11 ] * pA[ 121 ] + pA[ 17 ] * pA[ 122 ] + pA[ 23 ] * pA[ 123 ] + pA[ 29 ] * pA[ 124 ] + pA[ 35 ] * pA[ 125 ];
    }
    {
        pA[ 270 ] = pA[ 0 ] * pA[ 126 ] + pA[ 6 ] * pA[ 127 ] + pA[ 12 ] * pA[ 128 ] + pA[ 18 ] * pA[ 129 ] + pA[ 24 ] * pA[ 130 ] + pA[ 30 ] * pA[ 131 ];
        pA[ 271 ] = pA[ 1 ] * pA[ 126 ] + pA[ 7 ] * pA[ 127 ] + pA[ 13 ] * pA[ 128 ] + pA[ 19 ] * pA[ 129 ] + pA[ 25 ] * pA[ 130 ] + pA[ 31 ] * pA[ 131 ];
        pA[ 272 ] = pA[ 2 ] * pA[ 126 ] + pA[ 8 ] * pA[ 127 ] + pA[ 14 ] * pA[ 128 ] + pA[ 20 ] * pA[ 129 ] + pA[ 26 ] * pA[ 130 ] + pA[ 32 ] * pA[ 131 ];
        pA[ 273 ] = pA[ 3 ] * pA[ 126 ] + pA[ 9 ] * pA[ 127 ] + pA[ 15 ] * pA[ 128 ] + pA[ 21 ] * pA[ 129 ] + pA[ 27 ] * pA[ 130 ] + pA[ 33 ] * pA[ 131 ];
        pA[ 274 ] = pA[ 4 ] * pA[ 126 ] + pA[ 10 ] * pA[ 127 ] + pA[ 16 ] * pA[ 128 ] + pA[ 22 ] * pA[ 129 ] + pA[ 28 ] * pA[ 130 ] + pA[ 34 ] * pA[ 131 ];
        pA[ 275 ] = pA[ 5 ] * pA[ 126 ] + pA[ 11 ] * pA[ 127 ] + pA[ 17 ] * pA[ 128 ] + pA[ 23 ] * pA[ 129 ] + pA[ 29 ] * pA[ 130 ] + pA[ 35 ] * pA[ 131 ];
    }
    {
        pA[ 276 ] = pA[ 0 ] * pA[ 132 ] + pA[ 6 ] * pA[ 133 ] + pA[ 12 ] * pA[ 134 ] + pA[ 18 ] * pA[ 135 ] + pA[ 24 ] * pA[ 136 ] + pA[ 30 ] * pA[ 137 ];
        pA[ 277 ] = pA[ 1 ] * pA[ 132 ] + pA[ 7 ] * pA[ 133 ] + pA[ 13 ] * pA[ 134 ] + pA[ 19 ] * pA[ 135 ] + pA[ 25 ] * pA[ 136 ] + pA[ 31 ] * pA[ 137 ];
        pA[ 278 ] = pA[ 2 ] * pA[ 132 ] + pA[ 8 ] * pA[ 133 ] + pA[ 14 ] * pA[ 134 ] + pA[ 20 ] * pA[ 135 ] + pA[ 26 ] * pA[ 136 ] + pA[ 32 ] * pA[ 137 ];
        pA[ 279 ] = pA[ 3 ] * pA[ 132 ] + pA[ 9 ] * pA[ 133 ] + pA[ 15 ] * pA[ 134 ] + pA[ 21 ] * pA[ 135 ] + pA[ 27 ] * pA[ 136 ] + pA[ 33 ] * pA[ 137 ];
        pA[ 280 ] = pA[ 4 ] * pA[ 132 ] + pA[ 10 ] * pA[ 133 ] + pA[ 16 ] * pA[ 134 ] + pA[ 22 ] * pA[ 135 ] + pA[ 28 ] * pA[ 136 ] + pA[ 34 ] * pA[ 137 ];
        pA[ 281 ] = pA[ 5 ] * pA[ 132 ] + pA[ 11 ] * pA[ 133 ] + pA[ 17 ] * pA[ 134 ] + pA[ 23 ] * pA[ 135 ] + pA[ 29 ] * pA[ 136 ] + pA[ 35 ] * pA[ 137 ];
    }
    {
        pA[ 282 ] = pA[ 0 ] * pA[ 138 ] + pA[ 6 ] * pA[ 139 ] + pA[ 12 ] * pA[ 140 ] + pA[ 18 ] * pA[ 141 ] + pA[ 24 ] * pA[ 142 ] + pA[ 30 ] * pA[ 143 ];
        pA[ 283 ] = pA[ 1 ] * pA[ 138 ] + pA[ 7 ] * pA[ 139 ] + pA[ 13 ] * pA[ 140 ] + pA[ 19 ] * pA[ 141 ] + pA[ 25 ] * pA[ 142 ] + pA[ 31 ] * pA[ 143 ];
        pA[ 284 ] = pA[ 2 ] * pA[ 138 ] + pA[ 8 ] * pA[ 139 ] + pA[ 14 ] * pA[ 140 ] + pA[ 20 ] * pA[ 141 ] + pA[ 26 ] * pA[ 142 ] + pA[ 32 ] * pA[ 143 ];
        pA[ 285 ] = pA[ 3 ] * pA[ 138 ] + pA[ 9 ] * pA[ 139 ] + pA[ 15 ] * pA[ 140 ] + pA[ 21 ] * pA[ 141 ] + pA[ 27 ] * pA[ 142 ] + pA[ 33 ] * pA[ 143 ];
        pA[ 286 ] = pA[ 4 ] * pA[ 138 ] + pA[ 10 ] * pA[ 139 ] + pA[ 16 ] * pA[ 140 ] + pA[ 22 ] * pA[ 141 ] + pA[ 28 ] * pA[ 142 ] + pA[ 34 ] * pA[ 143 ];
        pA[ 287 ] = pA[ 5 ] * pA[ 138 ] + pA[ 11 ] * pA[ 139 ] + pA[ 17 ] * pA[ 140 ] + pA[ 23 ] * pA[ 141 ] + pA[ 29 ] * pA[ 142 ] + pA[ 35 ] * pA[ 143 ];
    }
    {
        pA[ 288 ] = pA[ 0 ] * pA[ 144 ] + pA[ 6 ] * pA[ 145 ] + pA[ 12 ] * pA[ 146 ] + pA[ 18 ] * pA[ 147 ] + pA[ 24 ] * pA[ 148 ] + pA[ 30 ] * pA[ 149 ];
        pA[ 289 ] = pA[ 1 ] * pA[ 144 ] + pA[ 7 ] * pA[ 145 ] + pA[ 13 ] * pA[ 146 ] + pA[ 19 ] * pA[ 147 ] + pA[ 25 ] * pA[ 148 ] + pA[ 31 ] * pA[ 149 ];
        pA[ 290 ] = pA[ 2 ] * pA[ 144 ] + pA[ 8 ] * pA[ 145 ] + pA[ 14 ] * pA[ 146 ] + pA[ 20 ] * pA[ 147 ] + pA[ 26 ] * pA[ 148 ] + pA[ 32 ] * pA[ 149 ];
        pA[ 291 ] = pA[ 3 ] * pA[ 144 ] + pA[ 9 ] * pA[ 145 ] + pA[ 15 ] * pA[ 146 ] + pA[ 21 ] * pA[ 147 ] + pA[ 27 ] * pA[ 148 ] + pA[ 33 ] * pA[ 149 ];
        pA[ 292 ] = pA[ 4 ] * pA[ 144 ] + pA[ 10 ] * pA[ 145 ] + pA[ 16 ] * pA[ 146 ] + pA[ 22 ] * pA[ 147 ] + pA[ 28 ] * pA[ 148 ] + pA[ 34 ] * pA[ 149 ];
        pA[ 293 ] = pA[ 5 ] * pA[ 144 ] + pA[ 11 ] * pA[ 145 ] + pA[ 17 ] * pA[ 146 ] + pA[ 23 ] * pA[ 147 ] + pA[ 29 ] * pA[ 148 ] + pA[ 35 ] * pA[ 149 ];
    }
    {
        pA[ 294 ] = pA[ 0 ] * pA[ 150 ] + pA[ 6 ] * pA[ 151 ] + pA[ 12 ] * pA[ 152 ] + pA[ 18 ] * pA[ 153 ] + pA[ 24 ] * pA[ 154 ] + pA[ 30 ] * pA[ 155 ];
        pA[ 295 ] = pA[ 1 ] * pA[ 150 ] + pA[ 7 ] * pA[ 151 ] + pA[ 13 ] * pA[ 152 ] + pA[ 19 ] * pA[ 153 ] + pA[ 25 ] * pA[ 154 ] + pA[ 31 ] * pA[ 155 ];
        pA[ 296 ] = pA[ 2 ] * pA[ 150 ] + pA[ 8 ] * pA[ 151 ] + pA[ 14 ] * pA[ 152 ] + pA[ 20 ] * pA[ 153 ] + pA[ 26 ] * pA[ 154 ] + pA[ 32 ] * pA[ 155 ];
        pA[ 297 ] = pA[ 3 ] * pA[ 150 ] + pA[ 9 ] * pA[ 151 ] + pA[ 15 ] * pA[ 152 ] + pA[ 21 ] * pA[ 153 ] + pA[ 27 ] * pA[ 154 ] + pA[ 33 ] * pA[ 155 ];
        pA[ 298 ] = pA[ 4 ] * pA[ 150 ] + pA[ 10 ] * pA[ 151 ] + pA[ 16 ] * pA[ 152 ] + pA[ 22 ] * pA[ 153 ] + pA[ 28 ] * pA[ 154 ] + pA[ 34 ] * pA[ 155 ];
        pA[ 299 ] = pA[ 5 ] * pA[ 150 ] + pA[ 11 ] * pA[ 151 ] + pA[ 17 ] * pA[ 152 ] + pA[ 23 ] * pA[ 153 ] + pA[ 29 ] * pA[ 154 ] + pA[ 35 ] * pA[ 155 ];
    }
    {
        pA[ 300 ] = pA[ 0 ] * pA[ 156 ] + pA[ 6 ] * pA[ 157 ] + pA[ 12 ] * pA[ 158 ] + pA[ 18 ] * pA[ 159 ] + pA[ 24 ] * pA[ 160 ] + pA[ 30 ] * pA[ 161 ];
        pA[ 301 ] = pA[ 1 ] * pA[ 156 ] + pA[ 7 ] * pA[ 157 ] + pA[ 13 ] * pA[ 158 ] + pA[ 19 ] * pA[ 159 ] + pA[ 25 ] * pA[ 160 ] + pA[ 31 ] * pA[ 161 ];
        pA[ 302 ] = pA[ 2 ] * pA[ 156 ] + pA[ 8 ] * pA[ 157 ] + pA[ 14 ] * pA[ 158 ] + pA[ 20 ] * pA[ 159 ] + pA[ 26 ] * pA[ 160 ] + pA[ 32 ] * pA[ 161 ];
        pA[ 303 ] = pA[ 3 ] * pA[ 156 ] + pA[ 9 ] * pA[ 157 ] + pA[ 15 ] * pA[ 158 ] + pA[ 21 ] * pA[ 159 ] + pA[ 27 ] * pA[ 160 ] + pA[ 33 ] * pA[ 161 ];
        pA[ 304 ] = pA[ 4 ] * pA[ 156 ] + pA[ 10 ] * pA[ 157 ] + pA[ 16 ] * pA[ 158 ] + pA[ 22 ] * pA[ 159 ] + pA[ 28 ] * pA[ 160 ] + pA[ 34 ] * pA[ 161 ];
        pA[ 305 ] = pA[ 5 ] * pA[ 156 ] + pA[ 11 ] * pA[ 157 ] + pA[ 17 ] * pA[ 158 ] + pA[ 23 ] * pA[ 159 ] + pA[ 29 ] * pA[ 160 ] + pA[ 35 ] * pA[ 161 ];
    }
    {
        pA[ 306 ] = pA[ 0 ] * pA[ 162 ] + pA[ 6 ] * pA[ 163 ] + pA[ 12 ] * pA[ 164 ] + pA[ 18 ] * pA[ 165 ] + pA[ 24 ] * pA[ 166 ] + pA[ 30 ] * pA[ 167 ];
        pA[ 307 ] = pA[ 1 ] * pA[ 162 ] + pA[ 7 ] * pA[ 163 ] + pA[ 13 ] * pA[ 164 ] + pA[ 19 ] * pA[ 165 ] + pA[ 25 ] * pA[ 166 ] + pA[ 31 ] * pA[ 167 ];
        pA[ 308 ] = pA[ 2 ] * pA[ 162 ] + pA[ 8 ] * pA[ 163 ] + pA[ 14 ] * pA[ 164 ] + pA[ 20 ] * pA[ 165 ] + pA[ 26 ] * pA[ 166 ] + pA[ 32 ] * pA[ 167 ];
        pA[ 309 ] = pA[ 3 ] * pA[ 162 ] + pA[ 9 ] * pA[ 163 ] + pA[ 15 ] * pA[ 164 ] + pA[ 21 ] * pA[ 165 ] + pA[ 27 ] * pA[ 166 ] + pA[ 33 ] * pA[ 167 ];
        pA[ 310 ] = pA[ 4 ] * pA[ 162 ] + pA[ 10 ] * pA[ 163 ] + pA[ 16 ] * pA[ 164 ] + pA[ 22 ] * pA[ 165 ] + pA[ 28 ] * pA[ 166 ] + pA[ 34 ] * pA[ 167 ];
        pA[ 311 ] = pA[ 5 ] * pA[ 162 ] + pA[ 11 ] * pA[ 163 ] + pA[ 17 ] * pA[ 164 ] + pA[ 23 ] * pA[ 165 ] + pA[ 29 ] * pA[ 166 ] + pA[ 35 ] * pA[ 167 ];
    }
    {
        pA[ 312 ] = pA[ 0 ] * pA[ 168 ] + pA[ 6 ] * pA[ 169 ] + pA[ 12 ] * pA[ 170 ] + pA[ 18 ] * pA[ 171 ] + pA[ 24 ] * pA[ 172 ] + pA[ 30 ] * pA[ 173 ];
        pA[ 313 ] = pA[ 1 ] * pA[ 168 ] + pA[ 7 ] * pA[ 169 ] + pA[ 13 ] * pA[ 170 ] + pA[ 19 ] * pA[ 171 ] + pA[ 25 ] * pA[ 172 ] + pA[ 31 ] * pA[ 173 ];
        pA[ 314 ] = pA[ 2 ] * pA[ 168 ] + pA[ 8 ] * pA[ 169 ] + pA[ 14 ] * pA[ 170 ] + pA[ 20 ] * pA[ 171 ] + pA[ 26 ] * pA[ 172 ] + pA[ 32 ] * pA[ 173 ];
        pA[ 315 ] = pA[ 3 ] * pA[ 168 ] + pA[ 9 ] * pA[ 169 ] + pA[ 15 ] * pA[ 170 ] + pA[ 21 ] * pA[ 171 ] + pA[ 27 ] * pA[ 172 ] + pA[ 33 ] * pA[ 173 ];
        pA[ 316 ] = pA[ 4 ] * pA[ 168 ] + pA[ 10 ] * pA[ 169 ] + pA[ 16 ] * pA[ 170 ] + pA[ 22 ] * pA[ 171 ] + pA[ 28 ] * pA[ 172 ] + pA[ 34 ] * pA[ 173 ];
        pA[ 317 ] = pA[ 5 ] * pA[ 168 ] + pA[ 11 ] * pA[ 169 ] + pA[ 17 ] * pA[ 170 ] + pA[ 23 ] * pA[ 171 ] + pA[ 29 ] * pA[ 172 ] + pA[ 35 ] * pA[ 173 ];
    }
    {
        pA[ 318 ] = pA[ 0 ] * pA[ 174 ] + pA[ 6 ] * pA[ 175 ] + pA[ 12 ] * pA[ 176 ] + pA[ 18 ] * pA[ 177 ] + pA[ 24 ] * pA[ 178 ] + pA[ 30 ] * pA[ 179 ];
        pA[ 319 ] = pA[ 1 ] * pA[ 174 ] + pA[ 7 ] * pA[ 175 ] + pA[ 13 ] * pA[ 176 ] + pA[ 19 ] * pA[ 177 ] + pA[ 25 ] * pA[ 178 ] + pA[ 31 ] * pA[ 179 ];
        pA[ 320 ] = pA[ 2 ] * pA[ 174 ] + pA[ 8 ] * pA[ 175 ] + pA[ 14 ] * pA[ 176 ] + pA[ 20 ] * pA[ 177 ] + pA[ 26 ] * pA[ 178 ] + pA[ 32 ] * pA[ 179 ];
        pA[ 321 ] = pA[ 3 ] * pA[ 174 ] + pA[ 9 ] * pA[ 175 ] + pA[ 15 ] * pA[ 176 ] + pA[ 21 ] * pA[ 177 ] + pA[ 27 ] * pA[ 178 ] + pA[ 33 ] * pA[ 179 ];
        pA[ 322 ] = pA[ 4 ] * pA[ 174 ] + pA[ 10 ] * pA[ 175 ] + pA[ 16 ] * pA[ 176 ] + pA[ 22 ] * pA[ 177 ] + pA[ 28 ] * pA[ 178 ] + pA[ 34 ] * pA[ 179 ];
        pA[ 323 ] = pA[ 5 ] * pA[ 174 ] + pA[ 11 ] * pA[ 175 ] + pA[ 17 ] * pA[ 176 ] + pA[ 23 ] * pA[ 177 ] + pA[ 29 ] * pA[ 178 ] + pA[ 35 ] * pA[ 179 ];
    }
    return  true;
}

(b)の行列掛け算の実装
inline bool mmultABX_b( const double* pA, const double* pB, double* pX )
{
    {
        pX[ 0 ] = pA[ 0 ] * pB[ 0 ] + pA[ 6 ] * pB[ 1 ] + pA[ 12 ] * pB[ 2 ] + pA[ 18 ] * pB[ 3 ] + pA[ 24 ] * pB[ 4 ] + pA[ 30 ] * pB[ 5 ];
        pX[ 1 ] = pA[ 1 ] * pB[ 0 ] + pA[ 7 ] * pB[ 1 ] + pA[ 13 ] * pB[ 2 ] + pA[ 19 ] * pB[ 3 ] + pA[ 25 ] * pB[ 4 ] + pA[ 31 ] * pB[ 5 ];
        pX[ 2 ] = pA[ 2 ] * pB[ 0 ] + pA[ 8 ] * pB[ 1 ] + pA[ 14 ] * pB[ 2 ] + pA[ 20 ] * pB[ 3 ] + pA[ 26 ] * pB[ 4 ] + pA[ 32 ] * pB[ 5 ];
        pX[ 3 ] = pA[ 3 ] * pB[ 0 ] + pA[ 9 ] * pB[ 1 ] + pA[ 15 ] * pB[ 2 ] + pA[ 21 ] * pB[ 3 ] + pA[ 27 ] * pB[ 4 ] + pA[ 33 ] * pB[ 5 ];
        pX[ 4 ] = pA[ 4 ] * pB[ 0 ] + pA[ 10 ] * pB[ 1 ] + pA[ 16 ] * pB[ 2 ] + pA[ 22 ] * pB[ 3 ] + pA[ 28 ] * pB[ 4 ] + pA[ 34 ] * pB[ 5 ];
        pX[ 5 ] = pA[ 5 ] * pB[ 0 ] + pA[ 11 ] * pB[ 1 ] + pA[ 17 ] * pB[ 2 ] + pA[ 23 ] * pB[ 3 ] + pA[ 29 ] * pB[ 4 ] + pA[ 35 ] * pB[ 5 ];
    }
    {
        pX[ 6 ] = pA[ 0 ] * pB[ 6 ] + pA[ 6 ] * pB[ 7 ] + pA[ 12 ] * pB[ 8 ] + pA[ 18 ] * pB[ 9 ] + pA[ 24 ] * pB[ 10 ] + pA[ 30 ] * pB[ 11 ];
        pX[ 7 ] = pA[ 1 ] * pB[ 6 ] + pA[ 7 ] * pB[ 7 ] + pA[ 13 ] * pB[ 8 ] + pA[ 19 ] * pB[ 9 ] + pA[ 25 ] * pB[ 10 ] + pA[ 31 ] * pB[ 11 ];
        pX[ 8 ] = pA[ 2 ] * pB[ 6 ] + pA[ 8 ] * pB[ 7 ] + pA[ 14 ] * pB[ 8 ] + pA[ 20 ] * pB[ 9 ] + pA[ 26 ] * pB[ 10 ] + pA[ 32 ] * pB[ 11 ];
        pX[ 9 ] = pA[ 3 ] * pB[ 6 ] + pA[ 9 ] * pB[ 7 ] + pA[ 15 ] * pB[ 8 ] + pA[ 21 ] * pB[ 9 ] + pA[ 27 ] * pB[ 10 ] + pA[ 33 ] * pB[ 11 ];
        pX[ 10 ] = pA[ 4 ] * pB[ 6 ] + pA[ 10 ] * pB[ 7 ] + pA[ 16 ] * pB[ 8 ] + pA[ 22 ] * pB[ 9 ] + pA[ 28 ] * pB[ 10 ] + pA[ 34 ] * pB[ 11 ];
        pX[ 11 ] = pA[ 5 ] * pB[ 6 ] + pA[ 11 ] * pB[ 7 ] + pA[ 17 ] * pB[ 8 ] + pA[ 23 ] * pB[ 9 ] + pA[ 29 ] * pB[ 10 ] + pA[ 35 ] * pB[ 11 ];
    }
    {
        pX[ 12 ] = pA[ 0 ] * pB[ 12 ] + pA[ 6 ] * pB[ 13 ] + pA[ 12 ] * pB[ 14 ] + pA[ 18 ] * pB[ 15 ] + pA[ 24 ] * pB[ 16 ] + pA[ 30 ] * pB[ 17 ];
        pX[ 13 ] = pA[ 1 ] * pB[ 12 ] + pA[ 7 ] * pB[ 13 ] + pA[ 13 ] * pB[ 14 ] + pA[ 19 ] * pB[ 15 ] + pA[ 25 ] * pB[ 16 ] + pA[ 31 ] * pB[ 17 ];
        pX[ 14 ] = pA[ 2 ] * pB[ 12 ] + pA[ 8 ] * pB[ 13 ] + pA[ 14 ] * pB[ 14 ] + pA[ 20 ] * pB[ 15 ] + pA[ 26 ] * pB[ 16 ] + pA[ 32 ] * pB[ 17 ];
        pX[ 15 ] = pA[ 3 ] * pB[ 12 ] + pA[ 9 ] * pB[ 13 ] + pA[ 15 ] * pB[ 14 ] + pA[ 21 ] * pB[ 15 ] + pA[ 27 ] * pB[ 16 ] + pA[ 33 ] * pB[ 17 ];
        pX[ 16 ] = pA[ 4 ] * pB[ 12 ] + pA[ 10 ] * pB[ 13 ] + pA[ 16 ] * pB[ 14 ] + pA[ 22 ] * pB[ 15 ] + pA[ 28 ] * pB[ 16 ] + pA[ 34 ] * pB[ 17 ];
        pX[ 17 ] = pA[ 5 ] * pB[ 12 ] + pA[ 11 ] * pB[ 13 ] + pA[ 17 ] * pB[ 14 ] + pA[ 23 ] * pB[ 15 ] + pA[ 29 ] * pB[ 16 ] + pA[ 35 ] * pB[ 17 ];
    }
    {
        pX[ 18 ] = pA[ 0 ] * pB[ 18 ] + pA[ 6 ] * pB[ 19 ] + pA[ 12 ] * pB[ 20 ] + pA[ 18 ] * pB[ 21 ] + pA[ 24 ] * pB[ 22 ] + pA[ 30 ] * pB[ 23 ];
        pX[ 19 ] = pA[ 1 ] * pB[ 18 ] + pA[ 7 ] * pB[ 19 ] + pA[ 13 ] * pB[ 20 ] + pA[ 19 ] * pB[ 21 ] + pA[ 25 ] * pB[ 22 ] + pA[ 31 ] * pB[ 23 ];
        pX[ 20 ] = pA[ 2 ] * pB[ 18 ] + pA[ 8 ] * pB[ 19 ] + pA[ 14 ] * pB[ 20 ] + pA[ 20 ] * pB[ 21 ] + pA[ 26 ] * pB[ 22 ] + pA[ 32 ] * pB[ 23 ];
        pX[ 21 ] = pA[ 3 ] * pB[ 18 ] + pA[ 9 ] * pB[ 19 ] + pA[ 15 ] * pB[ 20 ] + pA[ 21 ] * pB[ 21 ] + pA[ 27 ] * pB[ 22 ] + pA[ 33 ] * pB[ 23 ];
        pX[ 22 ] = pA[ 4 ] * pB[ 18 ] + pA[ 10 ] * pB[ 19 ] + pA[ 16 ] * pB[ 20 ] + pA[ 22 ] * pB[ 21 ] + pA[ 28 ] * pB[ 22 ] + pA[ 34 ] * pB[ 23 ];
        pX[ 23 ] = pA[ 5 ] * pB[ 18 ] + pA[ 11 ] * pB[ 19 ] + pA[ 17 ] * pB[ 20 ] + pA[ 23 ] * pB[ 21 ] + pA[ 29 ] * pB[ 22 ] + pA[ 35 ] * pB[ 23 ];
    }
    {
        pX[ 24 ] = pA[ 0 ] * pB[ 24 ] + pA[ 6 ] * pB[ 25 ] + pA[ 12 ] * pB[ 26 ] + pA[ 18 ] * pB[ 27 ] + pA[ 24 ] * pB[ 28 ] + pA[ 30 ] * pB[ 29 ];
        pX[ 25 ] = pA[ 1 ] * pB[ 24 ] + pA[ 7 ] * pB[ 25 ] + pA[ 13 ] * pB[ 26 ] + pA[ 19 ] * pB[ 27 ] + pA[ 25 ] * pB[ 28 ] + pA[ 31 ] * pB[ 29 ];
        pX[ 26 ] = pA[ 2 ] * pB[ 24 ] + pA[ 8 ] * pB[ 25 ] + pA[ 14 ] * pB[ 26 ] + pA[ 20 ] * pB[ 27 ] + pA[ 26 ] * pB[ 28 ] + pA[ 32 ] * pB[ 29 ];
        pX[ 27 ] = pA[ 3 ] * pB[ 24 ] + pA[ 9 ] * pB[ 25 ] + pA[ 15 ] * pB[ 26 ] + pA[ 21 ] * pB[ 27 ] + pA[ 27 ] * pB[ 28 ] + pA[ 33 ] * pB[ 29 ];
        pX[ 28 ] = pA[ 4 ] * pB[ 24 ] + pA[ 10 ] * pB[ 25 ] + pA[ 16 ] * pB[ 26 ] + pA[ 22 ] * pB[ 27 ] + pA[ 28 ] * pB[ 28 ] + pA[ 34 ] * pB[ 29 ];
        pX[ 29 ] = pA[ 5 ] * pB[ 24 ] + pA[ 11 ] * pB[ 25 ] + pA[ 17 ] * pB[ 26 ] + pA[ 23 ] * pB[ 27 ] + pA[ 29 ] * pB[ 28 ] + pA[ 35 ] * pB[ 29 ];
    }
    {
        pX[ 30 ] = pA[ 0 ] * pB[ 30 ] + pA[ 6 ] * pB[ 31 ] + pA[ 12 ] * pB[ 32 ] + pA[ 18 ] * pB[ 33 ] + pA[ 24 ] * pB[ 34 ] + pA[ 30 ] * pB[ 35 ];
        pX[ 31 ] = pA[ 1 ] * pB[ 30 ] + pA[ 7 ] * pB[ 31 ] + pA[ 13 ] * pB[ 32 ] + pA[ 19 ] * pB[ 33 ] + pA[ 25 ] * pB[ 34 ] + pA[ 31 ] * pB[ 35 ];
        pX[ 32 ] = pA[ 2 ] * pB[ 30 ] + pA[ 8 ] * pB[ 31 ] + pA[ 14 ] * pB[ 32 ] + pA[ 20 ] * pB[ 33 ] + pA[ 26 ] * pB[ 34 ] + pA[ 32 ] * pB[ 35 ];
        pX[ 33 ] = pA[ 3 ] * pB[ 30 ] + pA[ 9 ] * pB[ 31 ] + pA[ 15 ] * pB[ 32 ] + pA[ 21 ] * pB[ 33 ] + pA[ 27 ] * pB[ 34 ] + pA[ 33 ] * pB[ 35 ];
        pX[ 34 ] = pA[ 4 ] * pB[ 30 ] + pA[ 10 ] * pB[ 31 ] + pA[ 16 ] * pB[ 32 ] + pA[ 22 ] * pB[ 33 ] + pA[ 28 ] * pB[ 34 ] + pA[ 34 ] * pB[ 35 ];
        pX[ 35 ] = pA[ 5 ] * pB[ 30 ] + pA[ 11 ] * pB[ 31 ] + pA[ 17 ] * pB[ 32 ] + pA[ 23 ] * pB[ 33 ] + pA[ 29 ] * pB[ 34 ] + pA[ 35 ] * pB[ 35 ];
    }
    {
        pX[ 36 ] = pA[ 0 ] * pB[ 36 ] + pA[ 6 ] * pB[ 37 ] + pA[ 12 ] * pB[ 38 ] + pA[ 18 ] * pB[ 39 ] + pA[ 24 ] * pB[ 40 ] + pA[ 30 ] * pB[ 41 ];
        pX[ 37 ] = pA[ 1 ] * pB[ 36 ] + pA[ 7 ] * pB[ 37 ] + pA[ 13 ] * pB[ 38 ] + pA[ 19 ] * pB[ 39 ] + pA[ 25 ] * pB[ 40 ] + pA[ 31 ] * pB[ 41 ];
        pX[ 38 ] = pA[ 2 ] * pB[ 36 ] + pA[ 8 ] * pB[ 37 ] + pA[ 14 ] * pB[ 38 ] + pA[ 20 ] * pB[ 39 ] + pA[ 26 ] * pB[ 40 ] + pA[ 32 ] * pB[ 41 ];
        pX[ 39 ] = pA[ 3 ] * pB[ 36 ] + pA[ 9 ] * pB[ 37 ] + pA[ 15 ] * pB[ 38 ] + pA[ 21 ] * pB[ 39 ] + pA[ 27 ] * pB[ 40 ] + pA[ 33 ] * pB[ 41 ];
        pX[ 40 ] = pA[ 4 ] * pB[ 36 ] + pA[ 10 ] * pB[ 37 ] + pA[ 16 ] * pB[ 38 ] + pA[ 22 ] * pB[ 39 ] + pA[ 28 ] * pB[ 40 ] + pA[ 34 ] * pB[ 41 ];
        pX[ 41 ] = pA[ 5 ] * pB[ 36 ] + pA[ 11 ] * pB[ 37 ] + pA[ 17 ] * pB[ 38 ] + pA[ 23 ] * pB[ 39 ] + pA[ 29 ] * pB[ 40 ] + pA[ 35 ] * pB[ 41 ];
    }
    {
        pX[ 42 ] = pA[ 0 ] * pB[ 42 ] + pA[ 6 ] * pB[ 43 ] + pA[ 12 ] * pB[ 44 ] + pA[ 18 ] * pB[ 45 ] + pA[ 24 ] * pB[ 46 ] + pA[ 30 ] * pB[ 47 ];
        pX[ 43 ] = pA[ 1 ] * pB[ 42 ] + pA[ 7 ] * pB[ 43 ] + pA[ 13 ] * pB[ 44 ] + pA[ 19 ] * pB[ 45 ] + pA[ 25 ] * pB[ 46 ] + pA[ 31 ] * pB[ 47 ];
        pX[ 44 ] = pA[ 2 ] * pB[ 42 ] + pA[ 8 ] * pB[ 43 ] + pA[ 14 ] * pB[ 44 ] + pA[ 20 ] * pB[ 45 ] + pA[ 26 ] * pB[ 46 ] + pA[ 32 ] * pB[ 47 ];
        pX[ 45 ] = pA[ 3 ] * pB[ 42 ] + pA[ 9 ] * pB[ 43 ] + pA[ 15 ] * pB[ 44 ] + pA[ 21 ] * pB[ 45 ] + pA[ 27 ] * pB[ 46 ] + pA[ 33 ] * pB[ 47 ];
        pX[ 46 ] = pA[ 4 ] * pB[ 42 ] + pA[ 10 ] * pB[ 43 ] + pA[ 16 ] * pB[ 44 ] + pA[ 22 ] * pB[ 45 ] + pA[ 28 ] * pB[ 46 ] + pA[ 34 ] * pB[ 47 ];
        pX[ 47 ] = pA[ 5 ] * pB[ 42 ] + pA[ 11 ] * pB[ 43 ] + pA[ 17 ] * pB[ 44 ] + pA[ 23 ] * pB[ 45 ] + pA[ 29 ] * pB[ 46 ] + pA[ 35 ] * pB[ 47 ];
    }
    {
        pX[ 48 ] = pA[ 0 ] * pB[ 48 ] + pA[ 6 ] * pB[ 49 ] + pA[ 12 ] * pB[ 50 ] + pA[ 18 ] * pB[ 51 ] + pA[ 24 ] * pB[ 52 ] + pA[ 30 ] * pB[ 53 ];
        pX[ 49 ] = pA[ 1 ] * pB[ 48 ] + pA[ 7 ] * pB[ 49 ] + pA[ 13 ] * pB[ 50 ] + pA[ 19 ] * pB[ 51 ] + pA[ 25 ] * pB[ 52 ] + pA[ 31 ] * pB[ 53 ];
        pX[ 50 ] = pA[ 2 ] * pB[ 48 ] + pA[ 8 ] * pB[ 49 ] + pA[ 14 ] * pB[ 50 ] + pA[ 20 ] * pB[ 51 ] + pA[ 26 ] * pB[ 52 ] + pA[ 32 ] * pB[ 53 ];
        pX[ 51 ] = pA[ 3 ] * pB[ 48 ] + pA[ 9 ] * pB[ 49 ] + pA[ 15 ] * pB[ 50 ] + pA[ 21 ] * pB[ 51 ] + pA[ 27 ] * pB[ 52 ] + pA[ 33 ] * pB[ 53 ];
        pX[ 52 ] = pA[ 4 ] * pB[ 48 ] + pA[ 10 ] * pB[ 49 ] + pA[ 16 ] * pB[ 50 ] + pA[ 22 ] * pB[ 51 ] + pA[ 28 ] * pB[ 52 ] + pA[ 34 ] * pB[ 53 ];
        pX[ 53 ] = pA[ 5 ] * pB[ 48 ] + pA[ 11 ] * pB[ 49 ] + pA[ 17 ] * pB[ 50 ] + pA[ 23 ] * pB[ 51 ] + pA[ 29 ] * pB[ 52 ] + pA[ 35 ] * pB[ 53 ];
    }
    {
        pX[ 54 ] = pA[ 0 ] * pB[ 54 ] + pA[ 6 ] * pB[ 55 ] + pA[ 12 ] * pB[ 56 ] + pA[ 18 ] * pB[ 57 ] + pA[ 24 ] * pB[ 58 ] + pA[ 30 ] * pB[ 59 ];
        pX[ 55 ] = pA[ 1 ] * pB[ 54 ] + pA[ 7 ] * pB[ 55 ] + pA[ 13 ] * pB[ 56 ] + pA[ 19 ] * pB[ 57 ] + pA[ 25 ] * pB[ 58 ] + pA[ 31 ] * pB[ 59 ];
        pX[ 56 ] = pA[ 2 ] * pB[ 54 ] + pA[ 8 ] * pB[ 55 ] + pA[ 14 ] * pB[ 56 ] + pA[ 20 ] * pB[ 57 ] + pA[ 26 ] * pB[ 58 ] + pA[ 32 ] * pB[ 59 ];
        pX[ 57 ] = pA[ 3 ] * pB[ 54 ] + pA[ 9 ] * pB[ 55 ] + pA[ 15 ] * pB[ 56 ] + pA[ 21 ] * pB[ 57 ] + pA[ 27 ] * pB[ 58 ] + pA[ 33 ] * pB[ 59 ];
        pX[ 58 ] = pA[ 4 ] * pB[ 54 ] + pA[ 10 ] * pB[ 55 ] + pA[ 16 ] * pB[ 56 ] + pA[ 22 ] * pB[ 57 ] + pA[ 28 ] * pB[ 58 ] + pA[ 34 ] * pB[ 59 ];
        pX[ 59 ] = pA[ 5 ] * pB[ 54 ] + pA[ 11 ] * pB[ 55 ] + pA[ 17 ] * pB[ 56 ] + pA[ 23 ] * pB[ 57 ] + pA[ 29 ] * pB[ 58 ] + pA[ 35 ] * pB[ 59 ];
    }
    {
        pX[ 60 ] = pA[ 0 ] * pB[ 60 ] + pA[ 6 ] * pB[ 61 ] + pA[ 12 ] * pB[ 62 ] + pA[ 18 ] * pB[ 63 ] + pA[ 24 ] * pB[ 64 ] + pA[ 30 ] * pB[ 65 ];
        pX[ 61 ] = pA[ 1 ] * pB[ 60 ] + pA[ 7 ] * pB[ 61 ] + pA[ 13 ] * pB[ 62 ] + pA[ 19 ] * pB[ 63 ] + pA[ 25 ] * pB[ 64 ] + pA[ 31 ] * pB[ 65 ];
        pX[ 62 ] = pA[ 2 ] * pB[ 60 ] + pA[ 8 ] * pB[ 61 ] + pA[ 14 ] * pB[ 62 ] + pA[ 20 ] * pB[ 63 ] + pA[ 26 ] * pB[ 64 ] + pA[ 32 ] * pB[ 65 ];
        pX[ 63 ] = pA[ 3 ] * pB[ 60 ] + pA[ 9 ] * pB[ 61 ] + pA[ 15 ] * pB[ 62 ] + pA[ 21 ] * pB[ 63 ] + pA[ 27 ] * pB[ 64 ] + pA[ 33 ] * pB[ 65 ];
        pX[ 64 ] = pA[ 4 ] * pB[ 60 ] + pA[ 10 ] * pB[ 61 ] + pA[ 16 ] * pB[ 62 ] + pA[ 22 ] * pB[ 63 ] + pA[ 28 ] * pB[ 64 ] + pA[ 34 ] * pB[ 65 ];
        pX[ 65 ] = pA[ 5 ] * pB[ 60 ] + pA[ 11 ] * pB[ 61 ] + pA[ 17 ] * pB[ 62 ] + pA[ 23 ] * pB[ 63 ] + pA[ 29 ] * pB[ 64 ] + pA[ 35 ] * pB[ 65 ];
    }
    {
        pX[ 66 ] = pA[ 0 ] * pB[ 66 ] + pA[ 6 ] * pB[ 67 ] + pA[ 12 ] * pB[ 68 ] + pA[ 18 ] * pB[ 69 ] + pA[ 24 ] * pB[ 70 ] + pA[ 30 ] * pB[ 71 ];
        pX[ 67 ] = pA[ 1 ] * pB[ 66 ] + pA[ 7 ] * pB[ 67 ] + pA[ 13 ] * pB[ 68 ] + pA[ 19 ] * pB[ 69 ] + pA[ 25 ] * pB[ 70 ] + pA[ 31 ] * pB[ 71 ];
        pX[ 68 ] = pA[ 2 ] * pB[ 66 ] + pA[ 8 ] * pB[ 67 ] + pA[ 14 ] * pB[ 68 ] + pA[ 20 ] * pB[ 69 ] + pA[ 26 ] * pB[ 70 ] + pA[ 32 ] * pB[ 71 ];
        pX[ 69 ] = pA[ 3 ] * pB[ 66 ] + pA[ 9 ] * pB[ 67 ] + pA[ 15 ] * pB[ 68 ] + pA[ 21 ] * pB[ 69 ] + pA[ 27 ] * pB[ 70 ] + pA[ 33 ] * pB[ 71 ];
        pX[ 70 ] = pA[ 4 ] * pB[ 66 ] + pA[ 10 ] * pB[ 67 ] + pA[ 16 ] * pB[ 68 ] + pA[ 22 ] * pB[ 69 ] + pA[ 28 ] * pB[ 70 ] + pA[ 34 ] * pB[ 71 ];
        pX[ 71 ] = pA[ 5 ] * pB[ 66 ] + pA[ 11 ] * pB[ 67 ] + pA[ 17 ] * pB[ 68 ] + pA[ 23 ] * pB[ 69 ] + pA[ 29 ] * pB[ 70 ] + pA[ 35 ] * pB[ 71 ];
    }
    {
        pX[ 72 ] = pA[ 0 ] * pB[ 72 ] + pA[ 6 ] * pB[ 73 ] + pA[ 12 ] * pB[ 74 ] + pA[ 18 ] * pB[ 75 ] + pA[ 24 ] * pB[ 76 ] + pA[ 30 ] * pB[ 77 ];
        pX[ 73 ] = pA[ 1 ] * pB[ 72 ] + pA[ 7 ] * pB[ 73 ] + pA[ 13 ] * pB[ 74 ] + pA[ 19 ] * pB[ 75 ] + pA[ 25 ] * pB[ 76 ] + pA[ 31 ] * pB[ 77 ];
        pX[ 74 ] = pA[ 2 ] * pB[ 72 ] + pA[ 8 ] * pB[ 73 ] + pA[ 14 ] * pB[ 74 ] + pA[ 20 ] * pB[ 75 ] + pA[ 26 ] * pB[ 76 ] + pA[ 32 ] * pB[ 77 ];
        pX[ 75 ] = pA[ 3 ] * pB[ 72 ] + pA[ 9 ] * pB[ 73 ] + pA[ 15 ] * pB[ 74 ] + pA[ 21 ] * pB[ 75 ] + pA[ 27 ] * pB[ 76 ] + pA[ 33 ] * pB[ 77 ];
        pX[ 76 ] = pA[ 4 ] * pB[ 72 ] + pA[ 10 ] * pB[ 73 ] + pA[ 16 ] * pB[ 74 ] + pA[ 22 ] * pB[ 75 ] + pA[ 28 ] * pB[ 76 ] + pA[ 34 ] * pB[ 77 ];
        pX[ 77 ] = pA[ 5 ] * pB[ 72 ] + pA[ 11 ] * pB[ 73 ] + pA[ 17 ] * pB[ 74 ] + pA[ 23 ] * pB[ 75 ] + pA[ 29 ] * pB[ 76 ] + pA[ 35 ] * pB[ 77 ];
    }
    {
        pX[ 78 ] = pA[ 0 ] * pB[ 78 ] + pA[ 6 ] * pB[ 79 ] + pA[ 12 ] * pB[ 80 ] + pA[ 18 ] * pB[ 81 ] + pA[ 24 ] * pB[ 82 ] + pA[ 30 ] * pB[ 83 ];
        pX[ 79 ] = pA[ 1 ] * pB[ 78 ] + pA[ 7 ] * pB[ 79 ] + pA[ 13 ] * pB[ 80 ] + pA[ 19 ] * pB[ 81 ] + pA[ 25 ] * pB[ 82 ] + pA[ 31 ] * pB[ 83 ];
        pX[ 80 ] = pA[ 2 ] * pB[ 78 ] + pA[ 8 ] * pB[ 79 ] + pA[ 14 ] * pB[ 80 ] + pA[ 20 ] * pB[ 81 ] + pA[ 26 ] * pB[ 82 ] + pA[ 32 ] * pB[ 83 ];
        pX[ 81 ] = pA[ 3 ] * pB[ 78 ] + pA[ 9 ] * pB[ 79 ] + pA[ 15 ] * pB[ 80 ] + pA[ 21 ] * pB[ 81 ] + pA[ 27 ] * pB[ 82 ] + pA[ 33 ] * pB[ 83 ];
        pX[ 82 ] = pA[ 4 ] * pB[ 78 ] + pA[ 10 ] * pB[ 79 ] + pA[ 16 ] * pB[ 80 ] + pA[ 22 ] * pB[ 81 ] + pA[ 28 ] * pB[ 82 ] + pA[ 34 ] * pB[ 83 ];
        pX[ 83 ] = pA[ 5 ] * pB[ 78 ] + pA[ 11 ] * pB[ 79 ] + pA[ 17 ] * pB[ 80 ] + pA[ 23 ] * pB[ 81 ] + pA[ 29 ] * pB[ 82 ] + pA[ 35 ] * pB[ 83 ];
    }
    {
        pX[ 84 ] = pA[ 0 ] * pB[ 84 ] + pA[ 6 ] * pB[ 85 ] + pA[ 12 ] * pB[ 86 ] + pA[ 18 ] * pB[ 87 ] + pA[ 24 ] * pB[ 88 ] + pA[ 30 ] * pB[ 89 ];
        pX[ 85 ] = pA[ 1 ] * pB[ 84 ] + pA[ 7 ] * pB[ 85 ] + pA[ 13 ] * pB[ 86 ] + pA[ 19 ] * pB[ 87 ] + pA[ 25 ] * pB[ 88 ] + pA[ 31 ] * pB[ 89 ];
        pX[ 86 ] = pA[ 2 ] * pB[ 84 ] + pA[ 8 ] * pB[ 85 ] + pA[ 14 ] * pB[ 86 ] + pA[ 20 ] * pB[ 87 ] + pA[ 26 ] * pB[ 88 ] + pA[ 32 ] * pB[ 89 ];
        pX[ 87 ] = pA[ 3 ] * pB[ 84 ] + pA[ 9 ] * pB[ 85 ] + pA[ 15 ] * pB[ 86 ] + pA[ 21 ] * pB[ 87 ] + pA[ 27 ] * pB[ 88 ] + pA[ 33 ] * pB[ 89 ];
        pX[ 88 ] = pA[ 4 ] * pB[ 84 ] + pA[ 10 ] * pB[ 85 ] + pA[ 16 ] * pB[ 86 ] + pA[ 22 ] * pB[ 87 ] + pA[ 28 ] * pB[ 88 ] + pA[ 34 ] * pB[ 89 ];
        pX[ 89 ] = pA[ 5 ] * pB[ 84 ] + pA[ 11 ] * pB[ 85 ] + pA[ 17 ] * pB[ 86 ] + pA[ 23 ] * pB[ 87 ] + pA[ 29 ] * pB[ 88 ] + pA[ 35 ] * pB[ 89 ];
    }
    {
        pX[ 90 ] = pA[ 0 ] * pB[ 90 ] + pA[ 6 ] * pB[ 91 ] + pA[ 12 ] * pB[ 92 ] + pA[ 18 ] * pB[ 93 ] + pA[ 24 ] * pB[ 94 ] + pA[ 30 ] * pB[ 95 ];
        pX[ 91 ] = pA[ 1 ] * pB[ 90 ] + pA[ 7 ] * pB[ 91 ] + pA[ 13 ] * pB[ 92 ] + pA[ 19 ] * pB[ 93 ] + pA[ 25 ] * pB[ 94 ] + pA[ 31 ] * pB[ 95 ];
        pX[ 92 ] = pA[ 2 ] * pB[ 90 ] + pA[ 8 ] * pB[ 91 ] + pA[ 14 ] * pB[ 92 ] + pA[ 20 ] * pB[ 93 ] + pA[ 26 ] * pB[ 94 ] + pA[ 32 ] * pB[ 95 ];
        pX[ 93 ] = pA[ 3 ] * pB[ 90 ] + pA[ 9 ] * pB[ 91 ] + pA[ 15 ] * pB[ 92 ] + pA[ 21 ] * pB[ 93 ] + pA[ 27 ] * pB[ 94 ] + pA[ 33 ] * pB[ 95 ];
        pX[ 94 ] = pA[ 4 ] * pB[ 90 ] + pA[ 10 ] * pB[ 91 ] + pA[ 16 ] * pB[ 92 ] + pA[ 22 ] * pB[ 93 ] + pA[ 28 ] * pB[ 94 ] + pA[ 34 ] * pB[ 95 ];
        pX[ 95 ] = pA[ 5 ] * pB[ 90 ] + pA[ 11 ] * pB[ 91 ] + pA[ 17 ] * pB[ 92 ] + pA[ 23 ] * pB[ 93 ] + pA[ 29 ] * pB[ 94 ] + pA[ 35 ] * pB[ 95 ];
    }
    {
        pX[ 96 ] = pA[ 0 ] * pB[ 96 ] + pA[ 6 ] * pB[ 97 ] + pA[ 12 ] * pB[ 98 ] + pA[ 18 ] * pB[ 99 ] + pA[ 24 ] * pB[ 100 ] + pA[ 30 ] * pB[ 101 ];
        pX[ 97 ] = pA[ 1 ] * pB[ 96 ] + pA[ 7 ] * pB[ 97 ] + pA[ 13 ] * pB[ 98 ] + pA[ 19 ] * pB[ 99 ] + pA[ 25 ] * pB[ 100 ] + pA[ 31 ] * pB[ 101 ];
        pX[ 98 ] = pA[ 2 ] * pB[ 96 ] + pA[ 8 ] * pB[ 97 ] + pA[ 14 ] * pB[ 98 ] + pA[ 20 ] * pB[ 99 ] + pA[ 26 ] * pB[ 100 ] + pA[ 32 ] * pB[ 101 ];
        pX[ 99 ] = pA[ 3 ] * pB[ 96 ] + pA[ 9 ] * pB[ 97 ] + pA[ 15 ] * pB[ 98 ] + pA[ 21 ] * pB[ 99 ] + pA[ 27 ] * pB[ 100 ] + pA[ 33 ] * pB[ 101 ];
        pX[ 100 ] = pA[ 4 ] * pB[ 96 ] + pA[ 10 ] * pB[ 97 ] + pA[ 16 ] * pB[ 98 ] + pA[ 22 ] * pB[ 99 ] + pA[ 28 ] * pB[ 100 ] + pA[ 34 ] * pB[ 101 ];
        pX[ 101 ] = pA[ 5 ] * pB[ 96 ] + pA[ 11 ] * pB[ 97 ] + pA[ 17 ] * pB[ 98 ] + pA[ 23 ] * pB[ 99 ] + pA[ 29 ] * pB[ 100 ] + pA[ 35 ] * pB[ 101 ];
    }
    {
        pX[ 102 ] = pA[ 0 ] * pB[ 102 ] + pA[ 6 ] * pB[ 103 ] + pA[ 12 ] * pB[ 104 ] + pA[ 18 ] * pB[ 105 ] + pA[ 24 ] * pB[ 106 ] + pA[ 30 ] * pB[ 107 ];
        pX[ 103 ] = pA[ 1 ] * pB[ 102 ] + pA[ 7 ] * pB[ 103 ] + pA[ 13 ] * pB[ 104 ] + pA[ 19 ] * pB[ 105 ] + pA[ 25 ] * pB[ 106 ] + pA[ 31 ] * pB[ 107 ];
        pX[ 104 ] = pA[ 2 ] * pB[ 102 ] + pA[ 8 ] * pB[ 103 ] + pA[ 14 ] * pB[ 104 ] + pA[ 20 ] * pB[ 105 ] + pA[ 26 ] * pB[ 106 ] + pA[ 32 ] * pB[ 107 ];
        pX[ 105 ] = pA[ 3 ] * pB[ 102 ] + pA[ 9 ] * pB[ 103 ] + pA[ 15 ] * pB[ 104 ] + pA[ 21 ] * pB[ 105 ] + pA[ 27 ] * pB[ 106 ] + pA[ 33 ] * pB[ 107 ];
        pX[ 106 ] = pA[ 4 ] * pB[ 102 ] + pA[ 10 ] * pB[ 103 ] + pA[ 16 ] * pB[ 104 ] + pA[ 22 ] * pB[ 105 ] + pA[ 28 ] * pB[ 106 ] + pA[ 34 ] * pB[ 107 ];
        pX[ 107 ] = pA[ 5 ] * pB[ 102 ] + pA[ 11 ] * pB[ 103 ] + pA[ 17 ] * pB[ 104 ] + pA[ 23 ] * pB[ 105 ] + pA[ 29 ] * pB[ 106 ] + pA[ 35 ] * pB[ 107 ];
    }
    {
        pX[ 108 ] = pA[ 0 ] * pB[ 108 ] + pA[ 6 ] * pB[ 109 ] + pA[ 12 ] * pB[ 110 ] + pA[ 18 ] * pB[ 111 ] + pA[ 24 ] * pB[ 112 ] + pA[ 30 ] * pB[ 113 ];
        pX[ 109 ] = pA[ 1 ] * pB[ 108 ] + pA[ 7 ] * pB[ 109 ] + pA[ 13 ] * pB[ 110 ] + pA[ 19 ] * pB[ 111 ] + pA[ 25 ] * pB[ 112 ] + pA[ 31 ] * pB[ 113 ];
        pX[ 110 ] = pA[ 2 ] * pB[ 108 ] + pA[ 8 ] * pB[ 109 ] + pA[ 14 ] * pB[ 110 ] + pA[ 20 ] * pB[ 111 ] + pA[ 26 ] * pB[ 112 ] + pA[ 32 ] * pB[ 113 ];
        pX[ 111 ] = pA[ 3 ] * pB[ 108 ] + pA[ 9 ] * pB[ 109 ] + pA[ 15 ] * pB[ 110 ] + pA[ 21 ] * pB[ 111 ] + pA[ 27 ] * pB[ 112 ] + pA[ 33 ] * pB[ 113 ];
        pX[ 112 ] = pA[ 4 ] * pB[ 108 ] + pA[ 10 ] * pB[ 109 ] + pA[ 16 ] * pB[ 110 ] + pA[ 22 ] * pB[ 111 ] + pA[ 28 ] * pB[ 112 ] + pA[ 34 ] * pB[ 113 ];
        pX[ 113 ] = pA[ 5 ] * pB[ 108 ] + pA[ 11 ] * pB[ 109 ] + pA[ 17 ] * pB[ 110 ] + pA[ 23 ] * pB[ 111 ] + pA[ 29 ] * pB[ 112 ] + pA[ 35 ] * pB[ 113 ];
    }
    {
        pX[ 114 ] = pA[ 0 ] * pB[ 114 ] + pA[ 6 ] * pB[ 115 ] + pA[ 12 ] * pB[ 116 ] + pA[ 18 ] * pB[ 117 ] + pA[ 24 ] * pB[ 118 ] + pA[ 30 ] * pB[ 119 ];
        pX[ 115 ] = pA[ 1 ] * pB[ 114 ] + pA[ 7 ] * pB[ 115 ] + pA[ 13 ] * pB[ 116 ] + pA[ 19 ] * pB[ 117 ] + pA[ 25 ] * pB[ 118 ] + pA[ 31 ] * pB[ 119 ];
        pX[ 116 ] = pA[ 2 ] * pB[ 114 ] + pA[ 8 ] * pB[ 115 ] + pA[ 14 ] * pB[ 116 ] + pA[ 20 ] * pB[ 117 ] + pA[ 26 ] * pB[ 118 ] + pA[ 32 ] * pB[ 119 ];
        pX[ 117 ] = pA[ 3 ] * pB[ 114 ] + pA[ 9 ] * pB[ 115 ] + pA[ 15 ] * pB[ 116 ] + pA[ 21 ] * pB[ 117 ] + pA[ 27 ] * pB[ 118 ] + pA[ 33 ] * pB[ 119 ];
        pX[ 118 ] = pA[ 4 ] * pB[ 114 ] + pA[ 10 ] * pB[ 115 ] + pA[ 16 ] * pB[ 116 ] + pA[ 22 ] * pB[ 117 ] + pA[ 28 ] * pB[ 118 ] + pA[ 34 ] * pB[ 119 ];
        pX[ 119 ] = pA[ 5 ] * pB[ 114 ] + pA[ 11 ] * pB[ 115 ] + pA[ 17 ] * pB[ 116 ] + pA[ 23 ] * pB[ 117 ] + pA[ 29 ] * pB[ 118 ] + pA[ 35 ] * pB[ 119 ];
    }
    {
        pX[ 120 ] = pA[ 0 ] * pB[ 120 ] + pA[ 6 ] * pB[ 121 ] + pA[ 12 ] * pB[ 122 ] + pA[ 18 ] * pB[ 123 ] + pA[ 24 ] * pB[ 124 ] + pA[ 30 ] * pB[ 125 ];
        pX[ 121 ] = pA[ 1 ] * pB[ 120 ] + pA[ 7 ] * pB[ 121 ] + pA[ 13 ] * pB[ 122 ] + pA[ 19 ] * pB[ 123 ] + pA[ 25 ] * pB[ 124 ] + pA[ 31 ] * pB[ 125 ];
        pX[ 122 ] = pA[ 2 ] * pB[ 120 ] + pA[ 8 ] * pB[ 121 ] + pA[ 14 ] * pB[ 122 ] + pA[ 20 ] * pB[ 123 ] + pA[ 26 ] * pB[ 124 ] + pA[ 32 ] * pB[ 125 ];
        pX[ 123 ] = pA[ 3 ] * pB[ 120 ] + pA[ 9 ] * pB[ 121 ] + pA[ 15 ] * pB[ 122 ] + pA[ 21 ] * pB[ 123 ] + pA[ 27 ] * pB[ 124 ] + pA[ 33 ] * pB[ 125 ];
        pX[ 124 ] = pA[ 4 ] * pB[ 120 ] + pA[ 10 ] * pB[ 121 ] + pA[ 16 ] * pB[ 122 ] + pA[ 22 ] * pB[ 123 ] + pA[ 28 ] * pB[ 124 ] + pA[ 34 ] * pB[ 125 ];
        pX[ 125 ] = pA[ 5 ] * pB[ 120 ] + pA[ 11 ] * pB[ 121 ] + pA[ 17 ] * pB[ 122 ] + pA[ 23 ] * pB[ 123 ] + pA[ 29 ] * pB[ 124 ] + pA[ 35 ] * pB[ 125 ];
    }
    {
        pX[ 126 ] = pA[ 0 ] * pB[ 126 ] + pA[ 6 ] * pB[ 127 ] + pA[ 12 ] * pB[ 128 ] + pA[ 18 ] * pB[ 129 ] + pA[ 24 ] * pB[ 130 ] + pA[ 30 ] * pB[ 131 ];
        pX[ 127 ] = pA[ 1 ] * pB[ 126 ] + pA[ 7 ] * pB[ 127 ] + pA[ 13 ] * pB[ 128 ] + pA[ 19 ] * pB[ 129 ] + pA[ 25 ] * pB[ 130 ] + pA[ 31 ] * pB[ 131 ];
        pX[ 128 ] = pA[ 2 ] * pB[ 126 ] + pA[ 8 ] * pB[ 127 ] + pA[ 14 ] * pB[ 128 ] + pA[ 20 ] * pB[ 129 ] + pA[ 26 ] * pB[ 130 ] + pA[ 32 ] * pB[ 131 ];
        pX[ 129 ] = pA[ 3 ] * pB[ 126 ] + pA[ 9 ] * pB[ 127 ] + pA[ 15 ] * pB[ 128 ] + pA[ 21 ] * pB[ 129 ] + pA[ 27 ] * pB[ 130 ] + pA[ 33 ] * pB[ 131 ];
        pX[ 130 ] = pA[ 4 ] * pB[ 126 ] + pA[ 10 ] * pB[ 127 ] + pA[ 16 ] * pB[ 128 ] + pA[ 22 ] * pB[ 129 ] + pA[ 28 ] * pB[ 130 ] + pA[ 34 ] * pB[ 131 ];
        pX[ 131 ] = pA[ 5 ] * pB[ 126 ] + pA[ 11 ] * pB[ 127 ] + pA[ 17 ] * pB[ 128 ] + pA[ 23 ] * pB[ 129 ] + pA[ 29 ] * pB[ 130 ] + pA[ 35 ] * pB[ 131 ];
    }
    {
        pX[ 132 ] = pA[ 0 ] * pB[ 132 ] + pA[ 6 ] * pB[ 133 ] + pA[ 12 ] * pB[ 134 ] + pA[ 18 ] * pB[ 135 ] + pA[ 24 ] * pB[ 136 ] + pA[ 30 ] * pB[ 137 ];
        pX[ 133 ] = pA[ 1 ] * pB[ 132 ] + pA[ 7 ] * pB[ 133 ] + pA[ 13 ] * pB[ 134 ] + pA[ 19 ] * pB[ 135 ] + pA[ 25 ] * pB[ 136 ] + pA[ 31 ] * pB[ 137 ];
        pX[ 134 ] = pA[ 2 ] * pB[ 132 ] + pA[ 8 ] * pB[ 133 ] + pA[ 14 ] * pB[ 134 ] + pA[ 20 ] * pB[ 135 ] + pA[ 26 ] * pB[ 136 ] + pA[ 32 ] * pB[ 137 ];
        pX[ 135 ] = pA[ 3 ] * pB[ 132 ] + pA[ 9 ] * pB[ 133 ] + pA[ 15 ] * pB[ 134 ] + pA[ 21 ] * pB[ 135 ] + pA[ 27 ] * pB[ 136 ] + pA[ 33 ] * pB[ 137 ];
        pX[ 136 ] = pA[ 4 ] * pB[ 132 ] + pA[ 10 ] * pB[ 133 ] + pA[ 16 ] * pB[ 134 ] + pA[ 22 ] * pB[ 135 ] + pA[ 28 ] * pB[ 136 ] + pA[ 34 ] * pB[ 137 ];
        pX[ 137 ] = pA[ 5 ] * pB[ 132 ] + pA[ 11 ] * pB[ 133 ] + pA[ 17 ] * pB[ 134 ] + pA[ 23 ] * pB[ 135 ] + pA[ 29 ] * pB[ 136 ] + pA[ 35 ] * pB[ 137 ];
    }
    {
        pX[ 138 ] = pA[ 0 ] * pB[ 138 ] + pA[ 6 ] * pB[ 139 ] + pA[ 12 ] * pB[ 140 ] + pA[ 18 ] * pB[ 141 ] + pA[ 24 ] * pB[ 142 ] + pA[ 30 ] * pB[ 143 ];
        pX[ 139 ] = pA[ 1 ] * pB[ 138 ] + pA[ 7 ] * pB[ 139 ] + pA[ 13 ] * pB[ 140 ] + pA[ 19 ] * pB[ 141 ] + pA[ 25 ] * pB[ 142 ] + pA[ 31 ] * pB[ 143 ];
        pX[ 140 ] = pA[ 2 ] * pB[ 138 ] + pA[ 8 ] * pB[ 139 ] + pA[ 14 ] * pB[ 140 ] + pA[ 20 ] * pB[ 141 ] + pA[ 26 ] * pB[ 142 ] + pA[ 32 ] * pB[ 143 ];
        pX[ 141 ] = pA[ 3 ] * pB[ 138 ] + pA[ 9 ] * pB[ 139 ] + pA[ 15 ] * pB[ 140 ] + pA[ 21 ] * pB[ 141 ] + pA[ 27 ] * pB[ 142 ] + pA[ 33 ] * pB[ 143 ];
        pX[ 142 ] = pA[ 4 ] * pB[ 138 ] + pA[ 10 ] * pB[ 139 ] + pA[ 16 ] * pB[ 140 ] + pA[ 22 ] * pB[ 141 ] + pA[ 28 ] * pB[ 142 ] + pA[ 34 ] * pB[ 143 ];
        pX[ 143 ] = pA[ 5 ] * pB[ 138 ] + pA[ 11 ] * pB[ 139 ] + pA[ 17 ] * pB[ 140 ] + pA[ 23 ] * pB[ 141 ] + pA[ 29 ] * pB[ 142 ] + pA[ 35 ] * pB[ 143 ];
    }
    return true;
}

これを100万回ループします。OpenMPで並列化したソースを、(a)(b)それぞれ以下に示します。ループの並列化は、parallel forではなく、並列スレッドを作成してから配列をmax_thread数で分割し、それぞれのスレッドでループしております。こちらの方がメモリに連続配置されたデータへのアクセスがし易いであろう、との目論見からです。

(a)の100万回ループ実装

// ==============================================================================
//  (a)
// ==============================================================================
bool    loop_mmult_a( double* prA , const size_t nElm , const size_t nElmAll ){
#pragma omp parallel
    {
        //  スレッド数とこのスレッドindexを取得
        const int nThreadMax = omp_get_max_threads();
        const int iThread    = omp_get_thread_num();

        //  1スレッド当たりで計算する要素数
        const int64_t   nElm_local = nElm / nThreadMax;

        //  このスレッドで計算する[A]行列の配列の先頭ポインタ
        double* prA_local = prA + ( nElm_local * iThread * nElmAll );

        //  各スレッドで分担して100万回計算します
        for( int64_t i = 0; i < nElm_local / 10; i++ ){
            //  10回アンローリング
            mmultABX_a( prA_local );    prA_local += nElmAll;
            mmultABX_a( prA_local );    prA_local += nElmAll;
            mmultABX_a( prA_local );    prA_local += nElmAll;
            mmultABX_a( prA_local );    prA_local += nElmAll;
            mmultABX_a( prA_local );    prA_local += nElmAll;
            mmultABX_a( prA_local );    prA_local += nElmAll;
            mmultABX_a( prA_local );    prA_local += nElmAll;
            mmultABX_a( prA_local );    prA_local += nElmAll;
            mmultABX_a( prA_local );    prA_local += nElmAll;
            mmultABX_a( prA_local );    prA_local += nElmAll;
        }
    }

    return  true;
}

(b)の100万回ループ実装

// ==============================================================================
//  (b)
// ==============================================================================
bool    loop_mmult_b(
    double*         prA,    double*         prB,    double*         prX,
    const size_t    nAelm,  const size_t    nBelm,  const size_t    nXelm, 
    const size_t    nElm,
    const size_t    nElmAll ){
#pragma omp parallel
    {
        //  スレッド数とこのスレッドindexを取得
        const int nThreadMax = omp_get_max_threads();
        const int iThread    = omp_get_thread_num();

        //  1スレッド当たりで計算する要素数
        const int64_t   nElm_local = nElm / nThreadMax;

        //  このスレッドで計算する[A][B][X]行列の配列のそれぞれの先頭ポインタ
        double* prA_local = prA + ( nElm_local * iThread * nAelm );
        double* prB_local = prB + ( nElm_local * iThread * nBelm );
        double* prX_local = prX + ( nElm_local * iThread * nXelm );

        //  各スレッドで分担して100万回計算します
        for( int64_t i = 0; i < nElm_local / 10; i++ ){
            //  10回アンローリング
            mmultABX_b( prA_local, prB_local, prX_local );  prA_local += nAelm; prB_local += nBelm; prX_local += nXelm;
            mmultABX_b( prA_local, prB_local, prX_local );  prA_local += nAelm; prB_local += nBelm; prX_local += nXelm;
            mmultABX_b( prA_local, prB_local, prX_local );  prA_local += nAelm; prB_local += nBelm; prX_local += nXelm;
            mmultABX_b( prA_local, prB_local, prX_local );  prA_local += nAelm; prB_local += nBelm; prX_local += nXelm;
            mmultABX_b( prA_local, prB_local, prX_local );  prA_local += nAelm; prB_local += nBelm; prX_local += nXelm;
            mmultABX_b( prA_local, prB_local, prX_local );  prA_local += nAelm; prB_local += nBelm; prX_local += nXelm;
            mmultABX_b( prA_local, prB_local, prX_local );  prA_local += nAelm; prB_local += nBelm; prX_local += nXelm;
            mmultABX_b( prA_local, prB_local, prX_local );  prA_local += nAelm; prB_local += nBelm; prX_local += nXelm;
            mmultABX_b( prA_local, prB_local, prX_local );  prA_local += nAelm; prB_local += nBelm; prX_local += nXelm;
            mmultABX_b( prA_local, prB_local, prX_local );  prA_local += nAelm; prB_local += nBelm; prX_local += nXelm;
        }
    }

    return  true;
}

計算結果

(a)(b)を比較した結果を示します。下記はOpenMPにより並列化した場合の結果です。(シリアル処理でも行ったのですが、(a)(b)に差は出ませんでした。)

(a) (b)
計算時間 (s) 計算速度(GFLOPS) 計算時間 (s) 計算速度(GFLOPS)
1回目 0.225 7.10 0.156 10.2
2回目 0.229 6.91 0.153 10.4
3回目 0.221 7.16 0.152 10.4
平均 0.224 7.06 0.153 10.3

※計算速度(GFLOPS)は、2つ行列[A][B]の積における浮動小数点演算回数(1584x100万回)を計算時間で割って求めました。

(a)と(b)の比較は、期待に反して、(b)の方が(a)よりも約45%も高速という結果となりました。10GFLOPSというのは速すぎる気もしますが…(~_~;)。(CPUベンチマーク[2]をご参照)

当方の環境は下記です。

CPU Core i5 3550 ( ivy bridge )
Memory 32 GB(DDR3 1333 MHz)
Compiler Visual C++ 2017

結論

 ともあれ、長く疑問に思っていた有限要素法における小行列の格納形式について、2つの手法を計測により比較した結果、A,B,Xを各個に連続してメモリに配置する(b)の方法が計算速度に有利、という結果が得られました。これを持って最終判断をするには至りませんが、参考値として活かせれば良いかと考えております。

参考文献
[1] Fortran90/95による有限要素法プログラミング, 藤井文夫,田中真人,佐藤維美著,丸善出版
[2] https://asteroidsathome.net/boinc/cpu_list.php

zmtx_pp
C++とFORTRANを使ってます。今更ですがBoostを勉強したほうが良さそうだと思ってます。記事の内容、間違ってたら是非ご指摘ください。
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away