C++を使っておりますが、テーマが大分それてしまいました。
カレンダーで、このタイトルだけ趣旨がずれて浮いてる気がしてます…失礼いたします…
#本稿でやったこと
大量の小規模行列の掛け算に対し、2種類のメモリ配置を試し、どちらが計算速度が速いかを確認しました。
いつも通り設計していたのですが、速度向上を図って試行錯誤していくうちに、C言語のようなソースコードになってしまいました。いずれ整えます。
私はHPCの分野は素人ですので、ご意見頂けると有り難いです。
#背景
※これはテーマからかなり逸れるので、飛ばしても構いません。
多くの工学分野の解析実務で、FEM（有限要素法）という数値解析手法が広く用いられています。FEMは振動解析や郷土解析などの構造計算の分野でも使われます。私は構造計算でFEMを使うことが多いため、今回それを前提とします。
FEMによる数値解析においては、限られたリソースで高速に解く必要があり、効率の良いメモリ配置が重要と思います。FEMの書籍には、解説とともにサンプルプログラムが公開されているものがあります（例えば[1]）。しかしその解説において、計算速度に対するメモリ配置の説明が無いと私は思っております。そこで、現状、私が想定している2種類のメモリ配置の方法について計算速度を確認し、今後の機会に活かそうかと思った次第です。
ただし、HPCの分野でよく見られる大規模疎行列の計算の高速化ではなく、大量の小規模行列の演算に着目します。FEMの用途によっては、大量の小規模行列の演算が必要になるからです。
#FEMによる構造計算の手順の概説
※これもテーマからかなり逸れるので、飛ばしても構いません。
例えば周波数応答の振動解析の場合、一般的な計算手順はおおよそ下記ようになります。
１． 解析対象の物体を$n$個の要素に分割する。
２． 着目したい周波数$ \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}
$$
３つの行列の掛け算を積分した形になります。各行列が何を表すかの説明は省略します。
今回、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）は、２つ行列[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 |
#結論
ともあれ、長く疑問に思っていた有限要素法における小行列の格納形式について、２つの手法を計測により比較した結果、A,B,Xを各個に連続してメモリに配置する(b)の方法が計算速度に有利、という結果が得られました。これを持って最終判断をするには至りませんが、参考値として活かせれば良いかと考えております。
参考文献
[1] Fortran90/95による有限要素法プログラミング, 藤井文夫,田中真人,佐藤維美著,丸善出版
[2] https://asteroidsathome.net/boinc/cpu_list.php