陳富州 程晨 2) 羅洪剛 2)?

1)(蘭州大學(xué)物理科學(xué)與技術(shù)學(xué)院,蘭州 730000)

2)(北京計(jì)算科學(xué)研究中心,北京 100084)

1 引 言

密度矩陣重正化群方法(DMRG)[1,2]是研究相互作用量子系統(tǒng)最重要的多體數(shù)值方法之一.眾所周知,多體量子系統(tǒng)的復(fù)雜度隨著系統(tǒng)尺寸指數(shù)增長(zhǎng),給數(shù)值計(jì)算帶來(lái)了極大的困難.而DMRG給出了一種非常有效的截?cái)嘞柌乜臻g的方法,保留有限個(gè)狀態(tài)并通過(guò)掃描即可得到變分收斂的基態(tài)或低激發(fā)性質(zhì).應(yīng)用于自旋1的海森伯鏈時(shí),該方法僅保留數(shù)百狀態(tài)并通過(guò)數(shù)次掃描即可得到相對(duì)誤差 1 0-9左右的基態(tài)能量,而所需的計(jì)算量?jī)H隨格點(diǎn)尺寸線性增長(zhǎng)[3,4].作為求解一維格點(diǎn)模型基態(tài)的成熟方法,DMRG也不斷被應(yīng)用于其他各類(lèi)問(wèn)題并取得了一定的成功,如動(dòng)量空間中的哈密頓量[5]與量子化學(xué)問(wèn)題[6-8]、量子系統(tǒng)的時(shí)間演化問(wèn)題[9-11]、準(zhǔn)二維以及二維量子格點(diǎn)模型[12-14]等.這些嘗試極大擴(kuò)展了DMRG的應(yīng)用,但同時(shí)也對(duì)該數(shù)值算法的優(yōu)化提出了更高的要求.以二維相互作用電子系統(tǒng)為例,該系統(tǒng)包含非常豐富的物理,例如高溫超導(dǎo)條紋相[15-17]、量子自旋液體[18,19]等.然而,在面對(duì)二維或者準(zhǔn)二維格點(diǎn)模型時(shí),DMRG所需的保留狀態(tài)數(shù)大致隨格點(diǎn)寬度指數(shù)增長(zhǎng),得到收斂結(jié)果所需的掃描次數(shù)也大大增加.此時(shí),DMRG所需的計(jì)算量和存儲(chǔ)量較大,使用各種方法優(yōu)化算法顯得十分必要.

在實(shí)際應(yīng)用中,對(duì)DMRG算法的優(yōu)化與加速主要體現(xiàn)在兩個(gè)方面.一方面,人們不斷改進(jìn)DMRG算法本身以減少計(jì)算量,包括使用各種好量子數(shù)對(duì)角化小的希爾伯特子空間[20,21],使用動(dòng)態(tài)保留狀態(tài)數(shù)節(jié)約計(jì)算資源[22,23],使用好的初始預(yù)測(cè)波函數(shù)減少對(duì)角化方法迭代次數(shù)[24],使用單格點(diǎn)DMRG方法減少計(jì)算量與內(nèi)存[25,26]等.另一方面,人們結(jié)合DMRG算法的特性,發(fā)揮高性能計(jì)算機(jī)的并行計(jì)算能力進(jìn)一步縮短計(jì)算時(shí)間,包括實(shí)空間并行[27]、共享存儲(chǔ)的多核行[28]、分布式存儲(chǔ)的多節(jié)點(diǎn)行[29]以及CPU-GPU異構(gòu)并行[30]等.

相比CPU,GPU具有較高的浮點(diǎn)計(jì)算能力以及較大的存儲(chǔ)帶寬,因此在通用計(jì)算中得到了大量的應(yīng)用.當(dāng)前很多求解強(qiáng)關(guān)聯(lián)格點(diǎn)模型的數(shù)值方法(比如：精確對(duì)角化方法[31]、量子蒙特卡羅方法[32]以及張量乘積態(tài)方法[33]等)都實(shí)現(xiàn)了基于GPU的并行優(yōu)化.在CPU-GPU異構(gòu)并行環(huán)境中,異構(gòu)并行算法可以同時(shí)發(fā)揮CPU和GPU的計(jì)算性能,并且可以實(shí)現(xiàn)內(nèi)存和GPU顯存的分布式存儲(chǔ),在一定程度上減小了GPU顯存容量對(duì)計(jì)算規(guī)模的限制.最近,Nemes 等[30]提出了DMRG方法的一種CPU-GPU異構(gòu)并行實(shí)現(xiàn),并將其應(yīng)用到一維格點(diǎn)模型的基態(tài)能量計(jì)算.在保留狀態(tài)較多時(shí),最耗時(shí)的哈密頓量對(duì)角化部分在GPU中獲得了接近峰值的運(yùn)算性能.在對(duì)角化過(guò)程中,對(duì)占用存儲(chǔ)空間較多的Davidson方法實(shí)現(xiàn)了內(nèi)存和GPU顯存的分布式存儲(chǔ),減小了GPU顯存的容量限制.然而與此同時(shí),Davidson方法在計(jì)算中內(nèi)存和GPU顯存之間需要通信向量數(shù)據(jù),其性能會(huì)受到通信帶寬的限制.另一方面,該方案中哈密頓量和波函數(shù)基于兩子塊表示,其子塊算符需要的存儲(chǔ)量大致相當(dāng)于目前流行的四子塊表示的d2倍(d為單格點(diǎn)希爾伯特空間維數(shù)),在實(shí)際應(yīng)用中存在很大的局限性.

在前人工作的基礎(chǔ)上,本文提出了一種新的CPU-GPU異構(gòu)并行優(yōu)化算法,并使用四子塊表示的DMRG超塊與波函數(shù).為了說(shuō)明該方法的有效性,將其應(yīng)用到準(zhǔn)二維模型的求解,獲得了不同保留狀態(tài)數(shù)下的基態(tài)能量及其性能基準(zhǔn).在Hubbard梯子模型的例子中,得到了非均勻的電荷密度分布,與高溫超導(dǎo)銅氧面中觀測(cè)到的條紋相定性一致.我們希望CPU-GPU異構(gòu)并行算法能進(jìn)一步推進(jìn)DMRG在求解二維與準(zhǔn)二維模型、時(shí)間演化、量子化學(xué)等問(wèn)題中的應(yīng)用,同時(shí)引起強(qiáng)關(guān)聯(lián)領(lǐng)域?qū)PU算法的關(guān)注和重視.

本文的內(nèi)容包括以下幾部分：第2節(jié)回顧了有限D(zhuǎn)MRG算法,基于四個(gè)子塊表示的DMRG實(shí)現(xiàn),介紹了該工作采用的基準(zhǔn)模型,并得到了在CPU執(zhí)行時(shí)各個(gè)部分的計(jì)算時(shí)間占總時(shí)間的比例;第3節(jié)針對(duì)DMRG中最耗時(shí)的哈密頓量對(duì)角化部分給出了異構(gòu)并行優(yōu)化方法;第4節(jié)以計(jì)算4腿Hubbard模型基態(tài)為例,對(duì)比了異構(gòu)并行優(yōu)化方法與單個(gè)CPU中MKL并行的性能;最后給出本文的總結(jié).

2 有限D(zhuǎn)MRG方法和基準(zhǔn)模型

有限尺寸格點(diǎn)模型的DMRG計(jì)算分為兩部分：首先執(zhí)行無(wú)限D(zhuǎn)MRG方法,這可為后面的計(jì)算提供較好的初始波函數(shù);然后執(zhí)行有限D(zhuǎn)MRG掃描,其中每一步有限D(zhuǎn)MRG掃描會(huì)優(yōu)化系統(tǒng)塊的狀態(tài),直到收斂至基態(tài).在準(zhǔn)二維梯子模型的計(jì)算中,有限D(zhuǎn)MRG部分需要多次掃描,且保留很多的狀態(tài)才能收斂,幾乎占用整個(gè)DMRG計(jì)算的全部時(shí)間,因此本文主要考慮該部分的并行優(yōu)化.有限D(zhuǎn)MRG計(jì)算中向左掃描和向右掃描非常相似,本文以向右掃描為例介紹有限D(zhuǎn)MRG方法.每一步優(yōu)化的過(guò)程中整個(gè)格點(diǎn)系統(tǒng)由四個(gè)部分構(gòu)成 (圖 1)：S,s1,e1 和E,其中S和s1 構(gòu)成系統(tǒng)塊,E和e1 構(gòu)成環(huán)境塊,系統(tǒng)塊和環(huán)境塊構(gòu)成超塊.基于四個(gè)子塊,波函數(shù)有如下形式：

其中|s〉,|e〉,|σs1〉和|σe1〉分 別為子塊S,E,s1 和e1上的狀態(tài).同樣,基于四個(gè)子塊,一般形式的哈密頓量H可以表示為

其中OS,OE,Os1和Oe1分別為子塊S,E,s1 和e1上的算符.以求解系統(tǒng)的基態(tài)為例,給出基于四子塊表示的DMRG一步優(yōu)化過(guò)程,如算法1所示.重復(fù)執(zhí)行DMRG的一步優(yōu)化過(guò)程,進(jìn)行多次實(shí)空間掃描,即可得到收斂的結(jié)果.在算法1中,對(duì)角化超塊哈密頓量最為耗時(shí),通常人們采用稀疏矩陣對(duì)角化方法(Lanczos方法、Davidson方法等).此時(shí),不需要超塊哈密頓量的矩陣表示,僅需要算符在四個(gè)子塊中的矩陣表示,并迭代執(zhí)行|φ〉=H|ψ〉,文獻(xiàn)[4]中給出了該操作的高效算法,具體過(guò)程如算法2.在后面的描述中,我們稱(chēng)算法2中第2,5和8行對(duì)應(yīng)的循環(huán)分別為step1,step2和step3.可以看到該算法(st)ep1和step3中為矩陣乘法,對(duì)應(yīng)的計(jì)算量為OD3,其中D為DMRG保留狀態(tài)(數(shù).)step2中僅包含向量操作,對(duì)應(yīng)的計(jì)算量為OD2.在保留狀態(tài)較多時(shí),哈密頓量作用在波函數(shù)的計(jì)算量主要由step1和step3決定.而對(duì)角化哈密頓量的總時(shí)間線性依賴(lài)于對(duì)角化方法的迭代次數(shù),本文實(shí)現(xiàn)的有限D(zhuǎn)MRG中通過(guò)上一步優(yōu)化波函數(shù)給出一個(gè)較好的初始迭代波函數(shù),可以有效地加快對(duì)角化方法的收斂[24].

圖1 超塊中的四個(gè)子塊Fig.1.4 Sub-blocks of super-block.

?

?

為了獲得異構(gòu)并行優(yōu)化方法的性能基準(zhǔn),選取4腿梯子上的Hubbard模型為例子,并利用DMRG求解其基態(tài).具體地,系統(tǒng)哈密頓量為

這里t為兩個(gè)格點(diǎn)之間的電子躍遷能,U為單個(gè)格點(diǎn)上的電子在位庫(kù)侖排斥能,為格點(diǎn)i上自旋σ的費(fèi)米子產(chǎn)生(湮滅)算符,為自旋σ的粒子數(shù)算符,其中σ∈{↑,↓}.在梯子模型中,格點(diǎn)i包含x,y兩個(gè)方向的分量,以x標(biāo)記梯子長(zhǎng)邊坐標(biāo),y標(biāo)記梯子短邊坐標(biāo).在后面的測(cè)試計(jì)算中,選取t=1 為能量單位,沿著長(zhǎng)邊和短邊方向分別使用開(kāi)邊界和周期邊界條件,梯子長(zhǎng)度為16.對(duì)于其他參數(shù),取相互作用U=8,總電荷密度0.875,即 1 /8 空穴摻雜,這是在其他數(shù)值工作中觀測(cè)到條紋相的典型參數(shù)[13,34].另外本文的實(shí)現(xiàn)利用了該模型的總粒子數(shù)和總自旋在z方向投影兩個(gè)好量子數(shù),每個(gè)子塊的希爾伯特空間可以被劃分為多個(gè)子空間,每個(gè)子空間中的狀態(tài)對(duì)應(yīng)于相同的量子數(shù).此時(shí)模型(1)中算符的矩陣表示為分塊矩陣,相應(yīng)地,超塊的希爾伯特空間可以用其四個(gè)子塊的好量子數(shù)劃分為多個(gè)子空間.

模型(1)哈密頓量的矩陣表示為實(shí)對(duì)稱(chēng)矩陣,本文使用Davidson方法[35,36]對(duì)角化該哈密頓量.Davidson方法是一種使用預(yù)條件技術(shù)的子空間迭代方法,算法3給出了該方法每一步迭代的具體操作.可以看到,每一步迭代都需要作用哈密頓量到波函數(shù),并且包含多個(gè)向量操作.其中向量操作的計(jì)算量和存儲(chǔ)量均線性依賴(lài)于Davidson方法子空間中向量的個(gè)數(shù)和向量的維數(shù).在性能測(cè)試中,使用7次有限D(zhuǎn)MRG掃描收斂到基態(tài),Davidson方法子空間中向量個(gè)數(shù)最大為11.

?

本文所有計(jì)算使用的異構(gòu)并行環(huán)境包含1個(gè)CPU (Intel Xeon E5-2650 v3,10核,主頻2.3 GHz)和Nvidia Tesla K80卡中的1個(gè)GPU (含12 GB顯存).程序的實(shí)現(xiàn)基于CUDA環(huán)境,CPU和GPU中的矩陣向量操作分別調(diào)用了Intel MKL和CUBLAS中的子程序.為了描述矩陣操作的性能,估計(jì)了單位時(shí)間內(nèi)處理器執(zhí)行的浮點(diǎn)操作的次數(shù),其中總的浮點(diǎn)操作次數(shù)在保留狀態(tài)數(shù)較大時(shí)可以由所有矩陣乘法中浮點(diǎn)操作次數(shù)之和近似.為獲取異構(gòu)并行的優(yōu)化表現(xiàn),需要給出CPU上的并行性能作為基準(zhǔn).哈密頓量對(duì)角化中包含大量矩陣乘法,這里首先測(cè)試了隨機(jī)方陣乘法的性能,結(jié)果見(jiàn)圖2(a).可以看出,當(dāng)矩陣尺寸較大(大于400)時(shí),矩陣乘法的性能較高,并隨著矩陣尺寸增大而增大.對(duì)于DMRG算法,圖2(b)中的結(jié)果表明,保留狀態(tài)數(shù)越大并行計(jì)算性能越高,并逐漸接近峰值性能.另外也測(cè)試了對(duì)角化哈密頓量以及哈密頓量作用于波函數(shù)部分占總計(jì)算時(shí)間的比例,如圖3 所示.在我們所關(guān)心的保留狀態(tài)數(shù)范圍內(nèi),對(duì)角化哈密頓量的耗時(shí)比例超過(guò)總計(jì)算時(shí)間的90%,其中作用哈密頓量到波函數(shù)占總時(shí)間比例超過(guò)80%.因此,我們的工作主要針對(duì)該部分進(jìn)行異構(gòu)并行優(yōu)化.

圖2 CPU中作用哈密頓量在波函數(shù)上的性能 (a)矩陣乘法的浮點(diǎn)性能;(b)作用哈密頓量于波函數(shù)的性能,及矩陣乘法中的最大矩陣尺寸Fig.2.Performance of acting the Hamiltonian on the wave function in CPU：(a)The matrix multiplication performance;(b)the performance of acting the Hamiltonian on the wave function,and the maximum matrix size of the matrix multiplications.

圖3 對(duì)角化哈密頓量和作用哈密頓量到波函數(shù)操作占總計(jì)算時(shí)間的比例Fig.3.Time ratio of diagonalization of the Hamiltonian and acting the Hamiltonian on the wave function to the total time cost.

3 DMRG的異構(gòu)并行實(shí)現(xiàn)

主要考慮DMRG方法在準(zhǔn)二維模型中的應(yīng)用,并且針對(duì)有限D(zhuǎn)MRG中計(jì)算量較大的哈密頓量對(duì)角化部分進(jìn)行并行優(yōu)化.對(duì)于準(zhǔn)二維問(wèn)題,DMRG方法達(dá)到較高精度通常需要保留較多的狀態(tài),相應(yīng)地,對(duì)角化哈密頓量時(shí)會(huì)執(zhí)行一些大尺寸矩陣操作.類(lèi)似于CPU,矩陣乘法的性能在GPU中隨著矩陣尺寸增大而增大;同時(shí),GPU的浮點(diǎn)運(yùn)算能力一般遠(yuǎn)大于單個(gè)CPU.因此,在異構(gòu)并行優(yōu)化中,我們傾向于盡可能在GPU中執(zhí)行大尺寸的矩陣操作.從存儲(chǔ)方面考慮,為了避免GPU顯存和內(nèi)存之間頻繁的數(shù)據(jù)通信,多次參與GPU計(jì)算的數(shù)據(jù)需要存儲(chǔ)在GPU顯存中,主要包括Davidson方法中的向量、算符數(shù)據(jù)和臨時(shí)數(shù)據(jù)(算法2中但同時(shí),保留較多的狀態(tài)也導(dǎo)致計(jì)算中需要的存儲(chǔ)容量較大;考慮到當(dāng)前GPU顯存容量較小,在異構(gòu)并行方法中需要合理分配各個(gè)部分的存儲(chǔ)利用.

首先考慮Davidson方法的異構(gòu)并行實(shí)現(xiàn),如算法3所示,Davidson方法中主要執(zhí)行各種向量操作,相對(duì)而言運(yùn)算量較小,因此應(yīng)盡可能充分發(fā)揮內(nèi)存和GPU顯存的存儲(chǔ)帶寬.具體地,我們將所有向量以相同的方式按行劃分為兩部分,使得一部分波函數(shù)基矢對(duì)應(yīng)的分量存儲(chǔ)在GPU顯存中,另一部分存儲(chǔ)在內(nèi)存中.這樣任一向量的操作將由CPU和GPU共同完成,并且不需要向量數(shù)據(jù)的通信.當(dāng)GPU顯存足夠時(shí)(通常內(nèi)存容量遠(yuǎn)大于GPU顯存),內(nèi)存和GPU顯存中向量行數(shù)的比值為兩者存儲(chǔ)帶寬的比值,理論上此時(shí)可以獲得最高的性能.由于Davidson方法的存儲(chǔ)量線性依賴(lài)于子空間中向量的個(gè)數(shù)和超塊希爾伯特空間的維數(shù),當(dāng)向量個(gè)數(shù)或者DMRG保留狀態(tài)較大時(shí),GPU顯存中的向量數(shù)據(jù)會(huì)受到GPU顯存容量的限制.這會(huì)導(dǎo)致較多向量操作在CPU中執(zhí)行,此時(shí)Davidson方法獲得的性能較低.

接下來(lái)給出哈密頓量作用于波函數(shù)(即算法2)部分的異構(gòu)并行實(shí)現(xiàn),首先將其中的操作合理劃分到CPU和GPU中.將算法2中所有操作按照子塊S中的好量子數(shù)分為多個(gè)組,此時(shí)各個(gè)組中的運(yùn)算可以獨(dú)立進(jìn)行,僅在求和計(jì)算(即算法2中step3)時(shí)需要通信.此時(shí)波函數(shù)〉也 按S中好量子數(shù)劃分為兩個(gè)部分,其中一部分僅在GPU中計(jì)算(記為另一部分僅在CPU中計(jì)算(記為

算法2中耗時(shí)較多的運(yùn)算為一系列矩陣乘法,且S中每個(gè)好量子數(shù)對(duì)應(yīng)分組的計(jì)算量在執(zhí)行運(yùn)算前可以比較準(zhǔn)確地估計(jì).通常,GPU的浮點(diǎn)運(yùn)算能力強(qiáng)于CPU,適合處理大尺寸矩陣乘法運(yùn)算,因此在這一步盡量將大矩陣運(yùn)算分配至GPU,將相對(duì)較小的矩陣運(yùn)算分配至CPU執(zhí)行.為了實(shí)現(xiàn)這一目標(biāo),在具體操作中,我們將矩陣乘法平均運(yùn)算量較大的組分配到GPU中執(zhí)行,這里平均運(yùn)算量為組內(nèi)矩陣乘法總計(jì)算量與矩陣乘法個(gè)數(shù)的比值.進(jìn)一步,可根據(jù)GPU中計(jì)算量占的比例PGPU將相互獨(dú)立的分組計(jì)算分配給GPU,剩余組的計(jì)算分配給CPU執(zhí)行.

為了獲得較高的異構(gòu)并行效率,在執(zhí)行Davidson對(duì)角化方法時(shí)動(dòng)態(tài)調(diào)整PGPU以盡可能實(shí)現(xiàn)負(fù)載均衡.在具體操(作中,將)每一步迭代優(yōu)化比例所處的區(qū)間記為P0,P1.設(shè)定初始區(qū)間P0=0,P1=1,然后進(jìn)入Davidson(方法迭代)過(guò)程.令GPU中計(jì)算量比例PGPU=P0+P1/2,執(zhí)行一步Davidison迭代,可以獲得此時(shí)作用哈密頓量到波函數(shù)的CPU和GPU計(jì)算時(shí)間,分別記為T(mén)CPU和TGPU.以此為依據(jù)更新下一步迭代區(qū)間,使得

并開(kāi)始下一步Davidson迭代.如此按照類(lèi)似二分法的收斂思路,通過(guò)少數(shù)幾步迭代就可以收斂到優(yōu)化比例,之后的大部分迭代運(yùn)算都趨于負(fù)載均衡.

本文主要針對(duì)保留狀態(tài)數(shù)較大的問(wèn)題進(jìn)行優(yōu)化,此時(shí)涉及到的矩陣較大,GPU中矩陣運(yùn)算的并行效率較高.同時(shí),考慮到GPU顯存的限制,為了減小臨時(shí)數(shù)據(jù)的存儲(chǔ)量,根據(jù)子塊S分組后各個(gè)組的操作依次執(zhí)行.這種情況下,僅需要分配一段存儲(chǔ)空間,使其同時(shí)滿足任意一個(gè)組中的所有操作即可.在圖4中,給出了對(duì)角化哈密頓量部分運(yùn)算的存儲(chǔ)需求,并給出了與兩子塊表示所需存儲(chǔ)的對(duì)比.可以明顯看出其總體的顯存需求遠(yuǎn)遠(yuǎn)小于兩子塊表示.因此,相比于參考文獻(xiàn)[30]中的異構(gòu)并行,本文方案可以處理需要更大DMRG保留狀態(tài)數(shù)的問(wèn)題.

圖4 存儲(chǔ)臨時(shí)數(shù)據(jù),子塊算符需要的GPU顯存Fig.4.The GPU memory cost of temporary data and subblock operators.

4 數(shù)值結(jié)果

為了說(shuō)明本文優(yōu)化方法的有效性,我們分別保留4096,6144,8192,10240個(gè)狀態(tài)計(jì)算4腿Hubbard梯子的基態(tài)能量,并得到了相應(yīng)的性能基準(zhǔn).圖5(a)給出了DMRG優(yōu)化中各個(gè)部分相對(duì)于單個(gè)CPU的加速比(單個(gè)CPU計(jì)算時(shí)間與單個(gè)GPU計(jì)算時(shí)間的比值),可以看到作用哈密頓量在波函數(shù)部分獲得加速比最大,其加速比在保留狀態(tài)數(shù)大于4096后較為接近(不低于3.8).當(dāng)保留狀態(tài)數(shù)較大時(shí),由于GPU顯存總量的限制(如圖5(b)所示),大部分向量操作由CPU完成,這導(dǎo)致Davidson方法的加速比在保留狀態(tài)數(shù)較大時(shí)有所下降.然而,Davidson方法中向量操作占哈密頓量對(duì)角化的時(shí)間比例較少,因此對(duì)哈密頓量對(duì)角化部分加速比影響較小(不低于3.5).本文哈密頓量對(duì)角化之外的其他操作計(jì)算時(shí)間占總時(shí)間比例約10% (如圖3),該部分的GPU并行優(yōu)化還沒(méi)有被考慮,因此總并行加速比低于哈密頓量對(duì)角化部分的加速比,這里保留最大10240個(gè)狀態(tài)獲得的加速比為2.85.圖5(c)中給出了異構(gòu)并行實(shí)現(xiàn)中CPU和GPU分別貢獻(xiàn)的浮點(diǎn)性能,可以看到隨著保留狀態(tài)數(shù)的增大,CPU和GPU中執(zhí)行的大尺寸矩陣增多,兩者貢獻(xiàn)的浮點(diǎn)性能隨之增大.雖然本文數(shù)值計(jì)算保留狀態(tài)數(shù)較大,但是其中矩陣尺寸為兩子塊表示時(shí)的 1 /d(Hubbard模型中d=4),目前獲得的性能仍明顯小于GPU可以達(dá)到的峰值性能(1200 GFLOPS).對(duì) 1 6×4 的梯子,根據(jù)不同保留狀態(tài)數(shù)(4096,6144,8192,10240)得到的基態(tài)能量外推得到模型(1)在U=8.0 時(shí)的格點(diǎn)平均能量Eg=-0.75114(2),與文獻(xiàn)[34]結(jié)果一致,見(jiàn)圖6.進(jìn)一步給出基態(tài)的電荷密度分布,如圖7所示,可以觀察到明顯的電荷密度條紋,這是銅氧化物高溫超導(dǎo)體中經(jīng)常被觀測(cè)到的現(xiàn)象之一[13,14,37,38].

圖5 異構(gòu)并行的性能 (a)加速比;(b)Davidson方法中的向量占用GPU顯存;(c)作用哈密頓量到波函數(shù)部分的性能Fig.5.Performance of hybrid parallel strategy：(a)The speedup;(b)the GPU memory cost of vectors in Davidson;(c)the performance ofH|ψ〉.

圖6 基態(tài)能量關(guān)于截?cái)嗾`差的函數(shù)(直線表示對(duì)基態(tài)能量的線性外推,直至截?cái)嗾`差為0)Fig.6.Groundstate energy as a function of truncation error.The straight line gives a linear extrapolation of the ground energy until 0 truncation-error.

圖7 對(duì)于16×4 Hubbard模型,U=8.0時(shí)的基態(tài)電荷密度分布(可以觀察到明顯的電荷密度條紋)Fig.7.Ground state density profile for the 16×4 Hubbard ladder withU=8.0.Charge density stripes can be clearly observed.

5 結(jié) 論

本文主要考慮DMRG方法在準(zhǔn)二維格點(diǎn)模型中的應(yīng)用,針對(duì)其中最耗時(shí)的哈密頓量對(duì)角化部分實(shí)現(xiàn)了CPU-GPU異構(gòu)并行優(yōu)化,并且給出了負(fù)載平衡方法.為了減小準(zhǔn)二維格點(diǎn)模型計(jì)算中GPU顯存的限制,本文的實(shí)現(xiàn)中哈密頓量與波函數(shù)基于四子塊表示,其對(duì)角化時(shí)需要的GPU顯存占用遠(yuǎn)小于兩子塊表示,使得本文的異構(gòu)并行方法可以應(yīng)用于更多模型、更多問(wèn)題的研究.將該方法應(yīng)用到4腿Hubbard梯子模型的求解中,得到了不同保留狀態(tài)數(shù)時(shí)DMRG中各個(gè)部分的加速比.數(shù)值結(jié)果表明,本文的異構(gòu)并行方法適用于保留狀態(tài)數(shù)較大的準(zhǔn)二維模型計(jì)算,并且總性能隨著保留狀態(tài)數(shù)增大而增大.目前,強(qiáng)關(guān)聯(lián)物理問(wèn)題很大程度上依賴(lài)于多體數(shù)值計(jì)算,一些復(fù)雜問(wèn)題通常進(jìn)一步受制于計(jì)算方法的計(jì)算量與計(jì)算時(shí)間.在多體算法本身出現(xiàn)革命性發(fā)展之前,合理利用計(jì)算機(jī)技術(shù)的發(fā)展提升算法的效率能為研究強(qiáng)關(guān)聯(lián)系統(tǒng)提供很大的幫助.我們希望該并行方法可以在更多的復(fù)雜格點(diǎn)模型、更多問(wèn)題中得到應(yīng)用,并能夠進(jìn)一步引起強(qiáng)關(guān)聯(lián)領(lǐng)域?qū)τ谝訥PU為代表的新技術(shù)的關(guān)注和重視.