ToDo:
MADX->SAD変換した latticeを SAD->MADX変換した際に、vertical dumpの 符号が反転する謎現象が発生。
原因は、vkicker -> rotated BEND -> tilted hkickerと変換するのだが、 マニュアルの記載通りだとtilt = π/2な hkickerは vkickerになる はずなのだが、MADXの twissモジュールが計算するCODは符号が逆になる。
念のため、MADX trackモジュールをテストすると hkicker + tilt = π/2は vkickerと等価で、vkicker + tilt = π/2は -hkickerと等価になり こちらはマニュアル通りに動いているように見える。
MADXのtwissモジュールのバグか？
LHCの latticeを MADX->SAD->MADXと round trip変換したものを MADXで計算した結果が、オリジナルとほぼ一致するところまでは 実装出来た。(エレメントの情報量が全単射な変換ではないので round trip変換が恒等写像にならないため)
KEKB latticeの変換はアーク部に関しては、quadrupoleの fringeが MADX側に無いために betatron tuneがずれるようである。
IR部を含んだ変換は、twissパラメータの計算が収束しない模様
おそらく、DX/DYの処理と K0/SK0で軌道を変えずに dipole fieldを 入れるための小細工が必要と思われる。(MADXの multipoleの 0次項は bendと同様にデザイン軌道を曲げるらしい)
QUADのF1/F2サポートは、quadrupoleの前後にmatrix(2次までの写像を 記述できるエレメント)を貼り付けて、SADのlinear fringeの項を 加えるしかないか...
まあ、どんなSAD latticeでも変換できる変換エンジンとしては
2と3の方法は、本質的にはMADX互換のインターフェースを持ったSADの オプティクスモデルである。 多分、MADX版の latticeが欲しいと言う人間は、SADとMADXの信頼性や 精度を比較して MADXを使いたいと言っているのではなく、単に MADXしか使えないからだと思われるので、それなら UI部が MAD互換で ありさえすれば中身のコードは別に MAD由来でなくても良いはずだよね？
少なくとも、Makefileで使う perlスクリプトが source repositoryに 含まれて無かったり、テストセットが付属していないという意味では MADXがSADよりも良くテストされているとは思えないのだが...
コードの規模が小さいと言う意味では、本格的なテストを実施しやすい 素養はあるが、定期的に実施し結果を公開していなければ意味が無いし (外部ユーザーから見た信頼性という意味で)、コードの小ささは 機能の少なさの裏返しだし...
MADXのエレメントパラメータに一貫性が無いのはなぜだろう？ rbend/quadrupole等ではlength = 0が使えないのに、multipoleは lengthを指定できないとか... Orz
PEGのメモの状態やルールの適用順をダンプしながら調べることで 左再帰の実装のバグを特定＆修正出来た。
左再帰マークの付いたルールが、左再帰のループの根元のルールが 終了した後に再評価される際に、既に存在しない根元のルールを 探して再帰スタックをたどってしまうのが原因だった。 左再帰実装の元ネタ論文では考慮されていないケースと思われる。
おそらく、文法の書き方によっては避けられる(左再帰に含まれるルール集合 の上位ルールに唯一性が成り立つ場合とか)と思うが、パーサーの性能で 文法記述が制約されるのは気にくわないので修正。 左再帰の根元となるルールが適用中かどうかを記録するフラグを増設して、 根元のルールの適用が終了後に左再帰マークの付いたメモを見つけたら メモから削除するようにして解決した。
修正済みのパーサーで、Map演算子の左オペランドが Function演算子や Unset演算子で終端される構文が正しく解析できるようになった。
構文解析の進捗解析用に生成されるログは、約113MB、285万行に及んだ。 もう少し文法記述を改良してダイエットしたほうが良いのかもしれない。 もちろん、正しく動作する文法記述を完成させてからの話ではあるが...
a=.b=cのパースで再帰がうまくいかないようで、 (Unset a)までしか読めていない Orz
a=.*b=cは(Set (Times (Unset a) b) c)へパース出来るので、 可換積(Times)ルールの演算子省略形の定義が不味い？
1k/10k particlesの結果が出た。
肝心の ramping periodに対する依存性だが10ターンを越えると依存性が見えない
次は、crab cavity voltageへの依存性を見てみる
QUAD->quadrupole変換時にSADのtracking engineから linear fringeの写像から2次まで項を展開して MADXのmatrixエレメントとして挿入してみた。
KEKB arc cellを変換した際の betatron tuneの一致度がかなり改善して、 SADモデルの値をほぼ再現できた。
KEKBの変換結果に関しては、
| Ring | Section | Status |
| HER | Arc | OK |
| HER | 富士 | OK |
| HER | 大穂 | OK |
| HER | 日光 | OK |
| HER | 筑波 | NG |
| LER | Arc | OK |
| LER | 富士 | OK |
| LER | 大穂 | NG |
| LER | 日光 | NG |
| LER | 筑波 | NG |
筑波直線部に関しては、MULTエレメントの DX/DYサポートとSOLエレメントの近似を 入れないと動かないのは分かるのだが、LERの日光と大穂が なぜ変換できないかが謎である。 Bendの集団で記述している Wigglerが MADXと相性が悪いのか？
日曜日にinjection studyの様子を見に行ったコントロールルームで Rama達とデバックした結果をまとめておく
sbendのマニュアルにはk0パラメータは無視される的な記述があるので design orbitを曲げずにビームを蹴るエレメントがkickerか multipoleぐらいしか無いのだが、thickな記述はkickerでしか出来な い Orz
まともに変換出来るようにするたびに、エレメントの分割と matrixエレメントの挿入が増えている...
MADX変数経由のパラメータ埋め込みにすれば、一応マッチング等は可能だと 思うが一部に超越関数が含まれているので、この部分は MADXの実行エンジン にネイティブに実装されて無いとかなり無理っぽい(指数関数とか三角関数系 なので原理的にはマクローリン展開とかで書けるけど...)
こうした変換を行った latticeにどこまで価値があるかどうかは 正直疑問である。
getpgid/setpgid(2)システムコールへのインターフェースを実装
fork(2)での並列化を行うときに効率よく forkした子プロセスを待つための Process Group IDの設定と確認が行えるようになりました。
CVS MAIN trunkに入った汎用性が低そうなParallelize[]に対抗して、 ParallelMap[]の開発を開始したが、現時点でのプロトタイプによる 2 CPUノードでのベンチマークで、シングルの 85%増しの速度しか出ていないOrz
プロトタイプだけに性能に関するチューニングはしていないが、それ以上に 例外処理とかMap[]の階層指定子のサポートなどのコードが抜けているので プロダクション仕様だともう少し遅くなりそうな
ノード数に対するスケーラビリティが出ない理由は、推測の域を出ないが 評価関数自体の実行速度が落ちてる＆並列版ではMap[]での後半の評価ほど 遅くなる辺りから、fork(2)した後に評価関数を演算する際に発生する ヒープとインタプリタスタックをCopy on Writeで更新するコストが 高いように思える。後半の処理でヒープ等のフラッグメント化が 進行しているとすれば説明が付く。
後で、実行順序を制御したベンチマークで検証してみる。
フラグメント化が原因だとすると、リストの要素毎に fork()せずに 子プロセスを使いまわす工夫が必要になると思われる。
プロトタイプ実装した際の並列制御の実装メモと問題点への考察をまとめておく
Frok[]した子プロセスで評価した値をShared[]経由で親プロセスが回収。 子プロセスとの同期には、SetPGID[]で割り当てたプロセスグループIDを 使ってブロッキングモードのWait4[]で待ち受ける。(子プロセス達の どれかが計算終了するまで親プロセスは休眠する) 子プロセスを回収後に新たな子プロセスをFrok[]する。
Shared[]の使いまわし可。
制御は比較的簡単に実装できるが、評価毎に Fork[]のコストが必要。
Fork[]した子プロセスとShared[]経由で評価すべき内容と結果をやりとりする。 Fork[]した時点で、評価関数と評価対象のリスト全体が複製されているので 評価の指示は、リスト上の位置情報で十分である。 プロセス間の相互の同期は、SYSV IPCや POSIX Semaphoreへの インターフェースが未実装なので、Pipe[](pipe(2))と SelectUnit(select(2))を使って同期を実装する。
子プロセスは親プロセスからパイプに合図が送られてくるのを selectで ブロッキングしながら待機し、合図が来たら Shared[]から評価要求を 受け取り、評価結果を Shared[]へ書き込み親プロセスへのパイプに 1byteの合図を書き込んで待機状態に戻る。
親プロセスは、旗下の子プロセス達のパイプから合図がくるのを selectで ブロッキングしながら待機し、合図が来たら Shared[]から評価結果を 受け取り、次の評価要求を Shared[]へ書き込んだ後、子プロセスへのパイプに 1byteの合図を書き込んで待機状態に戻る。
エラーハンドリングとしては、子プロセスが異常終了した場合に備えて、 ノンブロッキングモードのWait4[]で子プロセスを指定して生存を確認し、 死亡している場合は、パイプを再生成して Fork[]する。
Flush[]と Close[]の動作に要注意！ パイプをFlush[]すると LFが送られることがあるので、 Flush[]の使用は避けること！
また、Close[]が、パイプを閉じるときに LFを送ることがある。 したがって、受信側のプロセスが死亡する前に送信側のパイプを Close[]しておかないと、Close[]実行時のLF送出に失敗して Broken pipeでアボートしてしまう。 (おそらく、Fortran I/Oライブラリ内で発生している)
数値リテラルの扱いをGMP/MPFRを使って実装しようと考えているのだが、 数値クラスを扱いやすくするために、C++上のクラスと 演算子オーバーロードを定義する方法を考察中。
単純な実装では、a=b+cを実行するにはb+cの結果を格納する 一時オブジェクトfooを生成してaに対する代入演算子の オペランドとして使用後に破棄する必要があるが、この一時オブジェクトが 多倍長精度整数や多倍長精度浮動小数点であった場合、生成破壊のコストが 高く付く(式が複雑な場合は特に)訳で、性能のボトルネックになるのは 確実である。既存の実装はどのようにしてそれを避けているかを調べてみた。
MPFR自身にはC++サポートは存在しないが、GMPにはC++サポートがある。 で、ヘッダーを読んでみて目から鱗が...
結論は、演算子が式を表現するProxyオブジェクトを返すことで、 実際の式評価を評価結果の格納先が確保されるまで遅延する(遅延評価)ことで、 全ての演算を3オペランド形式で実行し、余分な一時オブジェクトの 生成コストを回避していた。 GMPの実装では、演算子表現に関数オブジェクト、式のProxyオブジェクトは テンプレートクラスから実体化している。また、一部の演算は テンプレートクラスの特殊化を行い、演算量を削減している。
このテクニック自体は、ベクトルやマトリクス演算子にも応用できるはず。
実装上の注意としては、次の点があげられる。
遅延評価なので、最終的な結果が代入されない場合は Proxyオブジェクトの 生成破壊のみで、演算評価を行わないですませることが出来る
ただ、コンパイラの仕事は大変そうだ(w
フロントエンドラッパーともっとも単純なfork/wait4ベースの並列実装 (子プロセスを評価毎に生成して使い捨てる)をamorita branchの develに commitした。
Library@Require["Algorism/Parallel"]; (* ライブラリ要求 *) $ParallelMethod = "Fork"; (* デフォルトの並列化実装を Forkに設定 *) l = GaussRandom[100, 100]; result = Parallel@Map[Fourier, l, SharedSize->1024, Multiplicity->5];
通信バッファサイズを1024bytes、並列度を5で、 リストlの各要素にFourier[]並列に作用させる。
デフォルトオプションは、{SharedSize->256,Multiplicity->NPARA}
現時点の実装では、階数1へのMapのみをサポート (つまり、階数指定子無しのMapやMap演算子と同じ)
並列にすれば早くなる訳ではなく、 リストの各要素に作用させる評価関数が十分に重い＆返却データのコピーコストが十分に小さい等の 条件が成り立たないと早くならない。 また、fork/wait4システムコールの呼び出し負荷やスケジュラーへの負荷が あるので、混雑したシステム上での過度な並列化は全体の性能を低下させる。
(1)に関しては、部分的に実装済みで、子プロセスが共有メモリ経由で 結果を返送出来ない場合は必要な共有メモリサイズを親プロセスに通知、 親プロセスは受け取れなかった結果を自分で計算すると同時に 以降の通信バッファのサイズを通知に基づいて更新する。
最悪のケースは必要な通信バッファの量が評価順序にそって増えた場合で、 このケースは救えない。
結果を回収するには別チャンネルで、子プロセス-親プロセス間の通信を 行う必要がある。案としては...
(1)は、通信に際して親プロセスと子プロセスの相互同期が必要であり、 データエンコードスキームとしては Shared[]が内部的に使っている uni-byte stream encoder/decoderの相当品が必要 (ToString[]によるASCII変換は、ToExpression[ToString[#]]&が 恒等変換で無いので使用できない)であり、実装し難い。 (2)は、ファイルに格納するデータエンコードに(1)と同様の問題有り。 (fork/wait4モデルでは、ファイルシステムの名前空間の断絶を考慮する必要はない) (3)は、現状のOpenShared[]が無名の共有メモリしか扱えないので無理。
多分、性能面での筋が良さそうなのは、OpenShared[]とそのバックエンドとなる mapalloc内部APIを拡張して、ファイル上に共有メモリを確保 出来るようにして、共有メモリとして使う一時ファイル名を受け渡すのが 妥当と思われる。
別の応用(数値計算結果の正確な保存と復元)を考えると、Shared[]が内部的に 持っているuni-byte stream encoder/decoderの外部APIを提供することも 意味があると思われる。
昨晩、後からプロセス間共有メモリを作る方法を調査したのでメモメモ
open/mmapでファイルバックングストアを持った共有メモリをプロセスに貼れる。 問題は、共有メモリに対する排他制御。 SADに実装されている方式は、multi-coreシステムでは動作が保証できないと いうか、一般にメモリバリアもしくはアトミック操作 (read-modifier-write/compare-swap)無しに、不特定多数のプロセスからの クリティカルセクションは作れない。(特定N者間の排他制御を アトミック操作無しのスピンロックで実現するアルゴリズムは存在している)
問題は、mutexや semaphoreを共有メモリ内に生成するか、 カーネル内に生成する必要があるところ。
pthread_mutex_init(POSIX mutex)やsem_init(POSIX semaphore)は、 この時点で脱落確定。(実態の生成先を制御出来ない)
可能性としては、SYSV IPC semaphoreかsem_open(POSIX named semaphore)辺り が候補として残っている。 どちらにせよ、共有メモリ開放時に同期資源も開放する必要があるので、 共有メモリは自分自身への参照数と同期資源への参照を管理する必要がある。 ファイルバッキングストア経由で、attachすることを考えると header magic と check digitを含める必要がありそう。 (socket経由で file descripterを転送するという手もあるはずだが、 動作が分かりづらい＆共有メモリをattachするまでのコストが高そうだ)
と言う訳で、ファイルとして外部化した共有メモリが含むべき情報は
ということになる。
現在主力のGCC4.3は、Fortran2003ベースで 2008で導入される数学関数は 存在しない。しかし、hypotとかのようにあった方がコードが 意味論的にすっきりする数学関数がFortran2008では導入されている。 と言う訳で、増設された数学関数の内、実質的にユーザー定義可能なもの (COMPLEX/REALで多相化されていないもの)に関しては、 src/sim/fortran_math_foo_.cで実装を作って使う方向でコードを 整理してみる。 (hypot(x,y)を abs(dcmpx(x,y))で書いてあるのは読みにくいと思う)
制御変数FORTRAN_MISSING_MATH_FUNCSと 環境依存Fortranインクルードファイルinc/MATH.incを導入した。
Fortran2008で$\Gamma$関数と対数$\Gamma$関数が標準に入ったので、 使おうと思って調べていたのだが、最終的には libmを呼び出す訳で 呼び出し先のコードを追跡していたら... gamma(3m)は$\Gamma$関数ではなく、 対数$\Gamma$関数の値を返すんだね... Orz 何のために、lgamma(3m)が有るんだか src/lib/msum/src/e_gamma_r.cの __ieee754_gamma_rは$\Gamma$関数の対数を返すと言う コメントを読んでる当たりで、頭痛がしてきた。
$\Gamma$関数を返すのは、tgamma(3m)です。
はい、皆さん準備はいいですか？
右見て、左見て、半顔、見上げて...
絶望した！libmの紛らわしい仕様に絶望した！
本家のstable-MLでも悲鳴が挙がっていましたが、 DTraceのMFC直後のRELENG_7ツリーはコンパイル出来るけど、 再起動すると signal 11でPID 1のプロセスの立ち上げに 失敗して悲しいことになります。 (つまり、シングルユーザーモードも動かない)
取り合えず、Dual bootの WindowsでFreeSBIE 2.0.1を焼いて、 FreeSBIEでブートしてバックアップを置いてある外付けUSB HDD (中身は本体の 2.5inchと同じもの)から、 /bin /sbin /usr/bin /usr/sbin /usr/lib /usr/libexec辺りを 持ってきて仮復旧、cvs updateで修正済みツリーに更新して make buildworld buildkernelして再度更新かけました。
注意点は、FreeSBIE 2.0.xは FreeBSD 6.2ベースなので、 FreeBSD 7.0-STABLE系の正常なバックアップが無いと復旧が面倒。 バックアップを持ってないときは、FreeBSD 7.0の liveCDを 立ち上げてそこからシステムをコピーするのが早そうですね。
本家のstable-ML見る限り ZFSもデグレードしてるっぽいので、 使ってる人は要注意
一応、調査した結果をまとめておく
要するにSADとMADXのモデリングとバグのまめ知識
| MADX | SAD |
| betx, bety | BX, BY |
| alfx, alfy | AX, AY |
| dx, dy | PEX, PEY |
| dpx, dpy | PEPX, PEPY |
| x, y | DX, DY |
| px, py | DPX, DPX |
dlfunaloc APIとSAD_STATIC_DEF_FUNCS経由で機能単位の 組み込みを行っている部分に関して、機能毎での追加や禁止を行う フレームワークを導入。
どちらも、sad.confやconfig/*.specで行えます。
特定環境で不要、またはビルドできない機能を切り離したり、 特定環境のみに有用な機能を追加するのに使えます。
この機能は、ブランチ側での関連する修正が落ち着いたら、 CVS MAIN trunkへバックポートする予定
カテゴリー: Admin | Emacs | EPICS | Fortran | FreeBSD | GCC | hgsubversion | IPv6 | KEKB | LHC | Lisp | LLVM | MADX | Ryzen | SAD | samba | tDiary | unix | WWW | YaSAI | お仕事 | イベント | 出張 | 宴会 | 数学 | 艦これ | 買いもの | 追記 | 雑記