コード並列化 - Githubissues

shimomurakazuya commented 3 years ago

openMPを実装したバージョンをアップしました。並列化箇所は初期条件の設定関数と、拡散方程式計算関数です。 d269a2fe960d6431abdc39e8996558d2cf17f2b2

デバッグの結果メモ：実行の時dataディレクトリがないとセグフォで落ちる。

hasegawa-yuta-jaea commented 3 years ago

あまり検証しないでどんどん先に進んでいっている雰囲気がありますが大丈夫ですか？後で致命的なバグが見つかるなどして、大幅な差し戻しにならないか心配です。

（十分に検証が出来ていて、Issueに載せていないだけということなら、いいのですが。。。）

shimomurakazuya commented 3 years ago

処置しました。切りがいいところなので、検証を優先します。

並列化の今後の目標メモ：計算時間計測昨日の実装

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea 83c027fc9ac0ad2f7dd796ccc9bcb919eefc05b9

OpenMPを用いた並列化を行ったコードです。

define.hでスレッド数を決めています。並列部分はValueDiffusion.cppの初期条件の設定と方程式の部分のみです。また、main関数にて分布関数の時間発展の計算にかかる時間を計測しています。

shimomurakazuya commented 3 years ago

このコードで、各スレッド巣ごとのシミュレーションにかかる時間を調べたのが下の表で、４スレッド数で最速になっているのを確認しています。

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

ありがとうございます。

速度が線形で伸びないのと4スレッドでサチるのが良くわからないですね。

このコードはどこで走らせていますか？（CPUのコア数は幾つですか？）

あと、秒では速いのか遅いのか分からないですね。ルーフラインモデルで見積もって上限の何十%の性能が出ている、とかのデータが必要かと思います。 https://qiita.com/tomo0/items/05f761ab580b4a9b6458

少しめんどくさいですがルーフラインモデルはHPCの基本なので、やってみてください。

shimomurakazuya commented 3 years ago

走らせてるのはおそらくログインノードですね...。（ジョブスクリプトなどで投げてなく、makeからの./runでやってます）

ルーフラインモデルに関しては了解しました。勉強します。　

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

走らせてるのはおそらくログインノードですね...。（ジョブスクリプトなどで投げてなく、makeからの./runでやってます）

なるほど。性能を取りたいときはジョブにしたほうがいいですね。

とはいえ、いちいちスクリプト書いてジョブ投入するのがめんどくさいのは良くわかります。

インンタラクティブジョブ（ジョブだけどssh ログインのようにコマンドが叩けるもの）もご検討ください。

alias _qrsh='qsub -I -q pc2 -l select=1:ncpus=40:mpiprocs=2:ompthreads=20 -l walltime=24:0:0 -P interactive -N interactive'

などと bashrc に書いておけば、それほど手間でもないと思います。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea 延安性能の見積もりに関して質問なのですが、演算性能の指標とする上限は各計算機の1ノードあたりの理論速度でいいのでしょうか？

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

使われているハードウェアがどれなのかを理解した上で見積もるべきですね。

それで、SGI8600は 1 ノード 2 CPUですが、1 CPU を使うのか 2 CPU を使うのかは実装によるので、そこは調べる必要があります。（1 CPU しか使わない場合は、単純にflops 半分になりますよね）

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea

シミュレーションコードの演算性能を求め肩について質問があります。今拡散方程式のループ文での演算性能を計算しようとしているのですが（演算性能）= (ループ文中のFLOPS)/(経過時間)　であるのは理解できるのですが、

ループ文中のFLOP数についてどのように求めるのかがわかりません。・ループ分には複数の数式が存在しているのですが、全体のFLOP数はそれぞれの式のFLOP数を足しあげれば良いのでしょうか？

・計算式の"*"1個でFLOP数が増えることは説明されたのですが、他の演算子"/"や"+"についても同様の処理を行えばいいのでしょうか、または無視できるほど早く処理されるのでしょうか？

・ループ文中で関数の呼び出しを行なっていますが、これがFLOP数にどのように影響するのかがわかりません。（単純に関数内の式を考えるだけで良いのでしょうか？）

以下に考えているループ文を載せます。（codeにあげたものと同じものです）

void ValuesDiffusion::
  time_integrate(const ValuesDiffusion& valuesDiffusion) {
      real* fn = f_;
      const real* f = valuesDiffusion.f_;
  #pragma omp parallel
      {
  #pragma omp for 
          for(int j=0; j<defines::ny; j++) ;
                  const int jm = (j-1 + defines::ny) % defines::ny;     
                  const int jp = (j+1 + defines::ny) % defines::ny;
              for(int i=0; i<defines::nx; i++) {
                  const int im = (i-1 + defines::nx) % defines::nx;
                  const int ip = (i+1 + defines::nx) % defines::nx;

                  const int ji = index::index_xy(i,j);
                  const int jim = index::index_xy(im,j);
                  const int jip = index::index_xy(ip,j);
                  const int jim2 = index::index_xy(i,jm);
                  const int jip2 = index::index_xy(i,jp);

                  fn[ji] = f[ji] +
                      + defines::c_dif * defines::dt / defines::dx / defines::dx * (f[jim] - 2*f[ji] + f[jip] )
                      + defines::c_dif * defines::dt / defines::dx / defines::dx * (f[jim2] - 2*f[ji] + f[jip2] );
              }
          }
      }
 }

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

浮動小数点数の四則演算 1 つにつき 1 flop です。（int の演算は数えません）

演算子ごとの速度を気にしてルーフラインを出すようなことは、普通はありませんが、一応補足すると、

+, -, * は同速度（1 cycle）
a + b * c は FMA命令といって、 1 cycle でまとめて処理されるので少し速い（が、これは、 + と *があるので 2 flop としてカウントします）
/ は遅い。（アーキテクチャにもよるが 10 cycle くらい）

shimomurakazuya commented 3 years ago

ジョブ投入をしようとすると、以下のエラー文が出てきますが、これはSGIのリソース申請でCPUノードの時間を申請しなかったからできないという認識で良いですか。また、プログラムを実行するにはpg1などのGPUノードの方に投げなければならないということですよね？

qsub: Error: CPU/ISV nodes' CPU working overtime. user=a214026

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

見たことのないエラーですが、

これはSGIのリソース申請でCPUノードの時間を申請しなかったからできないという認識で良いですか。また、プログラムを実行するにはpg1などのGPUノードの方に投げなければならないということですよね？

多分この認識で良いかと思います。（GPUノードでCPUコードを動かすのはあまり好ましくないですが、ジョブ投入できないなら仕方ないですね）

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

すみません1つ忘れていましたが、

SGI8600でGPUノードを使う場合は、CPUの性能も変わってくるので（CPUノードに比べて上がるので）ルーフラインモデルのピークflopsとメモリバンド幅は修正してください。

shimomurakazuya commented 3 years ago

一応スレッド数=1(並列化無しを想定)での演算性能を計算してみました。

演算時間=9.9123e-05 拡散方程式の時間発展のループ分におけるFlop:10 格子点数:nx=ny=100

演算性能=10∗100∗100/(9.9123𝑒−05)=1.00884×10^9 𝐹𝑙𝑜𝑝𝑠 ~ 1𝑇𝐹𝑙𝑜𝑝𝑠

一方で、理論演算性能について、条件は次の通り・ジョブ投入時のステータスは以下の通り select=1:ncpus=40:mpiprocs=1:ompthreads=12:ngpus=1

・GPUライブラリは呼び出していないので無視する。・1クロックあたりの演算数 = 822= 32 （これが一番わからない。FMAが関係しているのは分かったが）・1ノードあたりのクロック数は13.056かつselectの値がノード数に対応しているものとする。

こられの仮定の元で、理論演算性能は 33213.056 = 1.253376 TFlops

であるので、この時点で 1.00884/1.253376=80.5%の性能が出てることになります。

理論演算性能の計算の仕方がおかしいと思うののですが、どうでしょうか？

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

演算時間=9.9123e-05

これの単位はなんですか？ sec/step ですか？あと、短すぎると誤差が乗るので、ステップ数を 100 〜 1000 くらいにして平均を取った方がいいです。

格子点数:nx=ny=100

格子点数は多め（メモリ使用量が数 GBになるくらい）にした方がいいです。でないと良く分からないキャッシュが作用してよく分からない結果が出てきます。

・1クロックあたりの演算数 = 32 （これが一番わからない。FMAが関係しているのは分かったが）・1ノードあたりのクロック数は13.056かつselectの値がノード数に対応しているものとする。

すみません数字の意味が全くわかりません。

http://www.hpctech.co.jp/hardware/cascadelake-sp.html を見ると、Intel 6248RはDP Peak = 2304 GFLOPS とありますが。。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea 演算時間はsec/時間step　です。

格子点数のサイズについては了解しました。平均で算出してみます。

演算数についてはもう無視してもらっていいです。公式の表の値に基づきます。

＞1ノードあたりのクロック数は13.056かつselectの値がノード数に対応しているものとする。これも忘れてください。SGIのマニュアル読み直してみたら全然違いました。CPGPU演算部の総コア数が13056個あるのを読み違えました。）

shimomurakazuya commented 3 years ago

格子点数を1000にしてかつ各ループでかかる時の平均を計算し直してみました。

select=1:ncpus=40:mpiprocs=1:ompthreads=12:ngpus=1 コード内のスレッド並列数=1 Flop数:10、ループ回数:𝑛𝑥∗𝑛𝑦=1000、最小経過時間: 0.00996283

演算性能=10∗1000∗1000/0.00996283=1.00373×10^9 𝐹𝑙𝑜𝑝𝑠 ~ 1𝐺𝐹𝑙𝑜𝑝𝑠

昼ごろのコメントで算出した演算性能は1TFlopsと書いてありますが、単位換算を間違えていただけです。

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

遅くまでご苦労様です（勤務時間外なのであまり無理しないでくださいね。私は寝れないで酒飲んでるだけですので。。。）

1𝐺𝐹𝑙𝑜𝑝𝑠

承知です。この数字なら、経験的には妥当な気がします。

拡散方程式なら、ゴリゴリにチューニングしても、（OpenMPをフルで使って、かつアーキテクチャに合わせてループ構造とか最適化しても） peak flops比だと一桁%以下、ルーフラインモデルで >60% が出るくらいが妥当だと思います。

それ以上の最適化は、まあ頑張れば出来る場合もあるんですが、コスパが悪いので私はあまりやらないですね。。。（まあ60%→90%に上げるのに命を賭けてる研究者も居るといえば居ますが、あまりに大変なので、あまり真似しないですね。。。）

shimomurakazuya commented 3 years ago

先ほど井戸村さんと議論しましたが、並列かの精度検証の課題として以下の点が挙げられましたので共有します。

1.上の演算性能計算時のFlop数について10は誤りで8であること 2.nx=ny=1000での使用メモリ容量は8MBで全ての格子点上の分布関数データがオンキャッシュ状態（一時保存されてる状態？） 3.理論演算性能を出すべき。 4.割り算を使うのは好ましくない。後、もっとFlop数は減らせる。 5.今ループ文で周辺の点をindex関数を使って（メモリアクセスして）呼び出しているが、直接ベタ打ち（メモリアクセス無し）した時との比較もしてみたらどうか？ 6.そもそも数値は拡散方程式の性質に則っているのか？の検証。井戸村さんがいうにはji,jim,jip等がsheared変数だからおかしいとのことです。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea すみません、午前中のミーティングで確認したいことがあるのですが、下図のfor文でij,ijp等の変数が不定にならないのは、各スレッドで変数を定義して代入する作業を行なっているから結果的にprivate変数の扱いになっているという認識でよろしいですか？

仮に並列部分の外で変数を定義した場合はijは各スレッドで共有されて値が各スレッドで共通になり、正しい計算ができなくなる。調べた結果、こうなることは理解したのですが、並列ブロック内で定義した変数の挙動がわからないです。

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

下図のfor文でij,ijp等の変数が不定にならないのは、各スレッドで変数を定義して代入する作業を行なっているから結果的にprivate変数の扱いになっているという認識でよろしいですか？仮に並列部分の外で変数を定義した場合はijは各スレッドで共有されて値が各スレッドで共通になり、正しい計算ができなくなる。

これでよいです。

調べた結果、こうなることは理解したのですが、並列ブロック内で定義した変数の挙動がわからないです。

雑な理解ですが、

#pragma omp parallel for
for(int j=0; j<2; j++) {
   const int jm = j-1;
}

は

#pragma omp parallel sections
{
    #pragma omp section
    {
        int j=0;
        const int jm = j-1;
    }
    #pragma omp section
    {
        int j=1;
        const int jm = j-1;
    }
}

になると思ってください。（sectionの文法がこれでいいのかはわかりませんが。。）

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea 金曜日に指摘されたコードについて最適化を施しました。 a334509c77ecc497c46e8a613cafba8fc8c26e

後最適化後のループ文の理想演算性能についてなのですが、キャッシュの影響を考えない場合、これの演算強度はFlop数=6、データ転送数:3 、倍精度なので6/8*3 = 0.25 で、ピークバンド幅は140GBより

140*0.25 = 35 GFlopsで良いのでしょうか？

 #pragma omp parallel
 54     {
 55 #pragma omp for 
 56         for(int j=0; j<defines::ny; j++) {
 57                 const int jm = (j-1 + defines::ny) % defines::ny;
 58                 const int jp = (j+1 + defines::ny) % defines::ny;
 59             for(int i=0; i<defines::nx; i++) {
 60                 const int im = (i-1 + defines::nx) % defines::nx;
 61                 const int ip = (i+1 + defines::nx) % defines::nx;
 62 
 63                 const int ji  = index::index_xy(i,j);
 64                 const int jim = index::index_xy(im,j);
 65                 const int jip = index::index_xy(ip,j);
 66                 const int jim2 = index::index_xy(i,jm);
 67                 const int jip2 = index::index_xy(i,jp);
 68 
 69                 fn[ji] = f[ji]
 70                     + defines::coef_diff * (f[jim] - 4*f[ji] + f[jip] + f[jim2] + f[jip2] );
 71             }
 72         }
 73     }
 74 }

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

導出の過程がよくわかりません。。。（単位が書かれてないのでこっちで補完する必要があります）

演算強度はFlop数=6 → 6 flop/grid ？データ転送数:3 、倍精度 → 3*8 byte/grid ？ 6/8*3 = 0.25 → 0.25 flop/byte ？ピークバンド幅は140GB → 140 GB/s ？

ですか？

それと、

データ転送数:3

これは違うような。。。隣接点すべてキャッシュに乗るなら 2, 全く乗らないなら 6 ではないでしょうか。

あと、

140*0.25 = 35 GFlopsで良いのでしょうか？

今回はどう考えてもメモリ律速なのでこれでも良いですが、厳密なルーフラインモデルではピーク性能とのminを取る必要があります。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea わかりづらくて、すみません。単位としては１グリッドあたりで考えるとその単位形で間違い無いです。

後データ転送数に関してはzx=ny=1000でやっていて全ての点てオンキャッシュになります。 (キャッシュメモリ35Mbyteに対して、アクセス範囲が240Kbyte) なので演算強度は140*3/8 = 52.5flop/byteとなります。

最後のピーク性能についてはハードウェアのピーク性能で、今回は仕様書記載の2304Gflopsとのminをとる必要があるということでしょうか？

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

承知です。

後データ転送数に関してはzx=ny=1000でやっていて全ての点てオンキャッシュになります。 (キャッシュメモリ35MByteに対して、アクセス範囲が240KByte)

そうですか。この場合は、メモリバンド幅ではなくキャッシュのアクセス速度で評価する必要があります。（キャッシュの速度はよく知りませんが。。。）というか、キャッシュに乗らないサイズまで格子点数を増やした方がいいですね。

最後のピーク性能についてはハードウェアのピーク性能で、今回は仕様書記載の2304Gflopsとのminをとる必要があるということでしょうか？

その認識で合っています。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea キャッシュの速度は自分も分かっていません。仕様書見ても書いてなかったです。

というか、キャッシュに乗らないサイズまで格子点数を増やした方がいいですね。

これ、nx =10^7　ny=100とかにしても問題ないですかね？　これだとアクセス範囲が240MBになってキャッシュに乗らなくなるはずなのですが、格子間隔がxとyで10^5程度異なってしまいますので、計算結果の妥当性が保証できないと思います。追記；とりあえず計算結果については無視して、一回試してみます。

前にnx=ny=10^6でやろうと思ったのですが、流石にシミュレーション時間が10^4分オーダーになるので断念しました。

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

これ、nx =10^7　ny=100とかにしても問題ないですかね？　追記；とりあえず計算結果については無視して、一回試してみます。

自分の中で整合性が取れてればなんでもいいとは思うのですが、計算結果が検証できない（解析解がない）と分かっていてやるのは、あまりよくないですね。

10^nで調整する必然性はないので、2^n でいい感じのサイズになるところを探せばいいんじゃないでしょうか。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea データ転送数について質問があります。今nx =2000000,　ny=100で計算を行なっています。これでメモリのアクセス範囲はx方向:24byte, y方向に48Mbyteあり、y方向はキャッシュに乗らないようになっています。（キャッシュメモリは35MB）この場合のデータ転送数に関してですが、

x方向の差分を取るためのデータ(f[ji],f[jim],[jip])はキャッシュロードである一方で、y方向差分をとるためのデータ(f[jim2] , f[jip2])はメモリアクセうとなります。

したがってメモリアクセス数はロード数が3、ストア数が1なのでデータ転送数は４という解釈でよろしいですか？

                  const int ji  = index::index_xy(i,j);
                  const int jim = index::index_xy(im,j);
                  const int jip = index::index_xy(ip,j);
                  const int jim2 = index::index_xy(i,jm);
                  const int jip2 = index::index_xy(i,jp);

fn[ji] = f[ji]  + defines::coef_diff * (f[jim] - 4*f[ji] + f[jip] + f[jim2] + f[jip2] );

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

y方向はキャッシュに乗らないとありますが必ずしもそうではありません。実装によります。（OpenMPのスレッド並列の立て方、コンパイラの自動最適化によるベクトル化の具合によっては、y方向もキャッシュに乗る場合があります）

正確なキャッシュの利用状況を知りたければプロファイラで測定するしかありません。が、CPUでのプロファイラの使い方は知らないですね。。。

次善策ですが、

最良（1 load, 1 store）と最悪（5 load, 1 store）の2つルーフラインモデルで見積もる

というのはよくやる手段ですね。たいていの場合、実行性能は2本のルーフラインモデルの間に乗るようになります。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea 計算の結果、実行性能が最悪よりも低い値になりました。

Flop数:6、格子点数:𝑛𝑥=2000000,𝑛𝑦=100、平均演算時間: 0.108218(s) 演算性能=6∗2∗10^8/0.108218 =11.0887𝐺𝑓𝑙𝑜𝑝𝑠　に対して、

max（1 load, 1 store）=140∗ 6/(8 ∗2)=52.5𝐺𝑓𝑙𝑜𝑝𝑠 min⁡（5 load, 1 store）=140∗6/(8∗6)=17.5 𝐺𝑓𝑙𝑜𝑝𝑠

となりました。最適化ができてないみたいですが、格子点数の比率が悪いからでしょうか？

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

最適化の話はそう単純でもないので多角的に見ていくしかないですね。（沼にはまっていくとも言いますが。。。）

とりあえず思いつくのは、

OpenMPの並列数に比例してFLOPSが伸びるでしょうか？比例より伸びが悪い場合は、OpenMPのスレッドの立て方が悪い可能性があります
ベクトル化の件で思い出しましたが、Intelコンパイラ使ってますか？あとベクトル最適化のオプションは付いてますか？

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea intelコンパイラ使ってなかったです

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

そうですか。それなら遅いのはしょうがないですね。

intelコンパイラのオプションなどは、金曜に井戸村さんが紹介していた STREAM Benchmark の Makefile が参考になるかと思うので、それを見ながら拡散方程式の Makefile を修正してみてください。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea stream bench ではgccでコンパイルしていたので、それに習ってgccを用いようと思ったのですが、stdに関するリンクが読まれずエラーで終わります。stdのリンクが指定されているのに、どうしてなのでしょうか？調べてみてもg++に戻すのが解決策らしく、gccのままの解決策もstdのリンクを明示する(つまり一行目に書いてあるオプションをつけろ)というものでした。

 1 CXXFLAGS += -std=c++11 -O3 -Wall -Wextra
  2 LDFLAGS += -std=c++11 -lm -Wall -Wextra
  3 OMPFLAGS += -fopenmp
  4 
  5 .PHONY: all clean resultclean tagfiles
  6 
  7 TARGET := run
  8 
  9 SRCS := $(wildcard *.cpp)
 10 OBJS := $(SRCS:.cpp=.o)
 11 
 12 all: $(TARGET)
 13 
 14 $(TARGET): $(OBJS)
 15 #       g++ -O3 -o $@ $^ $(LDFLAGS) $(OMPFLAGS)
 16         gcc -O3 -o  $@ $^ $(LDFLAGS) $(OMPFLAGS)
 17 
 18 .cpp.o:
 19 #       g++ -O3 -c $< $(CXXFLAGS) $(OMPFLAGS)
 20         gcc -O3 -c $< $(CXXFLAGS) $(OMPFLAGS)

これでコンパイルすると、

gcc -O3 -o  run ValuesDiffusion.o Output.o Index.o main.o -std=c++11 -lm -Wall -Wextra -fopenmp
ValuesDiffusion.o: 関数 `_GLOBAL__sub_I__ZN15ValuesDiffusion15allocate_valuesEv' 内:
ValuesDiffusion.cpp:(.text.startup+0xa): `std::ios_base::Init::Init()' に対する定義されていない参照です
ValuesDiffusion.cpp:(.text.startup+0x19): `std::ios_base::Init::~Init()' に対する定義されていない参照です
ValuesDiffusion.o:(.eh_frame+0x13): `__gxx_personality_v0' に対する定義されていない参照です
Output.o: 関数 `_GLOBAL__sub_I__ZN6Output19OutputDiffusionDataERK15ValuesDiffusioni' 内:
Output.cpp:(.text.startup+0xa): `std::ios_base::Init::Init()' に対する定義されていない参照です
Output.cpp:(.text.startup+0x19): `std::ios_base::Init::~Init()' に対する定義されていない参照です
collect2: エラー: ld はステータス 1 で終了しました
make: *** [run] エラー 1

というエラーが出てきます。

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

gcc は g++ の C版です。C++コードで使うのは不正です。

stream benchですが、stream.iccというターゲットで intel版コンパイルのオプションがベタ書きされています。そっちを参照してください。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea

あまり重要じゃないですが報告します。今朝のミーティングでstreambenchのコンパイラはitnelでやったと言いましたが、よくみたらgccでやっていました。

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

そうですか。まあなぜか gcc を使うことになったので、ちょうど良かったんじゃないでしょうか。。

別件ですが、拡散方程式コードが遅くなりそうなところを1つ見つけました。 index_xy()がインライン展開されない実装になっているようなので、

Index.cpp を削除、

Index.hを

#ifndef INDEX_H_
#define INDEX_H_

#include "defines.h"

namespace index {
   inline int index_xy(int i, int j) {
       return  ( i + defines::nx* j )
   }
}
#endif

としてみてください。

shimomurakazuya commented 3 years ago

fccb85797d1d34e4140a2150ccf1cbfdde319d47

streambench用に作成したコードを上げました。このコードは一次元拡散方程式のコードを一つのファイルで書いたものを元にしています。main.cppはc(gccコンパイル)環境でmain.cpp_cpはC++環境(g++コンパイル)で動きます。

これらで実行した時の結果は以下の通りです。設定としては、 Nx=10000000, thread=24, Flop数:2,データ演算数:3、各ループごとにかかった時間の平均を出力しています。

streambenchのサンプルコードでの演算性能は10.605Gflopsです。この数字は時間の最小値で計算をしていたので、若干こちらの数字が小さくなりますが、c++環境のものは、1Gflops程度小さいです。またC環境の方は逆に理論性能値を若干超えてしまっています。が、これも複数回実行を行なったせいでキャッシュが残っていた影響が考えられるので、実際の性能は理論値を下回ると考えられます。

また、CとC++の結果を比較すると、C++の方が性能が低いです。ただ、この違いが出るのは、変数宣言並びにメモリ確保の分であり、ループ文とのk直接的な関わりはないです。それなのに数値の違いが現れるのは不思議に感じています。

スクリーンショット 2021-06-08 18 23 08

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

3点コメントです。

streambenchがなんのことかよく分かりませんが、コードから察するに STREAM TRIAD ですか？

STREAM benchmarkなら、先例に倣ってFLOPSより GB/s で書いた方がわかりやすい気がします。その方がオリジナルとも比較しやすいですし、バンド幅ならルーフラインモデルの計算ミスなどを機にする必要がないですよね。

あと、

キャッシュの影響？

これは、井戸村さんが言っていたように、nx をさらに10倍すれば確認できますね。

また、CとC++の結果を比較すると、C++の方が性能が低いです。ただ、この違いが出るのは、変数宣言並びにメモリ確保の分であり、ループ文とのk直接的な関わりはないです。それなのに数値の違いが現れるのは不思議に感じています。

Cのコードを読みましたが、

real f[defines::nx];

などと、static配列になっていますね？この場合コンパイル時に配列サイズが決まるため、動的メモリ確保よりも最適化が簡単です。速くなるのも納得です。（STREAMとは別物ですが、流体分野で有名な姫野ベンチでも、static と dynamic で数倍の差が出ることがあります）

それで、純粋にCとC++を比較したいなら、Cでもmalloc()関数を使うべきです。（#include <stdlib.h>で使えると思います）

GitHubの使い方でコメントですが、

ディレクトリの中がごちゃごちゃしてきてよくわかりません。
- CとC++、STREAMと1次元/2次元拡散で共存不可能なら（簡単に後戻りできないくらいコードが変わるなら）、1つのディレクトリで管理するのは厳しいかもしれません。ディレクトリを分けて
- dif1d_C
- dif1d
- dif2d
などとした方がいいかもしれません。あと、2次元拡散コードがGitHubから消えたようですが大丈夫ですか？
拡張子のお行儀が良くないです。C言語なら*.c, C++言語なら*.cpp とするべきです。また main.cpp_cp などと、一般的でない拡張子を作るのも控えましょう。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea

承知しました。ソースコードの方を従来の拡散方程式は「src」で、straembenchの方は「src_streambench」のディレクトリに分けました。（ついでにlog.txt等の不要ファイルを削除しました）また、火k調子についても修正を実施しています。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea

using とnamespaveについて少し質問なのですが、これらはC++でしか使えないのでしょうか?

拡張子を.cppから.cに変えるとvim画面でこれらの文字に色がつかなくなって、以下のコンパイルエラーがでます。

main.cpp:13:7: エラー: expected nested-name-specifier before ‘real’　
using real  = double;

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

using とnamespave

いずれも C++ 専用です。

Cでは、

usingは、typedef double real; に置き換えてください
namespace の代わりになる機能はありません。namespaceなしで（グローバルにベタ書きで）対応してくだい。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea C++の環境(g++ -o Wall)にて静的メモリ確保(malloc不使用)と動的メモリ確保(malloc使用)で演算時間を調べていたのですが。メモリ確保の方法の違いで演算時間に影響はなく、設定や投入の投げ方である程度時間がばらつくことがわかりました。

まとめとしては、 1、最短時間で比較すると、動的メモリも静的メモリもほぼ等しい。 2、格子点数が10^7の時、複数回投げると後半ほど演算が早くなる傾向にある。３、格子点数が10^8の時そ、の傾向はほぼなくなる。しかし演算時間にばらつきが見られる。

　下の表は格子点数を10^7、10^8で静的または動的なメモリ確保した時のループの演算時間を計測した表です。10^7の場合はそれぞれ5回ジョブを投げました。その結果としては回を重ねるごとに演算時間が減少していく傾向にあることがわかりました。これについてはキャッシュが影響していると考えています。

　nx=10^8の場合に関しては静的メモリを5回、動的メモリの方を7回投入しました。基本的に個々のジョブが完了してから次のジョブを投入していたのですが、動的メモリの方の2〜5回目については完了を待たずに投入しています、その結果3~5回目にかけて演算時間が増加しています。　一方で、静的メモリの方は全て個々のジョブが完了してから次のジョブを投入していたのですが、2,5回目にかけて演算時間が増大しています。これについての原因に心当たりはなく、謎のままです。

　複数のジョブを同時に投入すると演算時間が遅くなる理由や、回を重ねるごとに演算時間が減少する理由についてはよくわかっていませんが、いずれにせよ演算性能の最高効率を見るためには、複数回の計測が必要であると考えています。

　また、昨日Cの環境(静的メモリ確保)の方が早いと言ったことについても、この傾向が関与していたと考えられます。表の数値においても最短時間の数値に差はほぼないため、静的メモリ確保と動的メモリ確保の違いはループの演算時間に影響を与えないと私は考えています。

追記(16時)表をバンド幅の数値を追加したものに差し替えました。10^7、10^8いずれの場合もハードウェアピーク値を超えているのでキャッシュの影響が少なからずあると考えられます。

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

メモリ確保の方法の違いで演算時間に影響はなく

そうですか。承知です。 C言語での結果が無いようですが、Cで書くのは辞めたのですか？（まあCityLBMはC++なので、Cが無いならないで良いです）

10^7の場合はそれぞれ5回ジョブを投げました。その結果としては回を重ねるごとに演算時間が減少していく傾向にあることがわかりました。これについてはキャッシュが影響していると考えています。

これはよくわからないです。ジョブと言っているのは、 qsub を5回叩いたということですか？だとしたら、計算ノードが変わるはずなのでキャッシュの影響というのは考えづらいです。いずれにせよ、10^8 の結果があるならそっちだけ評価すればいいと思いますが。。

複数のジョブを同時に投入すると演算時間が遅くなる理由や、回を重ねるごとに演算時間が減少する理由についてはよくわかっていませんが、いずれにせよ演算性能の最高効率を見るためには、複数回の計測が必要であると考えています。

複数ジョブもそうですが、1ジョブ内で複数回計測もやった方がいいです。例えばオリジナルのSTREAMベンチでは、NTIMESのマクロの分だけ for を回していますね。

あと、遅い方の結果は気にしてもしょうがないところがあります。たまたまネットワークが逼迫していたとか、一時的に電力不足になっていたとか、壊れかけの計算ノードを引いてしまったとか、我々の与り知らぬ理由で遅くなることもよくあるので。。。ですので、最適値だけ参照するようにしてください。（とはいえ、突然10倍以上遅くなるとか、あまりに酷い場合はバグを疑う必要がありますが。。。まあ2倍遅いくらいなら個人的には許容範囲です）

あと、余談ですが、

表でバンド幅が空白になっていますが、毎回手で計算するのが面倒くさいのであれば、コード中でバンド幅を計算して printf すれば良いかと思います。
表のスクリーンショットが大きくてみづらいので、markdownの表の機能で書いてもらえると助かります。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea Cへの書き換えは諦めました。あまり最適化において本質的ではなさそうなので。

ジョブの投げ方（qsubの叩き方）は qsub →ジョブ実行、完了→結果調査、数値表に貼り付け→qsub　という流れです。

また、演算時間の出力値は各時間ステップごとの演算時間を計測し、平均値を出力しているので、同じことをしているはずです（オリジナルの方は最小値でバンド幅を計算していますが）。

表のお大きさに関しては修正しておきます。

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

ジョブの投げ方（qsubの叩き方）は qsub →ジョブ実行、完了→結果調査、数値表に貼り付け→qsub　という流れです。

これは了解です。問題ないと思います。

演算時間の出力値は各時間ステップごとの演算時間を計測し、平均値を出力しているので、同じことをしているはずです

同じではないですね。

言いたかったのは、例えば、 1つのジョブ中で 1000 step の計算を 10回（計1万 step）繰り返せば、後半のどこかの 1000 step で最適な値が出るんじゃないでしょうか、ということです。

あと、最初の数秒間の計算は、初期化直後ということで速度が不安定になりやすいので計測に含めない方がいいです。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea

言いたかったのは、例えば、 1つのジョブ中で 1000 step の計算を 10回（計1万 step）繰り返せば、後半のどこかの 1000 step で最適な値が出るんじゃないでしょうか、ということです。

しょうちしました。初期化直後の速度が不安定なのは知らなかったです。コードを改良してみます。

shimomurakazuya commented 3 years ago

@hasegawa-yuta-jaea

すみません。ちょっと確認んしたいのですが下の画像の拡散方程式擬きを計算する関数のデータ演算数はストア１、ロード2の3で合ってますか？　(flpos数は4)

スクリーンショット 2021-06-10 15 58 57

hasegawa-yuta-jaea commented 3 years ago

@shimomurakazuya

ロード2回の根拠がよくわかりません。。。

fn[i] = ... の代入処理が load & store になるという意味ですか？

古いCPUだとそういうアーキテクチャもあるそうですが。。 SGI8600のCPUがそうなのかはよく分かりません（ググってそういう情報が出てくるなら、それでいいです）

あとすみません、完全に別件ですが、しばらくメールとかGitHubとかを確認する頻度が落ちそうです。（たぶん、12時前、17時頃、寝る前の3回くらい）返信が遅れ気味になりますがご了承ください。

急ぎの案件があったら、LINEか電話で呼び出してください。その時はZOOMか何かで対応します。

shimomurakazuya / diff_eq

コード並列化 #6