トップページへ戻る
過去の日記へ


予定


9月23日(金)

 学会でパワポを起動しようとしたら「認証しろ」というメッセージが。今回は マイクロソフトオープンライセンスと言うやつで、ダウンロードしてインストールしたもの。 その際「ログイン=認証してダウンロードしたから認証済みになってる」と思っていたのだが、 そうではなかったらしい。とりあえずボリュームライセンスサービスセンター(VLSC)にアクセスしようとしたら 「アカウントが凍結されています」・・・なんでだ? 凍結解除には携帯などの(アドレスではなく)電話番号の入力が必要で、 そこにショートメッセージで送られてくるアクセスコードが必要。 とにかくそれでログインして、該当プロダクトのライセンスキーを入力したら認証終了。

 その、この仕組「極めて」わかりづらい。特にVLSCのUIが最悪。 第一にこの画面を見て、ログインボタンが どこにあるかわかりづらい。他のログインが必要なサービスでは、「ログイン用のボタンだけ色を変える (twitter, Yahoo!メール)」 「ログイン/アカウント作成ボタンを右上に配置する (GitHub/Gmail)」等して「とりあえずそこを押せば良い」ことがわかりやすくなっている。対しVLSCでは、「サインインする」「認証をする」「自宅使用プログラムを確認する」「パートナーを見つける」という、論理的に同列でないボタンが別々の色で配置されている。しかも「認証」や「プログラムを確認」するためにはサインイン必須なので、まずは何はともあれサインインを促すべきである。ユーザに無用な選択肢を提示すべきではない。

 で、サインインを選んだ時の画面がこれ。 まずいところはいくつかあるんだけど、

 で、サインインした後、僕はMS Office 2016しか購入していないのに、製品が25個も登録されている。 自分が欲しいプロダクトキーがどれか探すのは素人にはかなり厳しいと思われる。原因は、Office 2016のライセンスに Office 2007、2010、2011、2013も含まれており、しかもSPが当たってるやつとか、Mac版とか、 CD版、DVD版など多数のプロダクトが「同列」に含まれてるから。

 これは、まず「最新のバージョンのみ提示」し、残りは「以前のバージョン」としてまとめるべき。 さらに「Mac版/Windows版」は、バージョンを指定してから選べるようにすべき。 もっと言えば、フィルタに「Windows版/Mac版」を付け加えるべき。もっと言えば「ダウンロード済み」という フィルタをつけるべき。

 もっと言えば、ダウンロード後認証が必要なら、ダウンロード画面でそれを促すべき。少なくとも僕はまったく気づかなかったぞ。

 その、MicrosoftのUIのまずさは、とにかく「論理的に同列にすべきでない情報を同列に提示してくる」ことに尽きる。 それにより「いま何をすべきか」が極めてわかりづらい。どうもMicrosoftには「ちゃんとしたUIの専門家」がいないのではないかと思われる。金はあるんだろうからちゃんとした人雇えばいいのに。

sent 4,632 bytes  received 31,873,582,975 bytes  3,348,241.78 bytes/sec
total size is 61,899,798,326  speedup is 1.94

 学会発表講演会。多数の講演があった。

 計算科学技術特論A、5頁ほど書いた。


9月22日(木)

 NOP


9月21日(水)

 いろいろあって眠くて仕方ないので、眠気覚ましのため、思想的に問題のあるSIMD化の続き。

 (1,1,1,0)、(2,2,2,0)、(3,3,3,0)、(4,4,4,0)という4つのベクトルが与えられた時、 (1,2,3,4)というベクトルを使いたい。三回演算は必要なはず。最初はunpackしてみる。 あとは適当にshuffleしてみて、所望の並びになった奴を採用。こんな感じ。

#include <immintrin.h>
#include <stdio.h>
//----------------------------------------------------------------------
void
printm256(__m256d r){
  double *a = (double*)(&r);
  printf("%f %f %f %f\n",a[0],a[1],a[2],a[3]);
}
//----------------------------------------------------------------------
int
main(void){
  __m256d r1 = _mm256_set_pd(0,1,1,1);
  __m256d r2 = _mm256_set_pd(0,2,2,2);
  __m256d r3 = _mm256_set_pd(0,3,3,3);
  __m256d r4 = _mm256_set_pd(0,4,4,4);
  printf("Before\n");
  printm256(r1);
  printm256(r2);
  printm256(r3);
  printm256(r4);
  printf("Unpack\n");
  __m256d r5 = _mm256_unpacklo_pd(r1,r3);
  printm256(r5);
  __m256d r6 = _mm256_unpacklo_pd(r2,r4);
  printm256(r6);
  printf("Shuffle\n");
  printm256(_mm256_shuffle_pd(r5,r6,12));
}
$ ./a.out
Before
1.000000 1.000000 1.000000 0.000000
2.000000 2.000000 2.000000 0.000000
3.000000 3.000000 3.000000 0.000000
4.000000 4.000000 4.000000 0.000000
Unpack
1.000000 3.000000 1.000000 3.000000
2.000000 4.000000 2.000000 4.000000
Shuffle
1.000000 2.000000 3.000000 4.000000

 できたっぽい。

 いや、やっぱり4倍展開して馬鹿SIMDせにゃならんかなぁ、と思って。 でも座標の読み込みのところをmovapdを使うことにこだわりつつ、 とりえあず距離の平方r2が(r2,r2,r2,0)の形で得られるので、 ループを4倍展開してこれを4つ作って、 (r2_a, r2_b, r2_c, r2_d)にしてから4つの力を同時に計算。 (f_a, f_b, f_c, f_d)にしてからまたバラす、という方向にしてみようかと。

   この方針って未来あるのか?えーと、力の計算って、r2ができたら、 r6を作るのに「r2*r2*r2」で積が2回。r14作るのに「r6*r6*r2」で積が2回。 (24 * r6 - 48)で積差が一回(時間刻み掛けるのは定数に含める)、除算が一回。 いま除算を何演算と数えるか知らんけど、とにかく命令としては6個しかない。 4倍展開したとして、パックに3命令かかる。で、4ペアの力計算ができたとして、 それを3つのレジスタにpermutationするので、やはり3命令か。 24命令を6命令に落とすかわりに、6命令増やすわけだから、実はわりといける?

 今調べたら、vdivpdのスループットが滅茶苦茶遅いということがわかったので、 馬鹿SIMD化なしのコードでvcvtpd2psで単精度に落としてvrcppsで逆数近似して、vcvtps2pdで倍精度に戻してから 一度精度補正するコードをためしてみた(1.356秒)。インテルコンパイラ(0.920秒)には勝てなかったが、GCC(1.748秒)には勝てた。 ちなみに精度補正付き逆数近似コードをフルアセンブラで書いてみたら、intrinsicより遅かった(1.588秒)。 さすがにこれだけ計算がややこしくなると、レジスタ割りつけがアホだとそれが見えてしまう。

 まー、フルアセンブラはネタとして、組み込み関数で逆数近似+精度補正+4倍アンロールによる馬鹿SIMD化、マスク処理による条件分岐を書く、というのが本命ですか。SoAにするかAoSにするかは正直微妙(この方針で性能に効くかわからない)・・・と。


9月20日(火)

 研究ミーティング。

 気泡論文の著者校正返した。

 学会発表講演会の概要を作って印刷をお願いした。本当は一週間前と前日にリマインダを出すらしいのだが、 今回ちょっといろいろ後手後手に。明日リマインドを出すだけにする。


9月19日(月)

 NOP


9月18日(日)

 NOP


9月17日(土)

 NOP


9月16日(金)

 学会最終日。自分の発表はそれなり。自分の研究である後半部分への質疑はわりとあったが、 気合を入れてしゃべった前半部分は聴衆の興味を惹かなかったようだ。conjectureだけ提示してもダメで、 やはり自分で結果を出せ、ということか。

 あと、前から思ってたんだけど、「研究の袋小路」にはまってる人がいる気がする。なんか「明らかに意味が無いこと」を 一生懸命やってしまっているみたいな。もちろん僕が間違っていて実はその研究には意味があるのかもしれないし、 僕も実は袋小路にはまりかけてるのかもしれないが。それにしたって「それは意味無いだろう・・・」と思う講演がいくつかあった。 やる前から「そうなるだろう」とほぼわかる計算をやって「そうなった」という奴とか、実験を説明するために理論にパラメータを入れるんだけど、そのパラメータはただ実験に合うように決めました、とか・・・。

 僕はあまり複数の人と共同研究したり、情報収集したり、というタイプじゃないんだけど、そういうタイプは 「袋小路」に陥りやすいんだよな。それは自覚しているのだが・・・。

 数名から「数独どうなった?」と聞かれたが、こちらは袋小路かどうかは知らないがストップ状態。 いま抱えてる9月末締め切りの原稿を脱稿するまではそちらまで手がまわらない・・・。


9月15日(木)

 水平SIMD化計画の続き。この計画の本質は、普通にSoAでMDをSIMD化しようとすると、アラインされていないロード/ストアにパック/アンパックが頻出し、 そのわりに計算が軽いので、SIMD化の加速よりもYMMレジスタへのload/storeのペナルティのほうが大きいのでは?という発想にある。 double q[N][4]みたいなAoS形式で持てば、3軸の座標のロード/ストアがアラインされてるので、vmovapdが発行できるから、 少しそこで無駄が生じても勝てるのでは?という思想。・・・というわけで、間接参照の練習として、こんなコードを書いてみる。

#include <immintrin.h>
#include <stdio.h>

const int N = 100;
const int D=4;
double q[N][D];
double p[N][D];

void
show256(__m256d x){
  double a[4];
  _mm256_store_pd(a,x);
  for(int i=0;i<D;i++){
    printf("%f ",a[i]);
  }
  printf("\n");
}

__m256d
func(int i, int j){
  __m256d vq = _mm256_load_pd((double*)(q+i));
  __m256d vp = _mm256_load_pd((double*)(p+j));
  return vq + vp;
}


int
main(void){
  for(int i=0;i<N;i++){
    for(int j=0;j<D;j++){
    p[i][j] = i;
    q[i][j] = i*2;
    }
  }
  show256(func(2,3));
}

 いや、こういう配列にしたとき、ちゃんとq[i][0]の先頭アドレスからロードさせるにはどうしたらいいか自信がなかったので。 で、実行するとちゃんとできている。ここで次元Dを4の倍数以外にしてみると、

$ icpc test.cpp
$ ./a.out
[1]    5029 segmentation fault  ./a.out

 うん、ちゃんとセグメった。funcの対応するアセンブリは

        movslq    %edi, %rdi                                    #20.19
        shlq      $5, %rdi                                      #21.41
        movslq    %esi, %rsi                                    #20.19
        shlq      $5, %rsi                                      #22.41
        vmovapd   q(%rdi), %ymm0                                #21.41
        vaddpd    p(%rsi), %ymm0, %ymm0                         #23.15
        ret

 あ、vaddpdって、片方メモリでもいいんだ。以前書いたやつは一度レジスタに落としてから書き戻してたから無駄だな。 あとshlqって何かとおもったらビットシフトか。 movslqってなんだ?と思ったら、signを意識した32bit→64bitmoveか。なんか似たようなのあったな、と思ったら、 ちょっと前にcltdって命令を調べたことがあった。 movslqとcltd/cltqって似たようなことやってない?と思ったら 同じことを思った人がStack Overflowで質問してた。 なんか歴史的経緯が書いてあるけどよくわからん。 あと、やはり僕はx86のアドレッシングについてよくわかってないな。rdiとrsiの使い方とか。

 で、ちょっとテストコード書いて、定数の扱いとかメモリアクセス少し工夫してみたけど、まだ インテルコンパイラに勝てない。ボトルネックが考えてたところと違うのかなぁ。

 いや、「水平方向SIMD」といっても、単に「movupdではなくmovapdを出す」 「パック/アンパックはしない」というだけのことなので、この状態でループ4倍展開して 馬鹿SIMD化をすれば、メモリまわりがスムーズなSIMD化ができるはず?途中で死ぬほどSIMDお手玉が必要になるけど。 で、それをやろうとしたんだけど、手でナイーブにやるにはレジスタがぜんっぜん足りなくて、 どうしたってスピルするんだけど、それを手で最適化するのは(少なくとも僕には)絶対無理。 っていうわけでintrinsicで書き直しですかね。

 っつーかインテルコンパイラは4倍展開してるのに、スピルを一度も起こさないようなコードを 吐いてるように見える。どうやってんだ?


9月14日(水)

 学会。バス混んでる。

 Cygwinのデフォルトシェルをzshにしたい。ググると、/etc/passwdを書き換えろとあるが、そんなものはない。 また、Cygwinターミナルであるminttyの実行オプションに /bin/zsh加えろというのもあるが、 そうするとパスがおかしくなったりして使えない。正解は/etc/nsswitch.confのdb_shellのところを以下のように書き換える。

db_shell: /bin/zsh

 デフォルトシェルを変更する前に、.zshrcのテストを十分やったほうがよかった。 今回適当にやったらログインできなくなって面倒なことになった。

 しかしCygwinのターミナルっていつのまにかminttyになってたのか。昔はcygwinってバッチファイルから起動してた気がする。 2011年に変わったとのことなので、もう5年前か。うーん。

 あと9日のキックオフミーティング中にTexLiveを入れた。標準インストールに2時間くらいかかった。いま試してみたら、何もしなくても日本語に対応しててすばらしい。

 gitの設定。

$ git config --global alias.st 'status -s' 
$ git config --global alias.ci 'commit -a' 

9月13日(火)

 そんな時間は全くないんですよ?無いんだけど、目の前の締め切り原稿からの逃避でSIMD。VPERMPDのテスト。

#include <immintrin.h>
#include <stdio.h>

void
printm256(__m256d r){
  double *a = (double*)(&r);
  printf("%f %f %f %f\n",a[0],a[1],a[2],a[3]);
}

int
main(void){
  __m256d r = _mm256_set_pd(3,2,1,0);
  printm256(r);
  printf("permute\n");
  for(unsigned char i=0;i<16;i++){
    printm256(_mm256_permute4x64_pd(r,i));
  }
}
$ ./a.out
0.000000 1.000000 2.000000 3.000000
permute
0.000000 0.000000 0.000000 0.000000
1.000000 0.000000 0.000000 0.000000
2.000000 0.000000 0.000000 0.000000
3.000000 0.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000
1.000000 1.000000 0.000000 0.000000
2.000000 1.000000 0.000000 0.000000
3.000000 1.000000 0.000000 0.000000
0.000000 2.000000 0.000000 0.000000
1.000000 2.000000 0.000000 0.000000
2.000000 2.000000 0.000000 0.000000
3.000000 2.000000 0.000000 0.000000
0.000000 3.000000 0.000000 0.000000
1.000000 3.000000 0.000000 0.000000
2.000000 3.000000 0.000000 0.000000
3.000000 3.000000 0.000000 0.000000

 うん、要するに4進法だと思えばいいんだな。ほしい並びは、「0123」から「1203」と「2013」なので、 それはそれぞれ「201」と「210」だ。

 次はインラインアセンブラで行く。ymmにローカル変数で触る方法がよくわからなかったのでグローバル変数で。

#include <stdio.h>
#include <emmintrin.h>
#include <immintrin.h>

typedef __m256d v4df;

double a[4] = {1,2,3,0};
double b[4] = {};

void
puts4(double a[4]){
  printf("%f %f %f %f\n",a[0],a[1],a[2],a[3]);
}

int
main(void){
  v4df *a4 = (v4df*)a;
  v4df *b4 = (v4df*)b;
  puts4(a);
  __asm__ (
      "vmovupd a(%rip), %ymm0\n\t"
      "vpermpd $201, %ymm0, %ymm1\n\t"
      "vpermpd $210, %ymm0, %ymm2\n\t"
      "vmovupd %ymm1, a(%rip)\n\t"
      "vmovupd %ymm2, b(%rip)\n\t"
  );
  puts4(a);
  puts4(b);
}
$ ./a.out
1.000000 2.000000 3.000000 0.000000
2.000000 3.000000 1.000000 0.000000
3.000000 1.000000 2.000000 0.000000

 よし、できてる。1230から、2310と、3120を作った。これでやりたいことをするための道具が揃った。

 やりたいことというのはLJの力計算の「水平SIMD化」。普通に「馬鹿SIMD」をするなら縦方向、つまりループをアンロールして ループがまわる方向にSIMD化すると思う。これは、規則的なメモリアクセスならいいんだけど、MDでは間接参照になってしまうため、 データのパック/アンパックが必要で、これがかなりのペナルティになる(気がする)。 さらに、if分岐のマスク処理も必要で、これも結構重い。さて、データ構造がAoSになっていれば x,y,zの情報にアラインされた状態にアクセスできるので、vmovapdで素直にymmに載せることができる。 この情報を使って横方向、つまりループ内SIMD化ができればパック/アンパックは不要で、3次元のベクトル演算が すなおに表現できる。というわけで組んでみた。

//------------------------------------------------------------------------
#include <stdio.h>
#include <emmintrin.h>
#include <immintrin.h>
//------------------------------------------------------------------------
typedef __m256d v4df;
const double dt = 0.01;
double a[4] = {1,2,3,0};
double b[4] = {1.5,2.6,3.7,0};
double c24[4] = {24*dt,24*dt,24*dt,0};
double c48[4] = {-48*dt,-48*dt,-48*dt,0};
//------------------------------------------------------------------------
void
puts4(double x[4]){
  printf("%f %f %f %f\n",x[0],x[1],x[2],x[3]);
}
void
check(double qi[4], double qj[4]){
  double dx = qj[0] - qi[0];
  double dy = qj[1] - qi[1];
  double dz = qj[2] - qi[2];
  double r2 = dx*dx + dy*dy +dz*dz;
  double r6 = r2*r2*r2;
  double df = (24.0 * r6 - 48.0)/(r6*r6*r2)*dt;
  printf("%f %f %f\n",df*dx,df*dy,df*dz);
}
//------------------------------------------------------------------------
int
main(void){
  v4df *c24v = (v4df*)c24;
  v4df *c48v = (v4df*)c48;
  check(a,b);
  __asm__ (
    "vmovapd a(%rip), %ymm0\n\t" // qi
    "vmovapd b(%rip), %ymm1\n\t" // qj
    "vmovapd c24(%rip), %ymm4\n\t" // ymm4 = (24,24,24,0)
    "vmovapd c48(%rip), %ymm5\n\t" // ymm4 = (48,48,48,0)
    "vsubpd %ymm0, %ymm1, %ymm0\n\t" // ymm0 = (dx, dy, dz,0 )
    "vmulpd %ymm0, %ymm0, %ymm1\n\t" // ymm1 = (dx**2, dy**2, dz**2,0)
    "vpermpd $201, %ymm1, %ymm2\n\t" // ymm2 = (dy**2, dz**2, dx**2,0)
    "vpermpd $210, %ymm1, %ymm3\n\t" // ymm3 = (dz**2, dx**2, dy**2,0)
    "vaddpd %ymm1, %ymm2, %ymm1\n\t" // ymm1 = (xy,yz,zx,0)
    "vaddpd %ymm1, %ymm3, %ymm1\n\t" // ymm1 = (r2,r2,r2,0)
    "vmulpd %ymm1, %ymm1, %ymm2\n\t" // ymm2 = (r4,r4,r4,0)
    "vmulpd %ymm1, %ymm2, %ymm2\n\t" // ymm2 = (r6,r6,r6,0)
    "vmulpd %ymm2, %ymm2, %ymm3\n\t" // ymm3 = (r12,r12,r12,0)
    "vmulpd %ymm1, %ymm3, %ymm3\n\t" // ymm3 = (r14,r14,r14,0)
    "vfmadd132pd %ymm4, %ymm5, %ymm2\n\t" // ymm2 = (24.0*r6-48)*dt
    "vdivpd %ymm3, %ymm2, %ymm1\n\t" //ymm1 = (24.0*r6-48)/r14*dt =df
    "vmulpd %ymm1, %ymm0, %ymm0\n\t" // ymm0 = (df*dx, df*dy, df*dz, 0)
    "vmovapd %ymm0, a(%rip)\n\t"
  );
  puts4(a);
}
//------------------------------------------------------------------------

 v4dfによる宣言は、インテルコンパイラが最適化でc24とかが使われていないと判断して消してしまうのを防ぐため。 で、テスト。

$ ./a.out
-0.041196 -0.049436 -0.057675
-0.041196 -0.049436 -0.057675 -nan

 うん、正しく計算できてるっぽいな。

 これは要するにフルアセンブラによるLJの力計算である。ポイントはvpermpdで(dx**2, dy**2, dz**2, 0)のベクトルを入れ替えてることで、 これを足すと(r**2, r**2, r**2,0)のベクトルになる。ここから工夫の余地はありそうだが、とりあえず 無駄に力を計算させてしまう。すると、力のベクトル(df,df,df,0)ができる。これに相対座標ベクトル (dx,dy,dz,0)をかければ、力積ベクトルが完成する。こいつはそのまま運動量にベクトル命令で足し込める。

 うん、一見して遅そうな気はする。でも、データの流れはスムーズだし、力の計算はまったくSIMD化できてない(どころか余計な組み換えが入ってる)けど、自明なベクトル演算(相対座標の計算、力積ベクトルの足し込み)は利用しているし、ループをまたいでないからif分岐によるジャンプも素直にかけるし、データのパック/アンパックが不要な分、レジスタ割り当てとか工夫したりしたら、わりといけないかな?無理かしら?

 ついでにループの練習でこんなの組んだ。

#include <stdio.h>

const int N = 10;
double a[N];
double b[N] ={};

void
func(void){
  asm volatile (
    "xor %eax, %eax\n\t"
    "mov $a, %ebx\n\t"
    "mov $b, %ecx\n\t"
    "LOOP:\n\t"
    "movq (%ebx), %rdx\n\t"
    "movq %rdx,(%ecx)\n\t"
    "add $8, %ebx\n\t"
    "add $8, %ecx\n\t"
    "inc %eax\n\t"
    "cmpl $10,%eax\n\t"
    "jl LOOP\n\t"
  );
}

int
main(void){
  for(int i=0;i<N;i++){
    a[i] = i;
  }
  func();
  for(int i=0;i<N;i++){
    printf("%f\n",b[i]);
  }

 これ、最適化レベルが低いうちは問題なく実行できるんだけど、最適化レベルを上げると「ラベルLOOPが多重定義されている」と怒られる。 しばらく原因がわからなかったんだけど、これ、funcがインライン展開されて、その際にラベルの名前を変えなくて、それで ラベルが二重定義になるんだな。なんかうまく回避する方法もあるらしいんだけど、簡単なのは関数funcにstaticをつけて inline展開を抑制してしまうこと。

 えーと、力の計算ができて、かつループの回し方がわかったので、これで本気でLJのフルアセンブラ(水平AVX化版)ができるようになったのか。 うーむ、やるのか?本当にやるのか?


9月12日(月)

 メーラにたまったメール数をみて、そっとじしたくなったが、そうもいかない。

sent 19,186 bytes  received 22,901,492,060 bytes  3,763,291.64 bytes/sec
total size is 43,664,793,931  speedup is 1.91

・・・的な(なにが?)

 スパコン使いとして生きていくつもりはあんまり無いが、やっぱりスパコンは好きだなぁ。 馬鹿でかい計算結果をローカルに落としてきて、解析スクリプトを走らせる瞬間っていつになってもちょっとアガる。 実験屋さんが合成した試料を取り出したり、はじめて何か測定する時もこういう気持ちなのかなぁ。

 v4dfのテスト。

#include <stdio.h>
#include <emmintrin.h>
#include <immintrin.h>

typedef __m256d v4df;

v4df func(v4df a, v4df b, v4df c){
  return a*b+c;
}

int
main(void){
  double a[4] = {1.0,2.0,3.0,4.0};
  double b[4] = {2.0,3.0,4.0,5.0};
  double c[4] = {3.0,4.0,5.0,6.0};
  v4df *a4  = (v4df*)a;
  v4df *b4  = (v4df*)b;
  v4df *c4  = (v4df*)c;
  (*a4) = (*a4) + (*b4) * (*c4);
  printf("%f %f %f %f\n",a[0],a[1],a[2],a[3]);
}
$ icpc -O3 -xHOST test.cpp
$ ./a.out

7.000000 14.000000 23.000000 34.000000

 あっさりできたぞ。配列を中途半端な大きさにしたり、変なの挟んだり、グローバル変数にしたりしたけど、 どうも勝手にアラインされているっぽい。__declspecとかやんなくていいのかしらん?

 力の計算をさせてみて、そのアセンブリを見る。vbroadcastsdは64bitを4つにコピーしてymmにつっこむ。 vmovupdはあらいんされていないデータをmoveする。

 ついでにSoA、AoSのチェック。コードはここの。

$ icpc -O3 -xHOST force.cc -o a.out 
$ icpc -O3 -xHOST force_soa.cc -o b.out     

$ ./a.out
N=20000, calcforce 0.869005 [sec]
N=20000, calcforce_hand 0.847749 [sec]
N=20000, calcforce_if 1.707234 [sec]
N=20000, calcforce_hand_if 1.535376 [sec]
N=20000, calcforce_hand_if_min 0.881039 [sec]
N=20000, calcforce_hand_if_condop 0.880944 [sec]

$ ./b.out
N=20000, calcforce 0.505722 [sec]
N=20000, calcforce_hand 0.458883 [sec]
N=20000, calcforce_if 1.696250 [sec]
N=20000, calcforce_hand_if 0.725826 [sec]
N=20000, calcforce_hand_if_min 0.466649 [sec]
N=20000, calcforce_hand_if_condop 0.466650 [sec]

 あ、もうダメだ。POWER以来ずっとAoSでだましだまし使ってきたけど、AoSはもうHaswellでは性能出ない。SoAへの書き換えが必要だ。


9月11日(日)

 NOP


9月10日(土)

 NOP


9月9日(金)

 ポスト「京」萌芽的課題キックオフミーティングに参加。 今後、またこういうのが増えるんだろうな。

 仙台は近い。「はやぶさ」で大宮から仙台まで一時間ちょっと。逆に寝るには短いくらい。 そういう意味では東京ー神戸とかだと、「あー、寝た」と思って目がさめてもまだ新大阪とかなのでちょうど良いよね、みたいな話を懇親会でしてた。あと、普段電車乗らないので、野田線でうっかり女性専用車に乗り込んで慌てて隣の車両に移動した。これ、前もやったんだよな・・・。

 帰りのはやぶさを予約するかちょっと迷ったのだが、懇親会がいつ終わるとかわからなかったので予約しなかった。 しかしはやぶさは混んでいて、しかも本数が(のぞみに比べて)あまり無いことを忘れてた。っていうか「のぞみ」の本数が異常なんだよな。 とりあえずダメ元でチケット確認したら「特急立席券」みたいなのがあって、それで「はやぶさ」にのれた。全席指定なのにそんなのがあるのか・・・。


9月8日(木)

 被扶養者の要件の確認、書類出した。非課税証明が必要だった。

 気泡論文がアクセプトされた。いや長かった、本当に長かったな。 「 Electronic Copyright License Agreement Form Due」というメールが来たので、 リンク先にとんで、I agreeにチェックを入れて送信。これで完了。

 論文のリポジトリ切ったのは去年の8月だ。そこから9月あたりまで執筆が続いている。 12月にリファレンスの追加をしているので、イントロ用に改めて論文のサーベイをしているっぽい。 そこからしばらくあいて、今年3月に集中的に書きあげ、論文を投稿。 予想以上に厳しい査読結果が5月に返ってきたので、追加計算を投げつつ、6月一杯は理論の理解が怪しかったところの補間に費やし、 古典的なLSW理論(というか、LSは簡単なんだけど、Wが面倒なのでそっち)を完全に理解し、7月に論文をほぼ丸々書きなおして、 長文の返事とともに再投稿。 Editorから「ちゃんと返事できてると思うけど英語がひどいから直せ」と言われたので英文校閲に出した。 7月末にレフェリーから二回目の査読レポート。最初のレポートで「この論文は全くダメだ」と書いてあったが「だいぶ書き換えてきたので、 まぁいいかと思い直した」とトーンが変わった。それに対してさらにていねいに返事をしたつもりだが、 8月24日に「返事に納得はいかないが、ここまでの努力はpublishに値するでしょう。ただこの文献だけもう一度読んでおけ」とほぼOKの返事をもらった。その文献とは、Hohenberg and Halperinの「Theory of dynamic critical phenomena」というレビューで、 もちろん読んではいたのだが、その分類との関連が不明瞭である、というのがもともとの査読者の指摘だったため、 それに対してコメントを追加して昨日再投稿。今日アクセプトが決まった。

 論文を書き始めてから一年ちょっとかかっており、 かつ論文に使ったデータを取得したのは平成25年度HPCI「京」一般利用枠の計算なので、2014年の1〜3月に計算したもの。 2年半前ですか・・・。その間に物性研のシステムが更新され、2.6PFのマシンが入った。 この論文のデータは「京」の4096ノードによる計算で、これは500TFとか。当時の物性研が180TFとかだからかなり大きな計算だったんだけど、 今は2.6PFあるので、同規模の追加計算が物性研で出来てしまった。

 これ、査読が長引いたのは、理論の理解が甘かったってのもあるんだけど、結局のところ最初の書き方がまずかったんだな。 ある古典的な理論があり、それは様々な発展を遂げるんだけど、今回の結果は古典的な理論で説明できましたよ、 という論旨なので、「様々な発展」は直接実験結果を説明しないから、引用しなくて良いかな、と省いてしまった。 それが、査読者に古い論文しか引用していないから「最近の様々な発展に全く無知である」と思われて、かなり印象が悪かったらしい。 「様々な発展はしたんだけど[引用]、今回は古典的な話を調べるよ」というイントロにすべきだった。反省。

 学会、今回は佐々さんのネーター発表の前座っぽい内容を話すつもりだったんだけど、発表がうまい具合に先で、 しかも僕の発表のセッションの座長も佐々さんだな。予定よりもっと前座よりに話を作り変えようかな。


9月7日(水)

 気泡論文、再投稿した。今度こそ間違いなく通るだろう。

 学会発表講演会の手配に不備があって慌てて対応する。

 計算科学技術特論A、書き始めた。30ページのうち、まだ2ページか。 参考文献とかつけなきゃいけないんだけど「闘うプログラマー」とか引用してやろうかな?


9月6日(火)

 学会のスライド作った。タイトル込み14頁か。ちょっと早口になるかなぁ。


9月6日(火)

 学会のスライド作った。タイトル込み14頁か。ちょっと早口になるかなぁ。


9月5日(月)

 査読レポートの返事書いた。三回目のレポートが来たのは8月24日か。ちょうど夏季休暇に入るところで、 その後発表が二件とかあったからなかなか手がまわらなかった。さすがにこれで通ると思われる。通るまで油断できないが。

 しかし今週末は仙台で萌芽のなんか、そして物理学会があって、さらに計算科学技術特論Aを30頁書かないといけない。 それで正直手一杯だよなぁ。ついでに萌芽は今月中に使いきらないといけないリソースがあるのだが、計算がちっとも流れない。 学会で発表する内容も論文にまとめないといけない。うーむ。

 そうだ、大規模ジョブがなぜ失敗したか考えないと。そもそも「京」フルノードで成功している計算なので、 32768ノードは確実に動くのだが、とある事情でハイブリッドではなくflat-MPIで投げたかった。 フルノード+ハイブリッドで走ったんだから、32768ノードflat-MPIは大丈夫だろ、と思って投げたら全然ダメだった。 もう一回ちゃんと計算しよう。

 今回の計算、途中で640の三乗のグリッドを用意する必要があった。つまり262144000個。 で、ワーク領域として、intを一つ、unsigned charを2つ使うので、1グリッドあたり6バイト。 つまり1.6GBくらい。これをルートプロセスに集める。

 ちなみに、フルノードの計算では、960×1080×640のグリッドを切った。これだけで4GBくらいか。これがノードあたり6GB程度しかくってないから大丈夫だと思ったのだが、今見ると、粒子の配列がベンチマーク用に大きくなってんな。これが原因か?

 えと、いまは訳あって粒子のデータを生配列で持っている。その大きさを1Mにしていた。 で、1粒子あたり48バイトだから48MB/プロセス。8倍しても384MB/ノードか。 えーと、さらに力の計算用にいくつかワークを持ってるけど、さらに倍にしても1GB行かないから大丈夫なはず。

 次、ペアリストサイズを最大粒子数の80倍にとってた。80M個で、この長さの整数の配列を3本持つので これは1GBオーダーか。これを8プロセス持ってしまうと8GBでメモリを圧迫する。ありゃ、こいつが原因だ。

 ベンチマーク用にとってた粒子数の最大数を1Mから0.3M個にした。そしたらペアリストのサイズが2.4GBになって問題なくなるはず。 だが、いまから32768ノードの計算は無理だろう。まずは4096ノードを流して、32768ノードの計算は後期に回す。

 なんか事前にコードが使う総メモリ量を簡単に計算する方法が欲しいなぁ。

 計算科学技術特論A、とにかくリポジトリは作ったぞ・・・。練習問題込みで今月中に30頁か。かなり厳しいな・・・。

 気泡論文がもっとスムーズに通ると思ったんだよなぁ。こんなに長くかかるとは誤算だった。おかげで得るところは大きかったけれども。

 光成さんのJITアセンブラ、Xbyakが Intelの解説サイトで紹介されている。

The JIT assembler Xbyak allowed us to use SIMD operations more efficiently.

 いや、すごいなマジで・・・。

 なんか、実際に作ったものがこうして使われるってのは嬉しいことだろうな。僕はもう長いこと「自分が作ったアプリケーションが他人に 使われる」という経験から遠ざかっている。僕はただプログラム片を書き散らかし、コンパイラについて文句を言い、 HPC用言語の将来について好き勝手御託を述べるが、自分で何かつくり上げることはできない。 そうして徐々に「自分がプログラマである」という意識が薄くなっていく・・・。


9月4日(日)

 NOP


9月3日(土)

 NOP


9月2日(金)

 最近スレッドまわりを調べているのだが、x86は命令数が多くて、知らない命令が多い。 今回もぽつんとcltdってのが出てきて「なんだこりゃ」と思ったら、32ビット整数を64ビット整数に 拡張するものらしい。具体的にはeaxの最上位ビットをedx全部にコピーする。 つまり、eaxが正ならedxは0、eaxが負ならedxの全てのビットが立つ。 なんだかよくわからんので、インラインアセンブラの練習を兼ねて確認する。こんなコードを書いた。

#include <stdio.h>
#include <stdlib.h>
void
bitdump(unsigned int a){
  for(int i=0;i<32;i++){
    printf("%d",a%2);
    a = a >> 1;
  }
}
int
main(int argc, char **argv){
  int a = atoi(argv[1]);
  int b = 0;
  __asm__(
  "xor %%edx, %%edx;"
  "mov %1, %%eax;"
  "cltd;"
  "mov %%edx, %0;"
  :"=r"(b)
  :"r"(a)
  );
  bitdump(b);
}

 整数aとbを用意し、a経由でeaxに値を設定してcltdを呼んだ後にedxをbに書き込んでビットダンプ。 僕の趣味でMSBが一番右に来るようにしている。

$ ./a.out 1
00000000000000000000000000000000

$ ./a.out 2
00000000000000000000000000000000

$ ./a.out -1
11111111111111111111111111111111

$ ./a.out -2
11111111111111111111111111111111

 なるほど、確かにそうなってる。問題はなんでraxとかあるのにこの命令が吐かれているかだが、よくわからん。

 あんまりチューニングとかゴリゴリするほうじゃないし、ましてアセンブリチューニングとかはしないんだけど、これまでPOWERとかSPARCとか、x86以外のマシンを使うと、なんかピーク性能比に比べて遅い場合が多くて、その原因を調べるためにx86以外のアセンブリは良く読むようになったんだよな。x86では基本的にインテルコンパイラがが僕より賢いし、アセンブリ読んでもよくわからんので、これまでx86のは読んでこなかった。でもまぁ、x86系人外(とは?)の影響もあって、少しずつ読んでいこうかと。

 スレッドローカル変数の説明を見つけて、だいぶわかった。 スレッドローカルな変数、Thread Local Storageの略でTLSと呼ばれるそうな。 スタックフレームに積まれる変数は自動的にスレッド固有になるから、ここで対象になる奴は グローバル変数かstaticな変数。グローバル変数やstatic変数に__threadキーワードをつけると、 特別なセクション.tdata(初期値がある場合)か,.tbss(ゼロクリアされる場合)に置かれる。 こんなコードを書いてみる。

int *
func(void){
  static __thread int a;
  return &a;
}

アセンブリを出すとこんな感じ。

func():
        pushq   %rbp
        movq    %rsp, %rbp
        movq    %fs:0, %rax
        addq    func()::a@TPOFF, %rax
        leave
        ret

fsはセグメントレジスタで、64bit環境におけるthread-self-pointerとして使われている。 コンパイルする時にg++ -m32をつけて32bitモードにすると、gsが使われる。

func():
        pushl   %ebp
        movl    %esp, %ebp
        movl    %gs:0, %eax
        leal    func()::a@NTPOFF(%eax), %eax
        popl    %ebp
        ret

 データをおくセクションは、初期化されていなければ.tbssに置かれる。

        .section        .tbss,"awT",@nobits
        .align 4        
        .type   func()::a, @object
        .size   func()::a, 4
func()::a:
        .zero   4

 「static __thread int a=1;」のように、変数に0以外の初期値を与えると、.tdataに置かれる。

        .section        .tdata,"awT",@progbits
        .align 4
        .type   func()::a, @object
        .size   func()::a, 4
func()::a:
        .long   1

 Mac上のgccでは、TLSをソフトウェアエミュレーションする。

func():
  pushq %rbp
  movq  %rsp, %rbp
  leaq  ___emutls_v._ZZ4funcvE1a(%rip), %rdi
  call  ___emutls_get_address
  popq  %rbp
  ret

 この___emutls_get_addressってのがアドレスを取得する関数だな。 ソースはここかな? clangはハードウェアでTLSを実現してるっぽい。


9月1日(木)

 ポスト京ミーティング。FDPSユーザがいたので驚いた。簡単に書けてわりとスケールするとのこと。 うーん、今考えている奴、フルで自分で組むのと、こういうのを利用するのとどちらが良いんだろうか。 個人的な好みとしては自分で組みたいんだけど、なにしろ時間が・・・


トップページへ戻る