T2K と E-Simulator による高層ビルの大規模非線形耐震解析
秋葉博(1)，橋詰和明(1)，宮村倫司(2)(1)株式会社アライドエンジニアリング(2)日本大学工学部，防災科学技術研究所

1. はじめに
著者らは，2007 年から，独立行政法人防災科学技術研究所・兵庫耐震工学センター (Eデ ィ フ ェ ン ス )([1]) の 研 究 の 一 環 と し て ， 数 値 震 動 台 開 発 分 科 会 の 下 ， 数 値 震 動 台(E-Simulator) の開発に携わってきた。
E-ディフェンスは，世界最大規模の震動台施設 (実大三次元震動破壊施設) を有している。
2005 年以来，実物大の戸建建物や低層ビルなどの震動・破壊実験が多数行われており，実用に供してきた実績がある。しかしながら，震動台実験では，実大実験の規模や破壊過程の細部の評価にも自ずと限界があり，また，1 つの実験は，実物のビルを建設し，震動台で加震・破壊し，これを解体して撤去するという，コストも時間もきわめて大がかりなものになる。E-ディフェンスの目標の一つは，膨大な実験データの蓄積をもとに，実大震動台を補完し，より大規模な建設構造物の震動・破壊シミュレーションを行うための「数値震動台」の開発である。
E-Simulator のフレームワークは商用構造解析コード ADVENTURECluster (以下ADVC と記す)([2]) である。ADVC は ADVENTURE プロジェクト([3])の ADVENTUREシステムを元に作られており，並列処理アルゴリズムの主要部やデータ構造はADVENTURE のものを引き継いでいる。ADVC は，その上で汎用構造解析コードとしての幅広い機能を備えており，国内では自動車，エレクトロニクス，エネルギ，精密機器メーカなど多くのユーザを有しているほか，中国の主要な企業にも輸出されている。
E-Simulator には，さらにこの上にコンクリートの構成則，精密な複合硬化を考慮できる鋼材の構成則，地盤モデル，コンクリートのクラックなどを表現できる破壊モデルなどが追加されている。ただし，本稿で述べる E-Simulator の解析ではこれら拡張機能を使っていないため，ADVC と読み替えても同じことがいえる。
T2K のインフラと E-Simulator (ADVC) のパフォーマンスが力を発揮するのは巨大モデルの解析である。特殊な構成則などの紹介は他の文献に譲り([4]～[11])，本稿では T2Kに お け る ソ ル バ の 並 列 パ フ ォ ー マ ン ス に つ い て 述 べ る 。 ADVC の 並 列 ソ ル バ は
ADVENTURE システムや他の並列構造解析システムが採用している，いわゆる反復型部分構造法ではなく，CG 法をベースにした CGCG 法 (Coarse Grid Conjugate Gradient) という解法を用いている([12])。次節では CGCG 法の概要を述べ，3 節では解析モデルと結果，4 節では T2K 用に実装したハイブリッド並列化とそのパフォーマンスについて述べる。
解析対象はソリッド要素に分割した 31 層高層ビルである。兵庫県南部地震で記録された観測地震波でビル基礎を加震し，10 秒間の非線形動解析を行う。当初は，モデル作りおよび解析に手探りが続き，標準的なマルチコア型の T2K の性能も活かしきれていない点でも計算効率は良くなかった。しかし，さまざまな工夫を加えた結果，プロジェクト開始当初の 2008 年に T2K で出していたパフォーマンスは，2010 年 5 月には当初の約 4 倍となった。
1 ステップに 4 時間半，500 ステップの計算に，4 倍のパフォーマンス向上は，計算上，3か月強かかる計算が，1 か月弱でできるようになったことになる。しかし，ノード占有の問題，使用できるディスクの容量の問題，解析結果のダウンロードなど他にボトルネックはあり，理想通りには行かない。

2. CGCG 法
2.1 MPC と許容変位場の空間解くべき問題には十分な拘束条件 (ディリクレ条件) が課せられており，有限要素法に基づいて離散化された方程式 (下に述べる方程式 (2)) に現れる剛性行列は正則であるものとする。
その上で，一般に問題には MPC (Multi-Point Constraint) 条件が課せられる。構造解析
において，MPC は部品と部品をつなぐための必須の機能である。たとえば，いくつかの接合された部品に属する自由度が同時に動いたり，指定された線形関係のもとで動いたりする。したがって，自由度全体は退化している。より詳しくは，横長の行列 B によって自由度間の関係が，方程式Bu + r = 0(1)によって拘束されている。この解空間に，もとの問題の拘束条件を加えたものが 2.2 節でのべる述べる許容変位場の空間 V である。
2.2 CGCG 法
CGCG 法は，要約するとおよそ以下のような手法である。構造解析でも非線形問題，動的問題，あるいは熱伝導解析などは，静弾性構造解析の微分方程式 (あるいは楕円型方程式)を有限要素法で離散化した剛性方程式

Ku = f(2)に帰する。 K は剛性行列で正定値対称， f は外力， u は変位で，許容変位場の空間 V の元である。方程式 (2) はシンプルでしかも K はスパースなのだが，手ごわいのは規模が大きいことである。
解析対象領域 Ω をオーバラップしない部分領域に分割する。まず，並列プロセスにわたる「親」分割を行い，その上で各部分領域をさらに細かい「子」の部分領域に分割する。
この分割を Ω1 , L , Ω m とする。並列プロセス数は，プロセッサコア (たとえばフラット MPIの場合は使用するノードの全プロセッサコアとなる) にわたるように取るので，親分割はたとえば 8, 32, 128 あるいは本稿での 384 など，小さな領域数となる。子分割の分割数はパラメータだが，解析対象の形状，要素数などにより，数万から数百万になることもある。

細かい領域分割はコースグリッドの取り方からの要請であり，親分割に基づくコース運動は子分割に基づくコース運動で表現できるため，理論上は区別の必要はない。領域分割はオーバラップしないが，領域関境界はいくつかの部分領域に共有される。この領域関境界をグローバルに見て，内部境界とよぶ。
許容変位場の空間 V を，分割領域の，粗いグローバルな運動 (たとえば剛体運動) からなるコース空間 W1 とその K に関する共役直交空間 V1 に分解する。

V = W1 ⊕ K V1(3)W1 と V1 の上の K の働きをそれぞれ K% (1) ， K (1) と書くと，これらは V から W1 への射影 P とその相補的な射影 1 − P で表現できる。

K% (1) = PKPTK (1) = (1 − P) K (1 − P)T

(4)
(5)

もともとの方程式 (2) はつぎの 2 つの方程式

K% (1) w = f%(1)(6)

K u = f(1)(7)

f%(1) = Pf

(8)

f (1) = (1 − P ) f(9)(1)に分解される。ここでである。
一般に， W1 は粗い運動なので自由度は小さい。たとえば粗い運動として剛体運動をとれ
ば，1 分割領域の自由度は並進，回転を合わせて 6 自由度からなり，領域数×6 が W1 の自由度となる。これに比べれば V1 の自由度は領域数×6 を差し引いたものでしかなく， V の自由度と同じくらい大きい。そこで，(6) は直接法で解き，(7) は CG 法で解くことにする。(6)および (7) の解をそれぞれ w , u とすれば， w + u がもとの方程式 (2) の解となる。
この議論を拡張する。 V1 の中で， W1 とは異なる (独立な) 粗い運動をとり，上と同じことを繰り返すことを考える。もとの空間 V (正しくは方程式 (2)) に CG 法を直接には適用しないのと同様に，以下のように拡張する場合には V1 (正しくは方程式 (7)) に CG 法を直接には適用しない。

Ω1 , L , Ω m のそれぞれをさらに細かく領域分割し，そのグローバルな運動 (たとえば剛体運動) からコース空間 W2 をつくる。W2 の運動からは最初の領域分割の内部境界は除いておくか，ディリクレ条件により完全固定する。これにより，W2 は W1 とは独立な運動からなり，したがって， W2 は V1 の部分空間である。そこで， V1 をコース空間 W2 とその共役直交空間に分解する。

V1 = W2 ⊕ K (1) V2
すると， V はつぎのように書ける。

(10)

V = W1 ⊕ K W2 ⊕ K (1) V2(11)

V1 上の方程式 (7) は， V1 から W2 への射影により (4)，(5) と同様に得られる剛性行列をもつ 2 つの方程式

K% (2) w = f%(2)(12)

K u = f (2)(13)(2)に分解される。W2 は W1 と同じく小さな空間なので (12) は直接法で解く。方程式 (6) の解は直接法ですでに得られている。(13) は反復法で解く。これら 3 つの方程式の解をそれぞれ w1 , w2 , u としたとき， w1 + w2 + u がもとの方程式 (2) の解となる。
この操作は必要なだけ続けることができる。粗い運動の空間 W1 , W2,W3 , L に直接法を適用し，最後に残った共役空間に反復法 (CG 法) を適用する。解はコース空間の解の全部と共役空間の解の和で与えられる。
なお，実際には，本稿で述べる解析においては，ソルバの改良が以上述べたようにはコース空間の階層の一般化までは進んでおらず，使用した「開発版」(ここでは，リリース版と対比して開発版と呼ぶ) において，階層の深さは 2 と固定したものと同等の実装が行われていた。
2.3 非線形解析
非線形解析 (材料非線形解析，幾何学的非線形解析) は，非線型方程式の解法 (たとえばニュートン法) の中で得られる，接線剛性 K と外力増分 ∆f に対する変位増分 ∆u に関する方程式

K ∆ u = ∆f(14)に CGCG 法を用いる。方程式 (14) は方程式 (2) と同じ形の方程式である。
2.4 動解析
動解析では，時刻 t における各自由度の変位 u (t ) ，速度 v(t ) ，加速度 a (t ) を既知として，時刻 t + ∆t における u (t + ∆t ) , v (t + ∆t ) , a (t + ∆t ) を求める。ニューマークのβ法を改良した α 法と呼ばれる方法による。3 つの方程式v(t + ∆t ) = v (t ) + ∆t ((1 − γ ) a (t ) + γ a (t + ∆t ))

(15)

 1


u (t + ∆t ) = u (t ) + ∆tv(t ) + ∆t 2   − β  a (t )  + β a (t + ∆t )
(16)
 2

(17)
Ma (t + ∆t ) + Cv (t + ∆t ) + (1 + α )( I (t + ∆t ) − f (t + ∆t )) − α ( I (t ) − f (t )) = 0を連立する。(15)，(16) において，α は減衰を表わすパラメータであり， β , γ は時間方向の離散化パラメータ，(17) において M は質量行列，C は減衰行列である。ここではレーリー減衰を仮定する。(17) において I (t ) は内力であり，線形の場合は Ku (t ) に等しい。(15)，(16)，(17) から v (t + ∆t ) ，a (t + ∆t ) を消去すると u (t + ∆t ) の自由度分の未知数に対する非線形連立方程式が得られる。この非線形方程式にニュートン法を適用する。ニュートン法の中に現れる，増分修正量が満たす方程式が方程式 (14) と同じ形となり，ここに CGCG法が用いられる。

3. 解析モデルと解析結果仮想的に設計した 31 層からなる超高層ビルを，T2K 上の E-Simulator で解析する。高さは 129.7m，平面形状は 50.4m×36.0m の長方形である。架構形式は，センターコア形式の鋼構造であり，コアには制震ブレースを配置がされている。スラブはエレベータ部分を除いて全面に配置されている。図 1 に 3 次元 CAD モデルを示す。

図 1 31 層ビル CAD モデル
メッシュは 6 面体 1 次ソリッド要素で分割され，コンポーネントの結合に剛体梁，制震ブレースにはトラス要素が用いられている。6 面体要素数は 15,592,786，節点数は24,765,275，自由度は 74,295,825 である。剛体梁は基礎部の加震のためおよびブレースとガセットの接続のために合計 78,686 本，トラス要素 (ブレース) は 372 本使用されている。
メッシュの全体像は細かすぎて表示できない。図 2 に一部をクローズアップしたものを示す。ガセット

制震ブレース(トラス要素)

図 2 メッシュモデルの一部
まず固有値解析を行う。求めた固有周波数は低い方から 6 個，0.3074, 0.3484, 0.3822,0.9688, 1.052, 1.176Hz である。3 次モードまでの固有モードを図 3 に示す。

1 次モード2 次モード3 次モード

図 3 固有モード
つぎに動解析を行う。レーリー減衰行列は，1 次固有値に対して 2 %の初期剛性比例とする。加震は，JR 鷹取駅で記録された観測地震波のうちの，最大加速度を含む 10 秒間の，東西，南北，上下方向の加速度をそれぞれ長辺方向，短辺方向，上下方向の加速度履歴として柱脚下端部にそれぞれ空間方向に一様に入力して行う。
使用した地震波を図 4 に示す。
計算時間幅 ∆t の最大値は 0.01 秒に取られており，ニュートン法の収束状況などで自動的に分割される設定になっているが，本解析では 10 秒間を通じて分割されることなく，10 秒間の解析が 1000 ステップで終わっている。図 5 に，時刻 10 秒における相当応力分布に，20倍の変形を重ねた図を示す。図 6 に時刻 5 秒における 22 層付近の相当応力分布に，20 倍の変形を重ねた図を示す。図 6 で，水平方向に張られた梁と垂直方向の柱の接合部に赤い高応力部分 (430MPa 超) が生じており，降伏していることがわかる。

Acceleration [cm/s2]

800
600NSmax: 605.5

4002000-200-400-600012345Time [s]678910800max: 657

Acceleration [cm/s2]

600EW

400
2000-200-400-600-8000123Acceleration [cm/s2 ]300

45Time [s]6789max: 279.310UD

200
1000-100-200-300012345Time [s]6

図 4 入力地震波

78910図 5 時刻 10 秒における変形と相当応力分布 (変形は 20 倍)

図 6 時刻 5 秒における 22 層付近の変形と相当応力 (変形は 20 倍)

4. T2K における E-Simulator のチューニング
当初は，モデル作りおよび解析に手探りが続き，計算効率も良くなかった。標準的なマルチコア型の T2K の性能も活かしきれていなかった。さまざまな工夫を加えた結果，プロジェクト開始当初の 2008 年に T2K で出していたパフォーマンスは，2010 年 5 月には当初の約 4 倍となった。本節ではこの経緯について述べる。
図 7 に T2K (日立 HA8000) のノード構成を示す([13])。与えられた資源は 24 ノード，96ソケット，384 コアである。

図 7 HA8000 のノード構成
表 1 に E-Simulator のパフォーマンスの向上の経緯をまとめる。おおまかに 4 つのステップに分けられ，これらを便宜上，チューニングのステージ A，B，C，D と呼ぶ。解析ス
テップは自重による静解析から出発し，動解析が引き続いて行われる。静解析にも時刻が設定されており，動解析の出発は 1.0 秒からである。この表のとおり，調べたものは動解析の数ステップに過ぎず，実際にはある程度の長時間の実行をもとに評価すべきだが，時間的制約などから本表をもとにする。
表 1 E-Simulator パフォーマンスの変遷
チューニングのステージ
ソルバ
並列化
コンパイラ
解析モデル
実行日
時刻[s]
ステップ
1
1.01
2
1.02
3
1.03
動解析2ステップを基準とする速度向上率

ABCD
開発4.0αソルバ
開発4.4αソルバ(旧)開発4.4αソルバ(2008/5/15)(2010/3/11)(2010/5/26)
MPI 96並列
MPI 192並列
MPI 384並列
OpenMP 4並列(計384並列)
OpenMPなし
OpenMPなし
gcc 4.43
日立製
日立製
MPCによる層間の結合
共有節点による層間結合(非共有節点モデル)
2008/10/3
2010/3/23
2010/3/25
2010/5/27
計算時間[s]16207.4
9066.15667.03742.916646.7
8236.15698.44233.718217.1
7892.53635.8
未計測
1.0

2.0

2.9

3.9

チューニングのステージ A. 2008 年春より解析モデルの準備に入った。当初の解析モデルは，時間の制約などもあり，別々に作られた各層が MPC で結合されるという，明らかに能率の悪い解析モデルから出発した。このモデルを用いて当時のソルバ開発版 4.0α (リリース版に対比して開発版と呼んでいる) で解析を開始した。最初は 24 ノード 384 コアにフラットに MPI を適用したが，ソケット間の通信が飽和し，論外のパフォーマンスだったため，各ソケットから 2 コアのみを取り，合計 192 コアにフラットに MPI を適用した。これがチューニングステージ A である。
チューニングのステージ B. 2009 年 10 月には，開発版 4.4αにおいて，OpenMP による並列化を行った。文献 [14] などを参考にした。このチューニングでは，2010 年 3 月時点で最新の gcc 4.4.3 を用いている (本稿執筆時点の最新版は 4.6)。gcc の OpenMP 対応は 4.2から正式サポートとなり，gcc のバージョンアップによる OpenMP 効率の向上は大きい。
MPI 並列は 1 ソケット 1 プロセスとし，1 プロセスに 4 スレッドを生成した。スレッドの割り当てを図 8 に示す。

図 8 チューニングのステージ B のスレッドの割り当て
この並列化は，主として粒度の粗い，部分領域についてのループに渡り以下のように行った。
① 並列化関係部のスレッドセーフ化
多くの並列化関係関数をスレッドセーフ化した。また，もともと使用している LAPACKをスレッドセーフ版と置き換えた。これは以下の②，③，④の準備である。
② CG 法前処理ループの一部 (部分領域ループ)
並列同時処理を禁止する (Thread-unsafe) バリア数が多く，必ずしも大きなパフォーマンス向上にはつながらなかった。

③ CG 法行列ベクトル積，応力積分ループの一部 (部分領域ループ)
最も大きな効果が得られた部分で，1 ソケット 4 スレッドを用い，ほぼ 4 倍のパフォーマンスを確認できた。
④ 各種ベクトル行列処理や最大値探索などの関数 (非部分領域ループ)大きなパフォーマンス向上にはつながらなかった。
以上に加え，チューニングステージ B (ソルバ開発版 4.4α) では CGCG 法ソルバにおけるコースグリッドにおける直接法の解法に MPI による並列化を適用し，また CG 法トレランス，ニュートン法トレランスなどの見直しを行った。
チューニングのステージ C. 前述のようにチューニングステージ A，B では，解析モデルは各層を別々に作り，それらを MPC で接続したものだが，計算上は能率が悪い。そこで，層間の MPC 接続部分を，MPC を使わずに共有節点を用いて作りなおしたのがチューニング
ステージ C である。ただし，3 節で述べたように，剛体梁が 78,686 本使われている。剛体梁要素は MPC を用いて作られており，解析モデルにはモデル構築上必須の多くの MPC が残っている。
チューニングのステージ D. 2010 年 5 月，ソルバ開発版 4.4αで MPC 処理における全体通信を減らすことにより通信量の削減を行った。表 1 に示す通り MPI プロセスが 1 ソケットの全コアに適用されている。この結果，OpenMP の並列化よりフラット MPI の並列化の方が勝る結果が得られた。
パフォーマンスの比較は動解析の第 2 時間ステップに対して行った。チューニングステ
ージ A を基準とするとステージ B，C，D でそれぞれ約 2 倍，3 倍，4 倍の速度となった。
これは，モデル A では 1 ステップに 4 時間半，500 ステップの計算に，ステージ A では計算上，3 か月強かかった計算が，ステージ D では 1 か月弱でできるようになったことになる。ただし，最初に述べたように，ノード占有の問題，ディスクの容量の制限，解析結果のダウンロードなど他にボトルネックはあり，これは理想であって，現実とは異なる。
パフォーマンス向上の効果を項目別に調べると表 2 のような結果になった。この表中，1項目目のニュートン法反復の減少は，ステージ B の項で述べたトレランスの見直しの効果である。
表 2 各項目別パフォーマンス向上の効果
項目

計算時間短縮1ニュートン法反復の減少

約 40%

2

OpenMP 並列化

約 20%

3

直接法並列化

約 5%

4

解析モデルの共有節点化

約 30%

5

MPC 通信量削減

未測定

5. まとめ
E-Simulator の概要と T2K におけるチューニングについて述べた。24 ノード，96 ソケ
ット，384 コア，384 並列プロセスでは，フラット MPI がハイブリッド並列化を上回る結果となった。しかしこれは，チューニングの中間結果である。現状のままでも，並列プロセス数をより多く取ると逆の結果が得られることはあり得る。また，OpenMP あるいは場合により Pthread チューニングを強化する余地はある([15])。一般に，コンパイラもスレッド並列化対応を強めている。
E-Simulator は今後さらにパフォーマンスを上げていくのと同時に，解析モデルの作成の自動化なども検討されている。

謝辞
本稿で述べた解析では，防災科学技術研究所・数値震動台開発分科会からの全般にわたる指導を得た。E-Simulator および ADVC の開発と本稿に述べた解析は，株式会社アライ
ドエンジニアリングのメンバーが行った。ここに謝意を表します。

参考文献,
[1]

兵庫耐震工学センターウェブサイト

http://www.bosai.go.jp/hyogo/[2]

アライドエンジニアリングウェブサイト
http://www.alde.co.jp/[3]

ADVENTUREプロジェクトウェブサイト
http://adventure.sys.t.u-tokyo.ac.jp/jp/[4]

Makoto Ohsaki, Tomoshi Miyamura, Masayuki Kohiyama, Muneo Hori, Hirohisa
Noguchi, Hiroshi Akiba, Koichi Kajiwara, Tatsuhiko Ine, High-precision finite
element analysis of elastoplastic dynamic responses of super-highrise steel frames,
Earthquake Engineering Structural Dynamics, Vol. 38, pp. 635-654, 2009[5]

大崎純，宮村倫司，小檜山雅之，磯部大吾郎，秋葉博，堀宗朗，梶原浩一，井根達比古，鋼構造建築骨組の弾塑性地震応答の高精度有限要素解析，第14回計算工学講演論
文集，2009年5月12日～14日，東京大学生産技術研究所[6]

Tomoshi MIYAMURA, Makoto OHSAKI, Masayuki KOHIYAMA, Daigoro ISOBE,
Kunizo ONDA, Hiroshi AKIBA, Muneo HORI, Koich KAJIWARA and TatsuhikoINE, Large-Scale FE-Analysis of Steel Building Frames Using E-Simulator,
Proceedings of SNA＋MC2010, October 17-21, 2010, Tokyo[7]

宮村倫司，大崎純，小檜山雅之，張景耀，磯部大吾郎，秋葉博，堀宗朗，山下拓三，
梶原浩一，崩壊解析のための鋼構造骨組の高精度有限要素モデル，第16回計算工学講演会論文集，2011年5月25日～5月27日，東京大学柏キャンパス，発表予定

[8]宮村倫司，大崎純，堀宗朗，山下拓三，磯部大吾郎，小檜山雅之，秋葉博，梶原浩一，
E-Simulator による４層鋼構造骨組の仮想震動実験と骨組解析の比較，日本建築学会
2011 年度大会 (関東) 学術講演梗概集，2011 年 8 月 23～25 日，早稲田大学，発表予
定[9]

山下拓三，宮村倫司，小檜山雅之，大崎純，堀宗朗，秋葉博，梶原浩一，E-Simulatorによる合成梁の高精度コンポーネント解析，日本建築学会 2011 年度大会 (関東) 学術講演梗概集，2011 年 8 月 23～25 日，早稲田大学，発表予定[10] T. Yamashita, T. Miyamura, M. Ohsaki, M. Kohiyama, D. Isobe, K. Onda, M. Hori,J.Y. Zhang, H. Akiba and K. Kajiwara, High-Precision FE-Analysis for Seismic
Collapse Simulation of Steel Building Frames, Proceedings of COMPDYN 2011,May 25-28, 2011 , Corfu, Greece, 発表予定
[11] Makoto Ohsaki, Tomoshi Miyamura, Masayuki Kohiyama, Jingyao Zhang,
Daigoro Isobe, Kunizo Onda, Takuzo Yamashita, Muneo Hori, Hroshi Akiba andKoichi Kajiwara, High-Precision FE-Analysis for Collapse Simulation of SteelFrames Considering Composite Beam Effect, Proceedings of EUROSTEEL 2011,
August 31 - September 2, 2011, Budapest, Hungary, 発表予定
[12] 秋葉博，大山知信，CGCG 法とその拡張アルゴリズム，2010 年度 第 23 回 日

本機械学会計算力学講演会講演論文集，2010 年 9 月 23 日～25 日，北見工業大学」概要[13] 「T2K オープンスパコン（東大）
http://nkl.cc.u-tokyo.ac.jp/10w/CW-intro/HA8000-overview.pdf[14] 中島研吾，T2K オープンスパコン (東大) 連載講座番外編：Hybrid プログラミング
モデルの評価http:// www.cc.u-tokyo.ac.jp/publication/news/VOL11/No4/200907tuning.pdf[15] 荻野正雄，MPI/Pthread ハイブリッド並列化による領域分割型有限要素法の高速化，先駆的計算科学に関するフォーラム 2010 資料，2010 年 9 月 29 日，九州大学伊都キャンパス
http:// www.cc.kyushu-u.ac.jp/scp/users/forum2010-09/6.pdf