HΦ におけるいくつかの工夫
三澤貴宏

東京大学物性研究所計算物質科学研究センター
はじめにI.

厳密対角化、とくに Lanczos 法を用いた大規模計算を行なう上での計算上の幾つかの工夫を解説する。行列の大きさが、行列要素が全て確保できないほど巨大になってくると、行列要素を逐次生成しながら行列-ベクトル積の計算を行なうことになるが、その計算量は莫大になるため、計算コストを削減することが、現実的な時間で計算を終えるために必要であり、様々な工夫を行なうことになる。厳密対角化の先駆的なパッケージである TITPACK ver.21 およびそのマニュアルでは計算上の工夫をどう行なうかを明確に解説してあり、その後の対角化のコードを作成する人にとって大きな参考となってきた。この解説では、TITPACK で用いている工夫に付け加えて、HΦ で用いている工夫について解説する。

記法などについてII.

この説明では、spin 1/2 の系を例にとって説明を行なう。Hilbert 空間は, S z に関して対角的な表示をとることで、Hilbert 空間の全ての要素を単一の実空間配置で表すことができる。例えば、| ↑, ↑, ↑, ↓⟩ などは ↑=1, ↓ =0 とすれば, 2 進数で [1110] と表すことができる。状態を２進数表示すな
わち bit で表現することで、交換相互作用などは bit の入れ替えとして bit 演算を用いることで計算を効率的に行なうことができる。また、Hubbard 模型でも同様に表現が行える。Hubbard 模型の場合は up,down の電子をそれぞれ 0,1 で表現することで、2 進数で表せる。例えば, | ↑↓, ↑, ↓, 0⟩ は[11,10,01,00] と表せる。

III.
A.

HILBERT 空間の制限について
2 次元サーチ法

2 次元サーチ法は、TITPACK で用いられている基本的な方法であり、Hilbert 空間を制限した場合にその逆引きを効率よく行なう方法である。まずは、この方法を解説する。この解説は TITPACK

1

ver.2 のマニュアルを参考にしている。例えば、全 S z が一定の空間で考える場合は、全実空間配置の中から、指定した全 S z を持つものを抜き出して、その一つ一つの数 (=実空間配置)list に格納して行けばよい。問題はある実空間配置が制限した Hilbert 空間の中で、何番目の要素なのかを逆引きする場合である。一番、愚直な実装は全ての実空間要素に対して逆引きのリストもっておくことだがこれはあまりにもコストが大きすぎる。この逆引きを効率よくおこなうのが 2 次元サーチ法である 1,2 。

jb

ja

ja+jb

000 111011001011001 011

112001101

001 101

123001110

001 110134010011

010 011415010101

010 101

426010110010 110

437iib000111ia

図 1: 2 次元サーチ法の図。i の bit を２つに分割して、分割した bit の情報から、何番目の要素かを再現する方法。

この方法の肝は２進表示した状態の bit を２つに分割することにある (これが 2 次元サーチ法の名前の由来)。図 1 に 6 site, 全 S z = 0 の場合を示している。この場合、i を２つに分割した、上位bit が ib, 下位 bit が ia である。jb, ja は以下のルールで逐次的に決めていく。
1. 最初は jb=0, ja=1
2. ib が変化しないなら、jb はそのままで ja を１つずつ増やしていく
3. ib が変化したときは、ja を１にして、jb はひとつ前の ja+jb の値にする
こうしておけば、ja+jb が i が何番目かの要素を示すようになる。それぞれの ib, ia に対して、jb,ja は一意に決めるので、list jb[ib], list ja[ia] という２つの配列 (配列の大きさは 2N/2 (N は全サイト数) なので小さい) をもっておけば、i から ja+jb を生成することができるようになる。

22 次元サーチ法のアルゴリズムの擬似コードを Algorithm 1 に示す。例として N site の spin-1/2の系を想定している。また, 全 S z =Nup としている (つまり, 1 の bit の個数が Nup)。list 1 は要件をみたす (つまり所定の S z を満たす) Hilbert 空間の元を格納する配列であり、get ib,get ia はそれぞれ i から上位 bit ib, 下位 bit ia をとり出す関数である。

Algorithm 1 Two-dimensional searchjb← 0, ja← 0, icnt← 0, ib old ← 0
for i=0 to 2N -1 do
if number of bit in i = Nup thenlist 1[icnt] ← iib ← get ib(i)//get ib from iia ← get ia(i)//get ia from i

if ib = ib old thenja ← ja+1
else
ib old ← ibja ← 1jb ← icnt-1
end iflist ja[ia] ← jalist jb[ib] ← jb
icnt ← icnt+1
end ifend forB.

2 次元サーチ法のスレッド並列化について

さて、次は 2 次元サーチ法のアルゴリズムの並列化について述べる。2 次元サーチ法には逐次的な操作が入ってしまっているので、通常は並列化が不可能なように見えるが、ループを ib, ia に関するループに分割することで、並列計算が可能である。まず, 2 次元サーチののアルゴリズムからib, ia のループを入れ子にすることが可能である。ib で外側のループを回し、ia で内側のループを回せばよい。そして, ib に依る部分が jb だけであり ib と jb は一対一に対応するので、ib に対応す3

Algorithm 2 Parallelization for two-dimensional search algorithmjb← 0for ib=0 to 2N/2 -1 dolist jb[ib] ← ibib Nup ← count bit(ib)jb ← jb+Binomial(N/2,Nup-ib Nup)end forfor ib=0 to 2N/2 -1 do
jb ← list jb[ib]ib Nup ← count bit(ib)ja ← 1for ia=0 to 2N/2 -1 doia Nup = count bit(ia)if ib Nup+ia Nup=Nup thenlist 1[ja+jb] ← ia+ib*ihfbit
list ja[ia] ← jalist jb[ib] ← jb
js ← ja+1
end ifend forend forる jb がわかれば、ib に関するループは独立に計算でき、並列計算が行える。ib に対応する jb は ibに対応する ia の数を二項係数を計算することで、前もって知ることができる。この ib に関するループの部分は並列化不可能だがループ長が 2N/2 なので、この部分の計算時間は実質無視できる。そのあとで、ib,ia の入れ子のループで外側の ib のループを openmp で thread 並列化すればよい。
擬似コードを Algorithme 2 に示す。ib, ia の 1 の bit の個数を数える関数 (count bit) と二項係数を計算する関数 (Binomial) が必要だが、それ以外はもとの 2 次元サーチ法のアルゴリズムと同じである。また、ihfbit=2N/2 で ia,ib からもとの i を復元するために必要な数である。

4C.

同じ 1 の bit の個数をもつ次の bit を生成する方法

2 次元サーチ法の並列化のアルゴリズムをみると、並列化はされるが全体のループ長は依然として 2N で計算コストのオーダーは変わっていないことがわかる。実はこの計算コストを縮めることも可能である。鍵となるアルゴリズムは「ある数と同じ 1 の bit の個数を持つ数のなかで次に大きい数を求める方法」である。これは「ハッカーの楽しみ」3,4 という本に乗っている方法である 1 このアルゴリズムの擬似コードを Algorithm 3 にのせる。注意しなければならないのは x が 0 より大きくある必要があるということである。わずか、6 行で次の数が生成できる簡潔なアルゴリズムである。

Algorithm 3 snoob: Get next larger bitsmallest ← x & (-x)
ripple ← x + smallest
ones ← x ˆ ripple
ones ← (ones >> 2)/smallest
next x ← ripple | onesD.

さらなるスレッド並列化について「ハッカーの楽しみ」の方法を用いて、さらなる効率化をしてみよう。まず、最初の 2 次元サーチ法のアルゴリズムそのものにこの方法をつかってもよいがそれだと並列化はできない。2 次元サーチ法の並列化バージョンに対して、ia のループの部分で「ハッカーの楽しみ」の方法を使うのことで並列化ができる。

Algorithm4 に擬似コードを示す。基本は 2 次元サーチ法の並列化版と同じで ia のループに関してのみ、snoob を用いている。snoob は 0 を引数にとれないことに注意してほしい。これによって、計算コストは減少して、おちて N=36 といった巨大な系でもヒルベルト空間作成のコストはほぼ 0になる。また、飽和磁化近傍のようにシウテムサイズ N 自体は大きいが、状態数が少ないところの計算も行えるようになっている。例えば, 48 サイト 2up,46down といった条件のヒルベルト空間の作成はすぐに行える。HΦ では、このアルゴリズムを用いた計算は modpara.def で CalcHS=1 とすることで実行できる。

1

この本の作者はこのアルゴリズムの紹介の前に、「いったい誰がそんなものを計算したがるのかという疑問が生じる」と述べているが、実際に役に立つ例がここにあった!

5

Algorithm 4 Parallelization for 2D search algorithm IIjb← 0for ib=0 to 2N/2 -1 dolist jb[ib] ← ib
ib Nup ← count bit(ib)jb ← jb+Binomial(N/2,Nup-ib Nup)end forfor ib=0 to 2N/2 -1 do
jb ← list jb[ib]ib Nup ← count bit(ib)ja ← 1
if ib Nup ≤ Nup thenia = 2Nup−ib Nup -1
if ia < ihfbit thenlist 1[ja+jb] ← ia+ib*ihfbit
list ja[ia] ← jalist jb[ib] ← jb
ja ← ja+1
if ia ̸= 0 thenia ← snoob(ia)while ia < ihfbit dolist 1[ja+jb] ← ia+ib*ihfbit
list ja[ia] ← jalist jb[ib] ← jb
ja ← ja+1ia ← snoob(ia)end while
end ifend ifend ifend for6

IV.

FERMION サインの数え方について

ハバード模型など遍歴電子系のフェルミオンサインについて説明する。つねに粒子が存在するスピン系とは違い、遍歴電子系の場合は電子がいないサイトがあるので、ホッピング項などを作用させるときにフェルミオンの交換関係から符号が生じる。例えば、以下のような場合である。

c†2↑ c0↑ |0, 0, ↓, ↑⟩ = |0, ↑, ↓, 0⟩ この場合、ホッピングする間の 1 の bit の偶奇を計算する必要がある。
ここでも「ハッカーの楽しみ」3,4 にのっている方法が役にたつ。与えらた数の 1 の bit のパリティをを求めるアルゴリズムを用いて、元の bit(orbbit) から ±1 の sgn を求めている。もとの orgbit は64bit の整数を仮定しているに注意してほしい。
Algorithm 5 parity counting
bit ← orgbit ˆ (orgbit>>1)
bit ← bitˆ (bit>>2)
bit ← bitˆ (bit>>4)
bit ← bitˆ (bit>>8)
bit ← bitˆ (bit>>16)
bit ← bitˆ (bit>>32)
sgn ← 1-2*(bit & 1); //sgn=±11

http://www.stat.phys.titech.ac.jp/ nishimori/titpack2 new/index-e.html.

2H. Q. Lin, Phys. Rev. B 42, 6561 (1990).

3

S. H. Warren, Hacker’s Delight (2nd Edition) (Addison-Wesley, 2012), ISBN 0321842685.

4

S. H. Warren, ハッカーのたのしみ―本物のプログラマはいかにして問題を解くか (エスアイビーアクセス, 2004), ISBN 4434046683.

7