張 亮，常 旭，秦志楷，沈 立

(國防科技大學計算機學院，湖南 長沙 410073)

1 引言

量子計算是一種全新的計算模式，它以量子位為信息單元，遵循量子力學的規(guī)律，通過調(diào)控量子信息單元來進行計算[1]。量子計算機以量子力學為基礎(chǔ)，其理論模型是量子圖靈機。由于量子力學態(tài)存在疊加原理，量子信息的處理更加迅速，與通用電子計算機相比具有更強大的信息處理能力[2]，在合適的量子算法支持下，能夠以很低的時間開銷來完成密碼學、多體系統(tǒng)量子力學和量子機器學習[3]等領(lǐng)域的諸多復(fù)雜算法。

然而，量子計算機的實現(xiàn)非常困難，目前仍沒有實用化的量子計算機可以用來求解大規(guī)模的科學計算問題[4]，因此使用電子計算機來模擬量子計算系統(tǒng)成為最常用的研究手段。許多量子線路模擬器應(yīng)運而生，用來支撐量子領(lǐng)域復(fù)雜計算問題的模擬，如ProjectQ、qHiPSTER和QuEST等。這些模擬器作為重要的研究工具，受到了量子計算研究者的高度重視。但是，目前開發(fā)出的量子線路模擬器往往會消耗更多的存儲資源，帶來更多的計算成本。

ProjectQ[5]是一款用于模擬量子計算的開源軟件，它的初始版本就具有針對不同硬件的編譯器框架，該框架可以對量子算法進行仿真測試，并可以連接到IBM Quantum Experience云服務(wù)的后端，在已有的量子硬件上運行。

qHiPSTER[6]是一款在電子計算機上實現(xiàn)的高性能分布式量子線路模擬器，可以模擬常規(guī)的單量子位門和雙量子位控制門，并針對單CPU和多CPU做出了優(yōu)化，包括向量化、多線程和高速緩存等。該模擬器既可以實現(xiàn)高性能模擬，又可以達到較高的硬件利用率，僅受到內(nèi)存的限制[7]。

QuEST量子線路模擬器[8]的全稱為量子精確模擬工具包(the Quantum Exact Simulation Toolkit)，是第1款開源、采用OpenMP和MPI混合編程開發(fā)、支持單GPU加速的通用量子線路模擬器；從日常生活中的筆記本電腦到世界上功能最強大的超級計算機，它在這些經(jīng)典計算平臺上都可以實現(xiàn)量子計算的模擬。QuEST能夠在純態(tài)和混合態(tài)上、在有退相干的情況下，用狀態(tài)向量和密度矩陣表示量子狀態(tài)，模擬由一般的單量子位門、雙量子位門和多量子位控制門組成的通用量子線路[9]。

雖然目前已經(jīng)有了幾款開源的量子線路模擬器，但它們普遍都存在資源開銷大、模擬時間長等問題，而且不能很好地利用GPU等計算加速器，有的完全沒有GPU版本，有的只能在單GPU[10]上運行,有的雖然有多GPU版本，但不具有擴展性和多平臺移植性[11]。本文以QuEST為例，分析了其代碼結(jié)構(gòu)和性能特性，并實現(xiàn)了QuEST模擬器在多GPU平臺上的性能優(yōu)化。測試結(jié)果表明，在4塊NIVIDIA V100 GPU構(gòu)成的平臺上，相較于單CPU版本可獲得7～9倍的性能加速，相較于多CPU版本可獲得3倍的性能加速。

2 QuEST分析

在量子系統(tǒng)中，量子狀態(tài)不僅可以是二進制1或0，而且可以是同時持有1和0的多種組合，從而形成一個“量子位”。如式(1)和式(2)所示，一個量子位的狀態(tài)可以用一個二維的狀態(tài)向量來表示，其中α和β都是復(fù)數(shù)。多個量子位的表示同理，只要增加狀態(tài)向量的維數(shù)即可。

|φ〉=α|0〉+β|1〉

(1)

(2)

量子門是在一組量子位之間進行操作的基本量子線路，是量子線路的基礎(chǔ)[12]。在進行量子線路模擬時，多個量子門的運算是串行完成的，每個量子門都會對全部的狀態(tài)向量進行一次遍歷，完成狀態(tài)向量的更新計算。QuEST模擬器實現(xiàn)了tGate、sGate、hadamard、RotateX/Y/Z和PauliX/Y/Z等14種量子門。

2.1 算法分析

QuEST模擬器記錄了被模擬量子系統(tǒng)的狀態(tài)，以及各個進程的相關(guān)信息。它用一個復(fù)數(shù)數(shù)組來存放全部狀態(tài)向量，每個狀態(tài)向量由2個double類型的浮點數(shù)構(gòu)成，分別表示狀態(tài)向量的實部和虛部，是QuEST模擬器進行計算的基本數(shù)據(jù)單元。狀態(tài)向量的總數(shù)與被模擬的量子位數(shù)有關(guān)，假設(shè)對擁有N個量子位的量子系統(tǒng)進行模擬，就需要使用2N個狀態(tài)向量來進行表示。例如，當N=30時，2N個狀態(tài)向量共需要16 GB內(nèi)存空間。被模擬的量子系統(tǒng)每增加1個量子位，狀態(tài)向量需要的存儲空間就會隨之翻1倍，可見QuEST模擬器的內(nèi)存消耗很大。

為了提高模擬速度，QuEST模擬器實現(xiàn)了CPU分布式版本，可以使用多個進程在CPU平臺上進行量子系統(tǒng)的并行模擬。在并行模擬的初始化階段會對狀態(tài)向量分組，組數(shù)與進程數(shù)相同，每個進程維護1組私有狀態(tài)向量，并負責它們的更新。對于分布式版本，每個進程還需要1個額外的狀態(tài)向量數(shù)組，專門用于進程間的數(shù)據(jù)通信，其大小與進程私有的狀態(tài)向量數(shù)組相同，這就導(dǎo)致使用多個進程進行并行量子線路模擬時，存儲空間消耗更加巨大。

QuEST源碼的核心部分是實現(xiàn)了14個基本量子門的模擬函數(shù)，這些模擬函數(shù)是基于2個模板實現(xiàn)的。tGate、sGate和PauliZ使用了同一個模板，其特點在于根據(jù)輸入?yún)?shù)來決定要更新的數(shù)據(jù)，不存在進程間數(shù)據(jù)交換。其余函數(shù)使用另一種模板，根據(jù)輸入?yún)?shù)找到對應(yīng)的進程，先進行數(shù)據(jù)交換，再完成數(shù)據(jù)更新。

函數(shù)statevec_compactUnitary()是RotateX/Y/Z、PauliX/Y等量子門的核心函數(shù)，下面以該函數(shù)為例來分析QuEST模擬器的核心計算過程，找出其性能上的不足。

該函數(shù)的工作是先將所有狀態(tài)向量劃分為2部分，索引標記分別為a和b，二者之間的差值由輸入?yún)?shù)t決定(滿足關(guān)系：2t=a-b)，然后利用a與b對應(yīng)的狀態(tài)向量值，使用規(guī)定的算法進行更新后得到新的狀態(tài)向量值，每個狀態(tài)向量值的更新都需要同時使用到a與b對應(yīng)的向量數(shù)據(jù)。

當使用多個進程進行分布式模擬時，如果a與b對應(yīng)的狀態(tài)向量屬于不同進程，那么先要找到與本地對應(yīng)的進程(這個過程稱為進程配對)，配對的進程間進行通信，將對方的狀態(tài)向量拷貝到本地，然后再繼續(xù)狀態(tài)向量的更新工作；如果a與b對應(yīng)的狀態(tài)向量都在本地，那么就不需要配對和通信，直接進行狀態(tài)向量的更新工作。該計算過程的時間開銷很大，原因在于每個門函數(shù)的計算都需要將所有狀態(tài)向量遍歷1次，數(shù)據(jù)局部性很差；同時，進程間的通信也會帶來很大的時間開銷。

值得一提的是，在執(zhí)行具體的門函數(shù)時，同一門函數(shù)的核心計算部分循環(huán)內(nèi)部并沒有數(shù)據(jù)相關(guān)，可以將循環(huán)完全并行化處理，這正是使用單GPU版本可以獲得明顯加速效果的原因。但是，由于顯存限制，單GPU版本并不能滿足量子位更多的(例如N>30時)量子系統(tǒng)的模擬。

2.2 QuEST的性能特征

通過對QuEST源碼的算法分析，可以發(fā)現(xiàn)QuEST是一款計算密集型的應(yīng)用軟件，具有良好的可擴展性和并行性，并且內(nèi)存消耗巨大。本節(jié)以30個量子位的量子系統(tǒng)為例，使用Random與GHZ_QFT 2個量子線路作為測試算例，分別選取單CPU、多CPU和單GPU這幾種不同的測試平臺，對QuEST模擬器的性能進行測試。

Random測試算例是隨機生成的，從Random測試算例的部分代碼中可以看出，Random程序調(diào)用的門函數(shù)是無序、隨機的，而GHZ_QFT測試算例模擬的是實際量子線路。

Random測試算例的部分代碼如下所示：

Quregq=createQureg(30,Env);

floatq_measure[30];

tGate(q,25);

controlledNot(q,28,21);

controlledRotateX(q,17,5,0.3293660327520663);

rotateX(q,27,3.9207427542347926);

tGate(q,3);

controlledRotateZ(q,27,19,5.459935259485407);

controlledRotateX(q,26,3,4.736990305652013);

controlledRotateZ(q,8,11,3.594728080156504);

rotateX(q,10,4.734238389048838);

rotateY(q,8,4.959946047271496);

rotateZ(q,5,1.0427019597472071);

?

對比Random測試算例的代碼與GHZ_QFT測試算例的代碼可以發(fā)現(xiàn)，GHZ_QFT程序中的門函數(shù)更加規(guī)整，相較于Random測試算例更具有實際的研究意義。

GHZ_QFT測試算例的部分代碼如下所示：

/* QFT starts */

hadamard(q,0);

controlledRotateZ(q,0,1,1.5708);

hadamard(q,1);

controlledRotateZ(q,0,2,0.785398);

controlledRotateZ(q,1,2,1.5708);

hadamard(q,2);

controlledRotateZ(q,0,3,0.392699);

?

測試環(huán)境如表1所示，每個計算結(jié)點有2塊Intel Xeon E5-2692 V2處理器(內(nèi)存64 GB)和2塊GPU。GPU分別選用NVIDIA Tesla K80(顯存12 GB)與NVIDIA Tesla V100(顯存32 GB)2種型號。測試結(jié)果分別如圖1與圖2所示。

Table 1 Hardware and environment of test experiment表1 測試實驗的硬件配置與運行環(huán)境

Figure 1 QUEST simulation time (30 bits)圖1 30位QUEST模擬時間

Figure 2 QUEST performance speedup (30 bits)圖2 30位QUEST性能加速比

圖1中橫軸表示每組測試所使用的計算平臺，縱軸表示每組測試的運行時間，單位為秒。

圖2中橫軸表示每組測試的加速比(以單CPU性能為基準)，縱軸表示每組測試所使用的計算平臺。

測試結(jié)果表明，對比單CPU與多CPU版本代碼的運行時間，可以確認QuEST可擴展性良好。受GPU內(nèi)存限制，30個量子位的量子系統(tǒng)無法在顯存僅有12 GB的NVIDIA Tesla K80上完成模擬，而使用32 GB顯存的NVIDIA Tesla V100運行單GPU版本，相較于單CPU版本可以獲得高達10.79倍和18.24倍的性能加速比，加速效果明顯。這說明QuEST的并行性極好，十分契合GPU的單指令多線程體系結(jié)構(gòu)，使用GPU設(shè)備來模擬運行再合適不過，但是對單個GPU的顯存容量要求很高。因此，本文設(shè)計并實現(xiàn)了QuEST模擬器的多GPU版本，突破了單個GPU內(nèi)存容量對QuEST模擬器的限制，并大大提高了模擬速度。

3 多GPU版本的實現(xiàn)

3.1 總體框架

1個計算結(jié)點通常包含多個GPU設(shè)備，QuEST模擬器在MPI啟動時可以指定結(jié)點數(shù)和每個結(jié)點的進程數(shù)，如果采用每個GPU設(shè)備綁定1個MPI進程的設(shè)置，當程序處理的數(shù)據(jù)規(guī)模變大，需要增加GPU時，只需要在程序啟動時修改MPI的初始化參數(shù)即可，這樣程序就具有較強的可擴展性。反之，如果每個結(jié)點綁定1個MPI進程，即1個MPI進程管理多個GPU設(shè)備，程序的可擴展性就會大幅度下降。當結(jié)點內(nèi)的GPU數(shù)量變化時，需要對應(yīng)地修改程序，程序的可擴展性較差。

本文對QuEST的單GPU實現(xiàn)進行了修改。首先為每個GPU設(shè)備綁定1個MPI進程；其次在初始化階段根據(jù)模擬的量子位數(shù)為GPU分配相應(yīng)的內(nèi)存空間，狀態(tài)向量按照從前至后的順序分配到多個GPU設(shè)備上；然后將單GPU版本的量子門模擬函數(shù)修改為多GPU版本；最后完成多GPU之間的數(shù)據(jù)規(guī)約操作，得到正確的計算結(jié)果。

門函數(shù)的多GPU版本實現(xiàn)參考了CPU分布式版本的代碼框架，這里仍舊以statevec_compactUnitary()這個核心計算函數(shù)為例，當使用多個GPU時，初始化階段完成狀態(tài)向量的分組，組數(shù)與進程數(shù)一致；核心函數(shù)將狀態(tài)向量分為α與β2部分，并判斷2部分數(shù)據(jù)是否屬于同一進程；若是，直接進行狀態(tài)向量的更新計算；若否，進程配對后進行GPU間通信，將彼此的數(shù)據(jù)拷貝到對方的GPU顯存中，然后進行更新計算。其他門函數(shù)的實現(xiàn)過程與此類似，不再贅述。

3.2 多GPU通信

使用多個GPU并行執(zhí)行應(yīng)用程序時，需要合理地選擇GPU間的通信方式(或組合)，這對GPU間的數(shù)據(jù)傳輸效率具有非常大的影響[13]。多GPU系統(tǒng)提供了2種GPU連接方式：(1)單個結(jié)點內(nèi)的所有GPU連接到該結(jié)點的PCIe總線上；(2)結(jié)點間的多個GPU通過集群中的網(wǎng)絡(luò)交換機連接在一起。這2種連接方式并不互斥，其拓撲結(jié)構(gòu)如圖3所示。

Figure 3 Topology of multi-GPU system圖3 多GPU系統(tǒng)的拓撲結(jié)構(gòu)

點對點(P2P)傳輸可以實現(xiàn)GPU設(shè)備間的直接通信，對于連接在同一條PCIe總線上的多個GPU而言，這要比使用主機內(nèi)存中轉(zhuǎn)的通信方式高效得多。不同結(jié)點上的GPU設(shè)備也可以通過互連網(wǎng)絡(luò)來完成點對點通信，現(xiàn)有的MPI庫支持CUDA-aware，實現(xiàn)了對GPU顯存的直接訪問。

由于采取1個進程綁定1個GPU的設(shè)計思想，所以進程間的通信問題就轉(zhuǎn)換為GPU間的數(shù)據(jù)交換問題。為每個進程設(shè)置了2個數(shù)組：一個用來存放本地的狀態(tài)向量數(shù)據(jù)；另一個用來存放從配對進程拷貝過來的狀態(tài)向量數(shù)據(jù)。并直接使用MPI庫中的消息傳遞函數(shù)，在互相配對的2個進程間傳輸向量數(shù)據(jù)，以這樣的方式來實現(xiàn)多個GPU間的通信。

此外，在實際對源碼修改的過程中，發(fā)現(xiàn)MPI庫中的消息傳遞函數(shù)無法使用“指針+偏移量”的方式來直接訪問GPU設(shè)備上的數(shù)據(jù)，因而無法在GPU間分段傳輸數(shù)據(jù)，因此本文做出了進一步改進，使用“以空間換時間”的方法解決了此問題。具體步驟為：(1)為每一個進程再分配2個緩沖區(qū)B1和B2，用于完成GPU間的消息傳遞，其中，B1用于向配對進程發(fā)送數(shù)據(jù)，B2用于接收來自配對進程的數(shù)據(jù)；(2)實現(xiàn)2個Kernel函數(shù)，一個將數(shù)組內(nèi)的部分數(shù)據(jù)寫入到B1內(nèi)，另一個將B2內(nèi)的數(shù)據(jù)寫回到數(shù)組內(nèi)，從而實現(xiàn)了本地與這2個緩沖區(qū)之間數(shù)據(jù)的分段傳遞。

4 性能測試與分析

本節(jié)仍然采用與第2節(jié)相同的環(huán)境和算例進行性能測試。

4.1 模擬時間

本節(jié)首先對30個量子位的量子系統(tǒng)設(shè)計了一組測試實驗，對多GPU版本的模擬時間進行測試。仍使用Random與GHZ_QFT 2個測試算例，分別在單CPU、多CPU、單GPU和多GPU 4個平臺上進行測試，多CPU和多GPU的數(shù)量均為4，實驗結(jié)果分別如圖4與圖5所示。

Figure 4 Comparison of simulation time between multi-GPU version and other versions (30 bits)圖4 多GPU版本與其它版本模擬時間比較(30位)

圖4中橫軸表示每組測試所使用的計算平臺，縱軸表示每組測試的運行時間，單位為s。

Figure 5 Comparison of performance speedup between multi-GPU version and other versions (30 bits)圖5 多GPU版本與其它版本性能加速比比較(30位)

圖5中橫軸表示每組測試的加速比(以單CPU性能為基準)，縱軸表示每組測試所使用的計算平臺。

通過數(shù)據(jù)對比可以發(fā)現(xiàn)，對于2個測試算例而言，多GPU版本的代碼均可以獲得明顯的加速效果。對于Random算例，多GPU版本的運行時間為34.1 s，相較于單CPU版本獲得的加速比高達9.09倍，相較于多CPU版本獲得的加速比大約為3.19倍；而對于GHZ_QFT算例，多GPU版本的運行時間為32.5 s，相較于單CPU版本獲得的加速比高達7.04倍，相較于多CPU版本獲得的加速比大約為3.76倍。

但是，對比單GPU版本與多GPU版本的運行時間可以發(fā)現(xiàn)，單GPU版本的加速效果要優(yōu)于多GPU版本的。這主要是因為多GPU版本在運行時，GPU之間需要進行通信，從而產(chǎn)生了額外的開銷，導(dǎo)致速度降低，這是不可避免的損失，而單GPU版本不涉及此類問題，所以加速效果更加明顯。換句話說，在模擬的量子位較少，GPU顯存充足的情況下，還是首選單GPU版本來進行量子系統(tǒng)的模擬。同時，雖然多CPU版本也存在通信問題，但其計算收益遠遠大于通信帶來的額外開銷，所以相較于單CPU版本可以獲得明顯的性能提升。

30個量子位的量子系統(tǒng)使用單GPU模擬需要消耗16 GB內(nèi)存空間，多GPU版本的總內(nèi)存消耗為48 GB，即每個GPU消耗12 GB。由于B1與B2緩沖區(qū)的實現(xiàn)，相較于單GPU版本，多GPU版本引入了4 GB的額外內(nèi)存開銷。

4.2 量子位擴展

本節(jié)設(shè)計了另一組測試實驗，仍舊使用Random與GHZ_QFT 2個測試算例，將被模擬系統(tǒng)的量子位提高到31位，同樣分別在單CPU、多CPU、單GPU和多GPU 4個平臺上進行測試，測試結(jié)果分別如圖6與圖7所示。

Figure 6 Comparison of simulation time between multi-GPU version and other versions (31 bits)圖6 多GPU版本與其它版本模擬時間比較(31位)

圖6中橫軸表示每組測試所使用的計算平臺，縱軸表示每組測試的運行時間，單位為s。

Figure 7 Comparison of performance speedup between multi-GPU version and other versions (31 bits)圖7 多GPU版本與其它版本性能加速比比較(31位)

圖7中橫軸表示每組測試的加速比(以單CPU性能為基準)，縱軸表示每組測試所使用的計算平臺

通過數(shù)據(jù)對比可以發(fā)現(xiàn)，對于2個測試算例而言，多GPU版本的代碼均可以獲得明顯的加速效果。對于Random算例，多GPU版本相較于單CPU版本獲得的加速比高達13.45倍，相較于多CPU版本獲得的加速比大約為5.04倍；而對于GHZ_QFT算例，多GPU版本的運行時間為10.7 s，相較于單CPU版本獲得的加速比高達42.61倍，相較于多CPU版本獲得的加速比大約為18.76倍。

然而，由于顯存不足，對于這2個算例，單個同樣型號的GPU均無法完成模擬測試，這說明改進后的多GPU版本代碼不僅能夠帶來明顯的性能加速，而且還確實解決了GPU顯存不足的問題。此外，GHZ_QFT算例模擬的是更加規(guī)整的實際量子線路，相對于隨機生成的Random算例，它取得的加速效果更加顯著，也更加具有實際研究意義。

4.3 結(jié)論

通過上述實驗，可以得出有關(guān)多GPU 版本QuEST量子線路模擬器的3個結(jié)論：

(1)多GPU版本代碼相對于CPU代碼能夠大大加速程序運行，加速效果顯著；

(2)多GPU版本代碼能夠模擬具有更多量子位的量子系統(tǒng)；

(3)多GPU版本代碼在通信較少的情況下，能夠得到比單GPU更好的加速效果。

5 結(jié)束語

本文首先詳細分析了QuEST量子線路模擬器的算法和可擴展性設(shè)計。QuEST目前提供單CPU、多CPU和單GPU 3個版本，雖然單GPU版本的QuEST模擬速度非?？?，但是當模擬的位數(shù)增加到31位時，模擬所需要的內(nèi)存已經(jīng)超過目前市面上內(nèi)存最大的GPU，所以單GPU版本最多只能模擬30個量子位的系統(tǒng)。當量子位數(shù)超過30位時，只能使用多CPU版本進行模擬。為了使QuEST能夠模擬更多的位數(shù)同時速度盡可能地快，本文開發(fā)了多GPU版本，能夠模擬超過30位的量子算法，解決了單GPU內(nèi)存不足的問題。測試結(jié)果表明，與多CPU版本相比，本文開發(fā)的多GPU版本相較于單CPU版本能夠獲取7～9倍的性能提升，相較于多CPU版本可以獲得3倍的性能提升。

未來將把工作重點放在優(yōu)化GPU代碼上，比如使用各種CUDA的優(yōu)化技術(shù)來對QuEST程序進行加速，使其獲得更好的加速效果，也會嘗試使用RDMA技術(shù)優(yōu)化結(jié)點間的通信，使其通信效率更高；同時，本文工作并沒有修改QuEST模擬器的框架，但目前已經(jīng)有一些工作利用張量網(wǎng)絡(luò)和張量積解決量子線路模擬器內(nèi)存消耗過大的問題[14-17]，如何使用這一方法優(yōu)化QuEST將是接下來的一個研究方向。