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


予定


6月24日(金)

 公式発表でた→ポスト「京」萌芽的課題アプリケーション開発の実施機関決定について

 僕は「萌芽的課題1 基礎科学のフロンティア - 極限への挑戦」に所属していることになっていて、大規模MDをやります。 大規模短距離古典MDという手法そのものの追究と、大規模MDを使って初めて見える物理の追究をする予定です(つまり、通常運転ということ)。

 ただし、せっかくこういうプロジェクトに入って、計算資源とかもいただくので(その分研究会とか報告会とか報告書とか増えるけど)、 このプロジェクト中にもういちどGordon Bellに挑戦しようと思っています。なにぶんそっち方面は素人なので、ご指導ご鞭撻のほど、よろしくおねがいします>そっち方面でプロな皆様


6月23日(木)

 並列MDコードを拡張したい。いまは粒子種も相互作用も一種類しか使えないのだが、これを多数の粒子が混ざるようにしたい。 また、二体相互作用だけでなく、三体相互作用もいれたい。これらを気軽に追加できるフレームワークが欲しいが、速度はそれなりに出て欲しい。 開発コストはなるべく抑えたい。さて、どうするか。

 ややこしいことはしないことにしよう。適当に要求要件を書く。以下はたわごと。これから数年でなんかする。

 面倒なのは、ボンドの情報の受け渡しだなぁ。 どの粒子とどの粒子がボンドで繋がれているかの情報をどうするか。一番簡単なのは、粒子番号をグローバルにしてしまって、 ボンド情報を全員で共有してしまうことだけれど、これは現実的ではない。うーむ。

 とりあえず似たようなコンセプトであるはずのFDPSをダウンロードして、使えそうな設計方針をパクることにしよう。触る時間は全く無いが、とりあえずDLだけしておく。

 なんかあっさり動いた。

 vdw-testが(おそらく一次元の)LJ粒子のサンプルっぽいが、見てもまだよくわからないな。

 どうでもいいけど、スレッドをプロセスのようにして、あたかもflat-MPIのようにハイブリッド並列する「いんちきflat-MPI」は、 コアあたりそれなりの粒子数が無いと性能が出ない。ということは、それなりの処理能力のコアと、それなりのコアあたりのメモリが無いと成立しない。 これは昨今の進化の方向と真逆なので、早晩使えなくなるかもしれない。そしたらガチハイブリッド書くのかなぁ。それはイヤなんだけどなぁ。

 なんだかその、いろいろなんだかなぁ。


6月22日(水)

 「るろうに剣心」という漫画に魚沼宇水という登場人物がいる。この人は志々雄に敗れて失明したため、「いつでも命を狙って良い」という 約束で志々雄に協力していることになっている・・・が、後に別の登場人物から、「すでに(圧倒的な実力差に)復讐はあきらめており、命を狙う振りをしていれば周囲には敗北をさとられずに済む」と看破される。

 こういうのってわりと多いというか、自分もそういうところがある。何かに挑戦したい。だが、挑戦してしまうと結果が出てしまう。 自分の力が足りないことが露呈してしまう。だから、挑戦できないふりをする。挑戦できない理由を探し、「挑戦したいができない」という構図を作りたがる。 こういうのを全て否定できるほど、人の心は強くないと思う。でもやっぱりそういうのはカッコ悪いので、挑戦しないと。

 というわけで、今後も「二つの世界一」を目指す所存。

 もう全く時間取れないんだけどさ(←言ってるそばから逃げ口上)。


6月21日(火)

 メール書いてばかりで精神衛生上良くない。コンパイラの吐いたアセンブリでも眺めて最適化に文句でもつけよう(←たち悪い)。

 以前、構造体を含んだ関数呼び出しの最適化について調べたが、最近Fortranにも構造体があることを知ったので、それを試してみる。 こんなコード。

program type_sample
  implicit none
  type point
    integer::x, y
  end type point
  type(point) :: p1,p2,p3
  p1%x = 1
  p1%y = 2
  p2%x = 3
  p2%y = 5
  p3 = add(p1,p2)
  print *,p3%x,p3%y
contains
  type(point) function add(p1,p2)
    type(point) ::p1,p2
    add%x = p1%x + p2%x
    add%y = p1%y + p2%y
    end function add
end program type_sample

 整数二つからなる構造体を作り、それを二つ受け取ってそれぞれの和をメンバに持つ構造体を返す関数addを呼ぶ。 ちゃんと解析できれば、function addはまるまる消えて、かつ足し算もせず、即値で返せるはず。 まずIntel Fortranの-Ofastの結果。

..B1.2:                         # Preds ..B1.7
        movl      $4, %eax                                      #11.3
        lea       (%rsp), %rdi                                  #12.3
        movl      %eax, type_sample_$P3.0.1(%rip)               #11.3
        movl      $-1, %esi                                     #12.3
        movl      %eax, 64(%rdi)                                #12.3
        movq      $0x1208384ff00, %rdx                          #12.3
        movl      $__STRLITPACK_0.0.1, %ecx                     #12.3
        lea       64(%rsp), %r8                                 #12.3
        xorl      %eax, %eax                                    #12.3
        movl      $1, type_sample_$P1.0.1(%rip)                 #7.3
        movl      $2, 4+type_sample_$P1.0.1(%rip)               #8.3
        movl      $3, type_sample_$P2.0.1(%rip)                 #9.3
        movl      $5, 4+type_sample_$P2.0.1(%rip)               #10.3
        movl      $7, 4+type_sample_$P3.0.1(%rip)               #11.3
        movq      $0, (%rdi)                                    #12.3
        call      for_write_seq_lis                             #12.3
                                # LOE rbx r12 r13 r14 r15
..B1.3:                         # Preds ..B1.2
        movl      $__STRLITPACK_1.0.1, %esi                     #12.3
        lea       (%rsp), %rdi                                  #12.3
        movl      $7, 72(%rdi)                                  #12.3
        lea       72(%rsp), %rdx                                #12.3
        call      for_write_seq_lis_xmit                        #12.3

 なんだ?なんかえらく非効率なコード吐いてるぞ。Fortranのprint文の処理ってよくわからないけれど、 おそらく最初の出力がfor_write_seq_lisで、次から毎回for_write_seq_lis_xmitが呼ばれるっぽい。 表示されるデータは、最初は64(%rdi)、次が72(%rsp)・・・と8バイトおき。ただ、文字列がまざるとなんか処理がややこしくなる。 まず、print文に渡される文字列は全て最初にまとめられる。で、64(%rdi)には、「文字列のサイズ」が入る。この時、オフセットは16バイト飛ばされる。 $__STRLITPACK_0.0.1とかはフォーマット指定かな。最初が%ecx、次が%esiに渡されてるのがよくわからんけど。

 そんなことはともかく、見た瞬間「print *,4,7」にできるコードなのに、構造体を全て作っている。 同等なC++コードはこんな感じ。

#include <stdio.h>
struct point{
  int x,y;
};
point add(point &p1, point &p2){
  point p3;
  p3.x = p1.x + p2.x;
  p3.y = p1.y + p2.y;
  return p3;
}
int
main(void){
  point p1,p2;
  p1.x = 1;
  p1.y = 2;
  p2.x = 3;
  p2.y = 5;
  point p3 = add(p1,p2);
  printf("%d %d\n",p3.x,p3.y);
}

 icpcは当然即値を吐く。

        movl      $.L_2__STRING.0, %edi                         #19.3
        movl      $4, %esi                                      #19.3
        orl       $32832, (%rsp)                                #12.11
        movl      $7, %edx                                      #19.3
        xorl      %eax, %eax                                    #19.3
        ldmxcsr   (%rsp)                                        #12.11
..___tag_value_main.6:
#       printf(const char *, ...)
        call      printf                                        #19.3

 g++ももちろん即値。

        movl    $7, %edx
        movl    $4, %esi
        movl    $.LC0, %edi
        xorl    %eax, %eax
        call    printf

 なんでFortranは即値吐けないの?もしかして変数がグローバル扱いになってて、どこか別の場所から参照されるかもしれないから?

 富士通コンパイラにいたっては関数addのインライン展開すらしない。

  sxar1 
  mov 1,%xg4  
  or  %l0,%m44(.LS1+4096),%l0
  sxar1 
  mov 2,%xg5  
  sllx  %l0,12,%l0  
  sxar1 
  mov 3,%xg6  
  or  %l0,%l44(.LS1+4096),%l0
  sxar1
  mov 5,%xg7  
  add %l0,-3960,%o0
  sxar1 
  stw %xg4,[%l0+-3960]
  mov %fp,%g2 
  sxar2 
  stw %xg5,[%l0+-3956]
  stw %xg6,[%l0+-3952]
  sxar1
  stw %xg7,[%l0+-3948]
  call  type_sample.add_

 全部メモリに構造体を作って、そのポインタをtype_sample.add_に渡している。うーん。 数値計算とかで、粒子の座標と運動量、粒子種その他を構造体にまとめて、それを配列にしたいことなんてよくあると思うんだけど、 こういうことされたら激遅になるよなぁ。僕がFortranわかってないだけで、ちゃんと回避策があるんだろうけど・・・。

 なんかアセンブリが読みづらくて仕方ないので、Syntax Highlighterのx86asm版作成。まだレジスタと命令の一部だけ。 $から始まる即値をちゃんと表示する正規表現がまだできてない。


6月20日(月)

 メールしたり相談したり電話したりしてたら一日が終わった。

 先月末から数独プロジェクトは全く手をつけることができず、自身の研究時間を確保することも難しい状況。多忙を言い訳にするのは無能の証というけれども・・・。 とにかくいま抱えているこの論文を再投稿しないことには次に進めない。


6月19日(日)

 泣きながらドイツ語の論文読んでた。多分理解したはず。査読レポートへの返事は4週間で返せと言われているが、完全にぶっちぎってる。


6月18日(土)

 NOP


6月17日(金)

 マナーや習慣が話題となると必ずあらわれる、いわゆる「マナー警察」が好きじゃない。 もちろんそのマナーや習慣にはある程度の意味はあるのだろうが、それが確定したのがわりと最近であることが多い、ということは知っておいて損はない。以前書き順でも話題にしたように、何かの混乱を避けるため、便利のために制定された「基準」が、いつのまにか「ルール」として一人歩きし、それ以外のものを排除しようとする動きになるのは興味深い。

 で、句読点の話。句読点が「、。」ではなく「,。」が正式だとか、横書きと縦書きで異なるとかいう話はよく聞く。 その根拠の一つが昭和21年3月に文部省(当時)から発行された「くぎり符号の使ひ方〔句読法〕(案)」であった。 その内容は文化庁のサイトから ダウンロードできる。正式には「国語表記の問題」という文部省から昭和38年に発行された本の付録として、 昭和21年に発行された「くぎり符号の使ひ方〔句読法〕(案)」が掲載されているを閲覧できる。

 その内容を参照すると、まず平成27年に策定された「公用文作成の要領」に、

句読点は,横書きでは「,」および「。」を用いる。

と指定されている。また、「国語の書き表わし方」という資料では、くぎり符号としては

  1. 。 まる
  2. 、 てん
  3. ・ なかてん
  4. () かっこ
  5. 「」『』 かぎ

の5種のみ定めており「、」が正式という印象だが、「横書きの場合の書き方」が別途設けられており、

くぎり符号の使い方は,縦書きの場合と同じである。ただし、横書きの場合は「、」を用いず「,」を用いる。

とあるので、これらを見る限り、読点は縦書きではてん(、)、横書きではカンマ(,)が正式のように見える。

 で、この資料を見ていると面白い。「マルの打ち方」にいろいろ例が出ている。 まず「」(カギ)の中でも文の終止にはうつ、とあり、例として 「どちらへ。」「上野まで。」とある。僕は「」で区切られた内部の最後はマルをつけていなかった。 手元の本を何冊か見たが、カギの中の文の終止にはマルは打っていない。

 また、「文の終止でカッコをへだててうつことがある」と定められている。こんな感じ。

このことは、すでに第三章で説明した(五七頁参照)。

 これはいいのかな?と思ってたが、これが正式らしい。ちなみに英語の引用符では、ピリオドを内側に書くと習った気がするが、 いま調べると外側っぽいなぁ。

 で、ここからが本論だけど、所謂三点リーダ「…」の話。先の「くぎり符号の使ひ方〔句読法〕(案)」によると、 まず呼び方が「テンテン」「テンセン」となっている。符号として「……」「…」「…………」と、2こ、1こ、4こが挙げられている。 で、「三点リーダは偶数個使う」ということが明示されているわけではないが、例文に出てくるのは2個バージョンと4個バージョンのみ。

 ネットでちょっとさらった限りでは、「三点リーダは偶数個で使用する」というのは出版の際のルールであって、特に国が定めた基準や標準に依拠するものではなさそう。っていうかなぜ偶数個にするのか納得のいく理由を見つけられなかった。いや、活字を組む都合とか、そういう物理的な要請由来なら納得できるんだけど、いまのところ「読みやすさ」とか「美しさ」とかそういう主観的なものしか見つけられていない。

 しかし、こういう見栄え系にはよく「警察」が湧きやすい。僕が見たことがある「警察」は、「スーツの一番下のボタンは止めない警察」とか。スーツにおいて一番下のボタンはそもそも留めるようにできておらず、留めるとスタイルが悪くなってしまうとか言うやつ。

 どちらも僕には全く理解不能で、ボタンが飾りで留めてほしくないならば対応する穴も飾りにして、実際に留められないようにしておけば良いだけだし、三点リーダも、どうしても偶数個しか使うべきでないなら、偶数個でしか入力できないようにすべき。

 「理解不能」と言いつつ、やっぱり理解できる気もする。人はこういう「踏み絵」が大好きなんだ。 知らないとうっかりやってしまいがちな「トラップ」を多数用意して、それで他人との差を付けたがる。 そのトラップは日常生活に支障をきたしてはいけない(例えば操作を誤ると爆発するような機械にはそういうトラップはしかけることができない)。 そんなわけで、「トラップ」は、失敗しても笑われるだけで済む礼節・マナーの分野に頻出するんだよな。 一度自分がひっかかって笑われた後、多くの人は「次の犠牲者を笑う側」に移行するんだろう。

 ・・・嫌な感じである(・・・とあえて三点リーダの使い方を間違える私)。

 ちなみに、いまのところ三点リーダが偶数であらねばならない理由で、僕が一番納得したのが、手書き原稿だとカタカナのミと紛らわしいから。これはわかりやすい。

 まぁ、もし僕が小説書いてどこかに投稿するなら、そういうルールは守ると思う。論文投稿する際に、リファレンスの形式をどうするかとか、そんなレベルの話。ただ、そういう「警察」、もっというと「自警団」が湧くのを見るのはイヤだな、という話。

 チェビシェフ多項式による多項式近似。x**nとT_n - x**nを表示している。ただしT_nの最高次の係数は1とする。

 次数があがるにつれて急速に近似がよくなる。そういう性質をもった関数がこんな簡単な漸化式で作れるとは不思議。


6月16日(木)

 一般公開向けのネタとして、棒の両端におもりをつけたものをモーターで回して「重力波発生装置」と題し、 その波の大きさのオーダーの見積もりと共に展示する、なんてことを考えたのだが、「まぁ、同じネタ考えている人がいるだろうな」と思ってちょっと先行研究を調べてみた。そしたら重力波関連の特許が出るわ出るわ・・・。飛行、発電、そして通信・・・。 その、ちょっとだけ「出願人」の名前でググッて、すぐに後悔した。これ系の話は深入りすると気持ち悪くなるからこれくらいにしておこう。

 昔はこういう話は好きだった(トルマリンとかタキオン健康法とかね)。最近だと水素水とか? だが、最近になって、こういう話は軽々しく話題にすべきではないと考えるようになった。 なんでだろう?ということを考えたのだが、結局のところ僕が「人の間の直接相互作用」だからなんだと思う。

 本でも仕組みでもソフトウェアでも芸術でもなんでもいいんだけど、人間の作ったものは好き。でも人間そのものは好きじゃない。 もっというと、人の間の直接の相互作用は好きじゃない。こういうトンデモ系の話には、必ず「それがトンデモである」と指摘する「トンデモ警察」があらわれる。だが、トンデモ系の人はもとより議論などするつもりは無いので(結論ありきで「議論」を行う)、全く「議論」は噛みあわない。結局お互いに感情的な何かを吐露して終わる。トンデモ系をちょっと突っ込んで調べると、すぐにこういう「人の間の摩擦」が見えてきて気持ち悪くなる。

 いや、トンデモをトンデモと指摘するのは良いと思うし、その活動には頭が下がるけれど、そういう活動に少なからず見られる「見下し感」がどうしても好きになれない。その「トンデモ理論」だけがその人の心の拠り所だったり、それが商売に直結していたりして、そうすると彼らは必死なんだ。それを抜きにして、ただ基礎理論の不勉強を馬鹿にするのは、あまり実り多い活動とは思えない。

 人間は、論文なりソフトウェアなりを介して相互作用をするのがちょうど良い。プログラマはプログラムで、研究者は論文で発言すればそれで良い。

 昨日来てた査読依頼、ちゃんと概要を読んでみると全く専門外。なので普通に断ろうと思ったが、なんとなくエディタの名前に心当たりがある(特徴的な名前だった)。で、過去の査読履歴調べてみると、これが二度目の依頼で、その時にも「専門外。以前ISSPにいた別のHiroshi Watanabeと間違えているのではないか?」と回答している。にもかかわらずまた来たので、今度は「専門外だって言ったでしょ!データベースから僕のアドレス削除しておいて!」と(もちろんソフトに)伝えておいた。

 想像するにこのエディタ、僕の「人違いでは?」という内容は無視して、単に別の候補に回したのだろう。今回も無視されそうな気がする。そうしたらまた回ってきそうな気がする。


6月15日(水)

 NOP

 千葉県民の日


6月14日(火)

 式変形の続き。分布関数は無次元なので、スケール変換しても(ヤコビアン以外に)ゴミはつかないが、 Kinetic Termは次元があるので、そのままではスケール変換でゴミがでてくる。 ある自由変数を固定して考える限りゴミは無視してもなんの問題もないのだが、 別の自由変数に取り直した時には、ゴミまでちゃんと考えないと変換がおかしくなる。 どういう表示を取ろうと全部コンシステントに計算はできるのだが、出てくる量を全て無次元化しておくと、ちゃんときれいで物理的に妥当な変換式が成り立つようになる。

 あと、理論構築では基本無視する比例係数だが、数値データと比較する際にはきちっと全部追わないといけない。 途中でπが出てきたりx乗されたりと、追うのが極めて面倒くさいが、やらないといけない。 有名な理論家であるK教授がようやく論文を書けるようになったのは、√2やπのようなクズどもを省略せずに計算をするようになってからだそうな。

 まぁとにかく混乱したが、収束しつつある。細かい計算の導出まで含めたノートを作ったので、論文再投稿したら公開しよう。ニーズあるかわからんけど。

 萱沼先生の定義集といえば、「D論は(物理的に)立つべきである」という言説を思い出す。以下引用。

筆者の前任校だったT北大学は重厚長大を重んずる学風だった。とりわけそれを感じたのは、ある院生の学位論文の副査になって審査委員会に臨んだときである。主査のE先生は「立たないじゃないか」と大いにご不満だった。長年の研鑽の結果である博士論文なのだから、机に立つくらい分厚いものであるべきだ、というご意見である。

 確か学生時代に助手のYさんに聞いたのだったと思うが、某大学には同様な思想の「磁石テスト」というものがあった。修士論文を印刷し、ホワイトボードに磁石一つで止めて落ちてこなかったら失格、というものである。ここで「両面印刷なのか片面印刷なのか」とか「紙質の指定はどうなのか」とかレギュレーションをどうこう言わないのが大人である。

 この話をした時、「修士論文を磁石一つで止めるということは、その磁石ってネオジムですよね。ネオジムはキュリー点が低いので、磁石テストの前に加熱しておけばいいんじゃないですか?」と言い出したのは当時川島研M2のK君であったが、彼はもうこの発言だけでも修士の学位に値すると思った。

 ちなみによく使われる磁石は、ネオジム(300度)、アルニコ(860度)、フェライト(450度)、サマリウムコバルト(750度)などがある。カッコ内はおよそのキュリー点(摂氏)。確かにネオジムは熱に弱い。


6月13日(月)

 ずっと計算が合わない問題、原因がわかった。ある等価な問題を別の変数で表現すると、一致すべき結果が一致しないことに悩んでいた。どちらの変数でも最終的に常微分方程式に落ちるのだが、同じ物理を表現しているのに、その最終的な結果が変数の取り方に依存してしまう。

 で、その原因は、微分方程式にあらわれる関数が分布関数であることを忘れていたから。 分布関数は、変数変換でヤコビアンがゴミとして出てくる。それを完全に失念していた。 言い訳すると、もともと偏微分方程式だったものを、変数をスケール変換することで常微分方程式に落とすのだが、 その計算がちょっとややこしくて、ある関数が「もともと分布関数であった」ということをすっかり忘れ、ただの常微分方程式としての変数変換だけを考えていた。これでそうとう時間を無駄にした。アホ過ぎる・・・。


6月12日(日)

 査読返した。

 ずっと悩んでいた問題、わかった気がする。


6月11日(土)

 NOP


6月10日(金)

 5月30日に行われたトークイベント「岩波データサイエンスⅡ ナンプレと魔方陣 解いたり作ったり数えたり」の様子を収めた動画が岩波データサイエンスサポートページに公開されました。 僕の話はこのあたり (@YouTube)から。

 CMSI配信講義の時も思っけど、僕ってしゃべるときに左右に揺れる癖があるな。完全に無自覚だった。動画で見るとかなり気になるので次から気をつけよう。あと、同じ話を何度録画&公開されても別に気にしないけど、ウケを取りに行くネタを録画されてしまうと厳しい。そんなにネタがあるわけでは無いので。

 このあたり、お笑い芸人とかどうなんだろ?コントとか漫才とかってそんな量産はできないはずだし、わりと動画サイトにアップロードされちゃったりしてるよね・・・(←立ち位置が完全にそっち寄り)。

 あ、そうそう、この講演ではいくつか不正確なことを言っているので注意。途中で紹介した「数独のサイズと解の数と最小ヒント数の関係に関するconjecture」は、より大きなサイズでは反例が見つかっているそうです。

 ちなみに岩波データサイエンスのサポートページはかなり充実しているので、本を持っていない人にもおすすめ。

 YouTubeと言えば、CMSI計算科学技術特論Cの講演動画もアップロードされてた。 これはわりと自分では意識的に芸風を変え、テンション低めに話したもの。自分では気に入っているが、公開後4ヶ月で視聴回数が21回とか・・・。

 その、途中のプログラム論はともかく、闘うプログラマーの感想あたりからは(自分で言うのもなんだけど)ちょっとおすすめというか、多くの人に聞いて欲しいというか、まぁ要するに愚痴みたいなもんですね。


6月9日(木)

 「FortranとCのコンパイラの最適化性能の違い」についてベンダーに投げたら、回答がきた。 force.cppについては、コンパイルオプション「-Kswp_strong」をつけることでforce.fと同等の性能になる。

「swp_strong!、そういうのもあるのか」

 なぜC++とFortranでコンパイラの解釈が異なるか、force.(f/cpp)とforce2.(f/cpp)の最適化内容の違いについても説明を受けた。こんなに早く回答が来るとは思ってなかったので驚いた。

 国産コンパイラだと、こういう対応の早さというメリットがあったりするなぁ。以前Intelコンパイラのバグを報告した時、その時の対応はえらく遅かった気がする。


6月8日(水)

 なんかもうなーんもできないうちに日々が過ぎていくぞ。

 この前のナンプレイベントで「数独の最小ヒントが17であることを示すのに、81C17を全探索するのは絶望的」という話に対して、 「81C17はどのくらいか?」という質問がきた。統計物理学者なら暗算できてしかるべきだが、とっさのオーダーもわからなかったので もういちど見積もってみよう。なお、実際に行われたのは「16ヒント問題全てが解が一意でないことの確認」なので、その文脈では81C16を見積もるべきだが、それはさておく。

 81C17 = 81!/(64! 17!)だが、とりあえずスターリングの公式を使おう。 まず 81 ln 81を見積もる。e〜3と近似すると81=3**4なので、 81 ln 81 〜 324。 次、64 ln 64だが、27<64<81なので、ざっくり64〜 3**3.5と近似しよう。すると 64 ln 64 〜 64 * 3.5 = 224。 17 ln 17だが、9 <17< 27なので、これも3**2.5と近似するか。17 ln 17 〜 17 * 2.5 = 42.5、四捨五入して43。 引き算すると、ln (81C17) 〜 324-224-43 = 57か。

 次に底の変換をしないといけない。e**57は10の何乗か。e=3と近似しているので、10=e**2と近似しちゃうか。 するとe**57 〜 10**28.5。答えは10**17.1なので、全然ダメですね。

 もうちょっとなんとかするか。ln(2) 〜 0.7、ln(10) 〜 2.3であることを知っていたとしよう。 また、ln(3)〜1はそのままにする。すると、

81 ln 81 〜 324
64 ln 64 〜 64*6*0.7 = 269
17 ln 17 〜 17 ln 18 = 17 (2 ln 3 + ln 2) = 17*2.7 〜46
ln (81C17) 〜 324-269-46 = 9
81C17 〜 10**4

 ありゃー、全然ダメだ。よく考えると81 ln 81ではe=3にしているのに、64 ln 64で正しいeを使ってるからコンシステントな近似になってないな。実際には

81 ln 81 〜 356
64 ln 64 〜 266
17 ln 17 〜 48
ln (81C17) 〜 356-266-48 = 42

 となり、これから計算すると42/2.3〜18となって、17とまぁまぁ合う。


6月7日(火)

 なにしてたっけ?確か査読の返事のための式変形をしてた気がする。また曜日とか日付とか間違ってるし。


6月6日(月)

 え?もうナンプレイベントから一週間たったの?


6月5日(日)

 NOP


6月4日(土)

 NOP


6月3日(金)

 論文の改訂をするために、レフェリーに指摘された大量の論文の読み込みをやっているのだが、なかなかすすまない。

 そんな中、以前査読した論文の返事が返ってきた。最近、査読した論文がRejectされることがおおくて(そういう論文ばっか来るからだけど)、ちゃんとした論文で、ちゃんとした返事が来るの久しぶり。著者から返事が来ると、もう一人の査読者のレポートが読めるんだよね。Rejectされちゃうと、もう一人のレポートが読めない。で、今回のもう一人の査読者のレポート見て「そうそう!そうだよね!そう思うよね!」とちょっとうれしくなってた。いや、今回の論文、内容はさほど問題ないんだけど、文章が極めて冗長で、あまり論文に関係ない文章が多かったから、「削れ!」って言ったんだけど、もう一人もそう言ってた。

 ・・・もう少し「もう一人の査読レポート」読んでたら、そのあまりに厳しい言葉にちょっとびっくりした。いや、たしかに僕もそう思ったけどさ、そんなに怒らなくても・・・。

 思うに、「もう一人の査読レポート」ってかなり貴重な気がする。同じ論文を詳細に読んだ感想をレポートにまとめるなんて、ゼミでもあまりやらない。

 いやでも厳しい。自分の論文の査読レポートの返事がかけてないのに、他人の論文の査読レポートの返事を読まないといけない。仕事が遅いのがいけないんだけどさ・・・。

 しばらく前からおかしくなったGreen 500のページ。トップページがエラー表示になってる。Drupalで作られてたのか。MySQLのエラーが出てたのが、Offline表示に変えられている。 で、/includesを見ると何も表示されないんだけど、HTMLソースにはコメントで「Placeholder file for previously hacked file that contained string HackeD By TiGER-M@TE」とある。ググってでてくる限りでは、バングラデシュのハッカーで、Googleバングラデシュ版をハックしたことがあるらしい→「Exclusive Interview with TiGER-M@TE (Bangladesh Google website Hacker)」。でも今調べても、TiGER-M@TEとGreen500を関連付ける情報は一つも出てこないなぁ。なんだろう?

 とりあえず、/includesにあったこの文字列が気になる。ここに置いておくので、誰か解読してくれません?

 論文だいぶ読んだ。査読レポートへの返事の方針は固まった気がする。


6月2日(木)

 大変不調。昨日の曜日間違ってるじゃん。これだから手書きは・・・。


6月1日(水)

 FortranとCのコンパイラの最適化性能の違いについて。いろんなところでよく言われていることだけれど、具体例を挙げる。 なんどか例に挙げているが、こんなコードを書く。

//----------------------------------------------------------------------
#include <stdio.h>
//----------------------------------------------------------------------
const int N = 20000;
const int D = 3;
const int X = 0;
const int Y = 1;
const int Z = 2;
//---------------------------------------------------------------------
double q[N][D],p[N][D];
//---------------------------------------------------------------------
void
calcforce(void){
  for(int i=0;i<N-1;i++){
    for(int j=i+1;j<N;j++){
      double dx = q[j][X] - q[i][X];
      double dy = q[j][Y] - q[i][Y];
      double dz = q[j][Z] - q[i][Z];
      double df = (dx*dx + dy*dy + dz*dz);
      p[i][X] += df*dx;
      p[i][Y] += df*dy;
      p[i][Z] += df*dz;
      p[j][X] -= df*dx;
      p[j][Y] -= df*dy;
      p[j][Z] -= df*dz;
    }
  }
}
//----------------------------------------------------------------------
int
main(void){
  for(int i=0;i<N;i++){
    q[i][X] = i;
    q[i][Y] = i;
    q[i][Z] = i;
    p[i][X] = 0.0;
    p[i][Y] = 0.0;
    p[i][Z] = 0.0;
  }
  calcforce();
}
//----------------------------------------------------------------------

 これは力の二重ループがソフトウェアパイプライニングされず、性能が出ない。実行すると手元のFX10で2.08sかかる。

 これを「そのまま」Fortranで書きなおそう。

       subroutine calcforce(q,p)
       implicit none
       integer,parameter:: N=20000
       integer,parameter:: D=3
       integer,parameter:: X=1,Y=2,Z=3
       double precision:: q(D,N), p(D,N)
       integer :: i,j
       double precision:: dx, dy, dz, df
       do i=1,N-1
         do j=i+1,N
           dx = q(X,j) - q(X,i)
           dy = q(Y,j) - q(Y,i)
           dz = q(Z,j) - q(Z,i)
           df = dx*dx + dy*dy + dz*dz
           p(X,i) = p(X,i) + df*dx
           p(Y,i) = p(Y,i) + df*dy
           p(Z,i) = p(Z,i) + df*dz
           p(X,j) = p(X,j) - df*dx
           p(Y,j) = p(Y,j) - df*dy
           p(Z,j) = p(Z,j) - df*dz
         enddo
       enddo
       end

       program force
       integer,parameter:: N=20000
       integer,parameter:: D=3
       integer,parameter:: X=1,Y=2,Z=3
       double precision:: q(D,N), p(D,N)
       integer :: i
       do i = 1, N
         p(X,i) = 0.0
         p(Y,i) = 0.0
         p(Z,i) = 0.0
         q(X,i) = i
         q(Y,i) = i
         q(Z,i) = i
       enddo
       call calcforce(q,p)
       end program force

 計算ロジックは全く同じ(CとFortranでは多次元配列の順序が違うので入れ替えてある)。これをFortranコンパイラでコンパイルすると、こちらはソフトウェアパイプライニングされる。計算時間は 1.65s。つまり、同じコードをFortranで書きなおすだけで26%高速化される。

 「じゃ、Fortranで書けば?」と言いたくなるのはグッとこらえてもらって、もう少し最適化(という言葉が適切かどうかはともかく)する。力の計算で、内側のループで値が変わらないi粒子の座標を外に出し、かつi粒子の運動量を毎回書き戻さないように外に出そう。

void
calcforce(void){
  for(int i=0;i<N-1;i++){
    const double qx = q[i][X];
    const double qy = q[i][Y];
    const double qz = q[i][Z];
    double px = p[i][X];
    double py = p[i][Y];
    double pz = p[i][Z];
    for(int j=i+1;j<N;j++){
      double dx = q[j][X] - qx;
      double dy = q[j][Y] - qy;
      double dz = q[j][Z] - qz;
      double df = (dx*dx + dy*dy + dz*dz);
      px += df*dx;
      py += df*dy;
      pz += df*dz;
      p[j][X] -= df*dx;
      p[j][Y] -= df*dy;
      p[j][Z] -= df*dz;
    }
    p[i][X] = px;
    p[i][Y] = py;
    p[i][Z] = pz;
  }
}

 こちらはソフトウェアパイプライニングされた上、総和演算の最適化もかかって、実行時間は1.21sになり、速度がFortranコードを超える。 Fortranコードも同様な最適化をかけてみよう。

       subroutine calcforce(q,p)
       implicit none
       integer,parameter:: N=20000
       integer,parameter:: D=3
       integer,parameter:: X=1,Y=2,Z=3
       double precision:: q(D,N), p(D,N)
       integer :: i,j
       double precision:: dx, dy, dz, df
       double precision:: qx,qy,qz,px,py,pz
       do i=1,N-1
           qx = q(X,i)
           qy = q(Y,i)
           qz = q(Z,i)
           px = p(X,i)
           py = p(Y,i)
           pz = p(Z,i)
         do j=i+1,N
           dx = q(X,j) - qx
           dy = q(Y,j) - qy
           dz = q(Z,j) - qz
           df = dx*dx + dy*dy + dz*dz
           px = px + df*dx
           py = py + df*dy
           pz = pz + df*dz
           p(X,j) = p(X,j) - df*dx
           p(Y,j) = p(Y,j) - df*dy
           p(Z,j) = p(Z,j) - df*dz
         enddo
         p(X,i) = px
         p(Y,i) = py
         p(Z,i) = pz
       enddo
       end

 同じことをやっているのだから当然なのだが、同じ最適化がかかって1.20sと、ほぼ同じ実行時間になる。Cで書いたオリジナルコードにたいして70%の高速化、Fortranのオリジナルコードに比べても38%の高速化となる。

 まとめると「同じコードをFortranで書き直すだけでかなり高速化される」という事実はある。だけど、ちゃんと書きなおせばどちらも同じ速度になる。どこまでチューニングするかはもちろん場合によるけれど、Fortranで全部書き直す手間と、なぜ性能が出ないかを確認してCのコードを直す手間はどっちがどうかなぁ、という気もする。少なくとも今回の例ではFortranで書き直すだけでは最適化は明らかに不十分なので、「CコンパイラがダメならFortranで書けばいいじゃん」は一種の思考停止だと思う。Fortranで書いたコードだって、チューニングするわけだし、ちゃんとチューニングした際に到達し得る性能は言語に関係ないはずだから。

 もちろん、「自分は最適化しない」という立場を取ることも考えられる。その場合は言語選択としてFortranを使う、というのは合理性がある。

 まぁ、文句言うばかりじゃなくて、ちゃんとベンダーに報告もしている。またこのコンパイラとお付き合いすることになりそうだし・・・。


トップページへ戻る