基于GPU的約束網絡模型和并行弧相容算法

李 哲 李占山 李 穎

(吉林大學計算機科學與技術學院 長春 130012) (符號計算與知識工程教育部重點實驗室(吉林大學) 長春 130012)(zslizsli@163.com)

基于GPU的約束網絡模型和并行弧相容算法

李 哲 李占山 李 穎

(吉林大學計算機科學與技術學院 長春 130012) (符號計算與知識工程教育部重點實驗室(吉林大學) 長春 130012)(zslizsli@163.com)

弧相容算法是約束滿足問題的基本壓縮求解空間算法之一，很多優秀的高級算法都以高性能的弧相容算法作為核心.近年來，以GPU為計算工具加速并行計算被用來嘗試解決許多問題.基于GPU和基本的并行算法，提出一種適合GPU運算的約束網絡表示模型N-E，給出其生成算法BuildNE.結合細粒度的弧相容算法——AC4，基于N-E模型提出AC4的并行化算法AC4GPU與改進算法AC4GPU+，使弧相容算法得以擴展到GPU上執行.實驗結果驗證了該算法的可行性，與AC4算法的比較，其在一些規模較小的問題上取得了10%～50%的加速，在一些規模較大的問題上則加速1～2個數量級.為今后進一步在GPU上以并行形式解決其他約束滿足問題提供了一種核心算法方案.

人工智能；約束滿足問題；弧相容；圖形處理器；計算統一設備架構

約束滿足問題(constraint satisfaction problem, CSP)[1]是人工智能領域的一個重要分支，現實生活中的很多問題，都可以約束建模轉化為約束滿足問題，使用約束求解技術進行求解.如生產調度問題、路由選擇問題、產品配置問題、診斷等.眾所周知，弧相容(arc consistency, AC)是最著名的約束傳播算法之一，是求解約束滿足問題的核心，在維持弧相容(maintain arc consistency, MAC)算法中廣泛使用，也是執行單弧相容(singleton arc consistency, SAC)算法[2]的底層算法,它確保每個變量值域中的每個值與每個約束是相容的.Mohr和Henderson[3]在1986年提出了AC4[4]，是第1個細粒度的弧相容算法，已經證明AC4具有最優的時間復雜度.

伴隨多核的硬件大量出現，并行理論迅速發展，逐漸產生出一大批并行算法，通過對原有算法的并行化或者提出全新的并行算法都大幅提升了算法的性能.從早先Nguyen[5]，Yokoo等人[6]提出利用分布式來解決特定的約束滿足問題，到 Campeotto等人[7]利用混合CPU和GPU進行約束傳播，在某些測試用例上取得了優勢.很多學者對并行約束滿足問題進行了嘗試，涵蓋并行、分布式[5-7]、多核、GPU等方向[8-9].圖形處理器(GPU)由于在處理大規模的并行計算領域的優點，給計算領域帶來巨大變革.現有的GPU并行計算框架，如CUDA，OpenCL，C++ AMP等已建立起一套成熟的編程模型，可以用來快速地解決目標問題.

推理和搜索是處理約束滿足問題的2種技術，現有基于GPU的約束滿足問題的研究主要集中在搜索方向，通過細粒度的GPU級別的并行機制處理約束傳播[8]；而在推理方向的研究則較少，以現有的弧相容算法為例，它們大都基于串行的思想提出，難以在大規模細粒度的并發程序中取得令人滿意的效果.本文嘗試采用推理的方法對約束網絡的轉化問題進行研究，提出將約束網絡轉換為更易于并行算法的新模型——N-E模型；并參考AC4算法，提出在該模型上基于GPU運行的并行弧相容算法AC4GPU及其改進算法AC4GPU+，進而分析了GPU并行算法的特點.最后實驗結果表明：這2個算法在一些規模較小的問題較之AC4取得了10%～50%的加速，在一些規模較大的問題上，取得了1～2個數量級的性能提高，并為單弧相容算法和求解約束滿足問題提供了一種新的可行底層算法.

本文提出的在GPU上解決弧相容的方案包括3部分：

1) 提出一種約束網絡(constraint network)在GPU上表示的模型——N-E模型；

2) 提出基于GPU創建N-E模型的并行算法BuildNE；

3) 提出在N-E模型上基于GPU的細粒度弧相容算法AC4GPU及其的改進算法AC4GPU+.

1 背景知識

定義1. 變量和約束.變量(variable)[10]是具有不同取值的對象.變量x擁有與其相關聯的論域，用dom(x)表示，該論域包括被x允許值的集合，dominit(x)表示初始時(未刪除變量值前)變量x的論域.約束(constraint)[10]是定義在變量集合上的關系.一個約束c表示變量取值之間的制約關系，c所包含的變量是整個變量集合的子集，用scp(c)表示；與c相關聯的關系用rel(c)表示.

定義2. 約束網絡[1].約束網絡P=(X,C)，其中X={x0,x1,…,xn-1}是n個變量的集合，C={c0,c1,…,ce-1}是e個約束的集合.

給定約束網絡P，對于任意c∈C，滿足|scp(c)|=2，我們稱該約束網絡為二元約束網絡，若存在c∈C，滿足|scp(c)|>2,則稱該網絡為多元約束網絡，多元約束網絡可以轉化為等價的二元約束網絡.因此，二元約束滿足問題一直是該領域研究的熱點，本文主要討論二元約束滿足問題.

定義3. 弧相容[4,10].對于任意約束網絡P=(X,C)，弧相容可以被遞歸定義如下：

1) 1個點(x,a)是弧相容的，當且僅當在每一個包含x的約束c中，都存在(x,a)的支持，其中c∈C,x∈scp(c),a∈dom(x);

2) 1個變量x是弧相容的，當且僅當?a∈dom(x)，(x,a)是弧相容的；

3) 1個約束網絡P是弧相容的，當且僅當該網絡的每一個變量都是弧相容的.

例1. 圖1(a)所示為一個簡單的二元約束網絡P=(X,C)，其中：

1) 變量集合X={x1,x2,x3};

2) 約束集合C={c1,c2}，c1:x1=x2，c2:x1

圖1(b)刪去了約束網絡P中不滿足弧相容的點，此時整個約束網絡P滿足弧相容.

Fig. 1 Binary constraint network圖1 二元約束網絡

1.1 AC4

AC4是Mohr和Henderson[3]在1986年提出的一種弧相容算法.該算法分為2個階段：初始化階段和傳播階段.

1) 初始化階段計算計數器counter[xi,vi,xj]，其中xi,xj∈ci j,vi∈dom(xi)，表示根據約束ci j包含的一個變量xi中的點(xi,vi)在另一個變量xj上的支持數.算法同時根據ci j建立由(xi,vi)支持的所有點的列表S[xi,vi].本階段在約束集合上執行所有可能的約束檢查，每當點(xi,vi)在ci j上找到1個支持vj∈dom(xj).計數器counter[xi,vi,xj]加1，同時(xi,vi)被加到S[xi,vi]列表中.若在某約束上一個點在另一個變量的論域中沒有支持時，此時對應counter三元組的值為0，算法就將這個點從論域中刪去，并放入隊列Q中，等待在傳播階段中進行約束傳播.

2) 傳播階段針對Q中每個待刪除的點(xj,vj)，為每一個點(xi,vi)∈S[xj,vj]的計數器counter[xi,vi,xj]減1.若counter[xi,vi,xj]減到0，表明(xi,vi)在變量xj沒有支持，(xi,vi)被放入隊列Q中.當Q為空時表明沒有不相容的點存在，達成弧相容.若在刪點過程中發現某一變量論域為空，則不相容.

算法1. AC4.

輸入：P=(X,C)；

輸出：P滿足AC輸出真，反之輸出假.

①Q←?;S[xi,vi]←?,vj∈dom(xj),?xj∈X；

② foreachxi∈X,ci j∈C,vi∈dom(xi) do

③ 初始化counter[xi,vi,xj] to |{vj∈

dom(xj)|(vi,vj)∈ci j}|；

④ ifcounter[xi,vi,xj]=0 then

⑤ 刪除(xi,vi)并將之加入Q中；

⑥ end if

⑦ 將(xi,vi)加入S[xj,vj]中s.t.(vi,vj)∈ci j；

⑧ ifdom(xi)=? then

⑨ return false；

⑩ end if

2 基本并行算法及CUDA編程技術

2.1 歸 約

定義4. 歸約(reduce).歸約[11-12]是一類基本的并行算法，對傳入的O(n)個輸入數據，使用一個滿足結合率的二元運算⊕，生成1個結果.歸約是其他高級算法中重要的基礎算法.

例2. 由滿足結合律的二元運算⊕連接的8個元素的歸約結果可以用如下形式表示：

由于該運算滿足結合律，故其輸入元素可以任意順序進行二元計算.圖2展示了處理8元素數組的不同方式.為了方便對比，我們在圖2(a)中給出了串行實現，現在只需要具備可以計算⊕操作的計算單元，完成計算需要執行7步，性能較差.

在例2并行執行的過程中(如圖2(b)(c))，步驟1并發啟動4個線程計算相鄰或交替的2個元素，將計算結果寫入內存中；步驟2啟動2個線程計算第1步的結果；步驟3再啟動1個線程計算步驟2得到的2個結果.假設不考慮線程啟動的代價，對數步長的歸約得到計算結果只需要O(logn)步(這里是3步).

Fig. 2 Reduction of set elements圖2 集合元素的歸約

按鍵歸約(reduce by key)對于一組用〈k,v〉形式表示的鍵值對序列：

{〈k0,v00〉,…,〈k0,v0m〉,〈k1,v10〉,…,

〈k1,v1l〉,…,〈kp,vp0〉,…,〈kp,vp n〉}.

對具有相同鍵(key)的若干值(value)進行歸約求和，其結果為

對于輸入規模為O(n)的數組，對其進行按鍵歸約操作的最壞步驟復雜度為O(logn).

2.2 掃 描

掃描(scan)[12-13]，或稱前綴掃描(prefix scan)、前綴求和(prefix sum)、并行前綴求和(parallel prefix sum)是并行編程的重要術語，并作為基本模塊用于許多不同算法.

定義5. 掃描.掃描可以分為包容性掃描(inclusive scan)和排他性掃描(exclusive scan)兩種情況.對傳入的n個輸入數據{a0,a1,…,an-1}，使用滿足結合率的二元運算⊕，生成輸出序列：

{a0,(a0⊕a1),…,(a0⊕a1⊕…⊕an-1)}，

稱為包容性掃描(inclusive scan).

現定義恒等式元素I，使得I⊕a=a，生成輸入序列：

{I,a0,(a0⊕a1),…,(a0⊕a1⊕…⊕an-2)}，

稱為排他性掃描(exclusive scan).

可以看出，包容性和排他性掃描可以通過添加或刪去輸入數組的元素而相互轉化.

流壓縮(stream compaction)是根據一定標準篩選數組元素的操作，對輸入序列{a0,a1,…,an-1}的每個元素建立判定P，用來判定對應輸入元素是否應該包含在輸出數組中，這些判定P按順組成判定序列{P0,P1,…,Pn-1}，對此判定序列執行1次排他性掃描計算出滿足條件的輸入元素在輸出數組的位置.對于輸入數據空間復雜度為n的數組，流壓縮的步驟復雜度為O(logn).

2.3 CUDA編程技術

計算統一設備架構(compute unified device architecture, CUDA)是顯卡廠商NVIDIA推出的并行計算架構.該架構通過利用GPU的并發處理能力，可大幅提升計算性能.

NVIDIA提供一套完整的CUDA編程模型[12]，通過編程語言的接口訪問GPU的計算和存儲資源.內核(kernel)是在GPU上執行的程序，問題可以被分成若干部分，以線程塊的形式被GPU上不同的流處理器(stream multiprocessor, SM)并發執行.1個線程塊由多個線程組成，1個流處理器可以順序或并發地運行多個線程塊.線程塊內的線程具有相同的生命周期，所有的線程都可以訪問GPU內部的全局內存.

在CUDA中，流(stream)[12,14]可以使主機與設備(GPU)之間的內存操作與內核操作并發執行，并且支持在GPU上的多個內核并發執行，并支持多GPU并發執行.在下面AC4GPU算法使用了CUDA流，并發地執行2個不相干的任務，將約束傳播與約束網絡的相容性驗證并發執行，從而提升算法性能.CUDA也擁有優秀的并行模板庫——Thrust[15]，它采用STL的編程風格，用精簡的代碼實現了高性能的并行程序.本文實驗程序采用Thrust庫函數編寫.

3 并行算法的約束網絡模型

V=v_vals(P);

scp(c)={x1,x2,…,xr}∧(a1,a2,…,ar)∈

support(c)}.

V=v_vals(P);

scp(c)={x1,x2,…,xr}∧(a1,a2,…,ar)∈

conflict(c)}.

3.1 適用于GPU運算的約束網絡模型：N-E模型

定義6中我們引入了相容超圖和不相容超圖的概念，它們各自記錄約束網絡的一部分.根據定義6，本節提出一種在GPU上表達約束網絡的數據結構——N-E模型，該模型可以完整地表達約束網絡，并利用GPU擅于大規模并發處理數組元素的特點，使得約束網絡可以在此基礎上進行并行化推理.

定義7. 點集(nodes).在約束網絡P=(X,C)中，nodes表示網絡中所有點組成的集合及它們在約束傳播過程中的存在情況，可以形式化表示為

nodes={〈x,v,e〉|x∈X,v∈dom(x),

e=(v∈dom(x)dom(x))},

其中,e表示(x,v)在約束傳播過程中是否被刪除.nodes中的元素稱作node，是1個三元組，即node=〈x,v,e〉.

定義8.edge元組.在標準二元約束網絡P=(X,C)中，?c∈C，有?i,j使得xi,xj∈scp(c)，?k,l使得vk∈dom(xi),vl∈dom(xj)；現定義1個五元組edge=〈xi,vk,xj,vl,s〉，其中s為2個點(xi,vk)與(xj,vl)對于c的滿足情況.由此可以將edge的概念推廣到標準約束網絡.

在標準約束網絡P=(X,C)中，?c∈C，有?i,…,k,使得scp(c)={xi,…,xk}，?l,…,n使得vl∈dom(xi),…,vn∈dom(xk)；現定義1個元組edge=〈xi,vl,…,xk,vn,s〉，其中s為一組點(xi,vl),…,(xk,vn)對于c的滿足情況.

定義9. 邊集(edges).edges是P中所有edge元組組成的集合.

edge∈edges∧edge.s=⊥}.

xk,vn}∈edge,edge∈edges}.

故edge是邊集E中對應元組擴展1個標識位s后得到的，這個標識位用來表示該元組的相容情況；node是點集中V對應元組擴展1個標識位e后得到的，這個標識位表示該點是否存在.

對于如例1所示的二元約束網絡，用N-E表示的形式如表1、表2所示:

Table 1 Data in nodes

Table 2 Data in edges

原有的弧相容算法大多基于串行思想提出，具有一般串行算法的特點：1)算法流程中含有大量復雜的邏輯判斷，具有分支多、邏輯復雜的特點；2)對附加數據結構的處理大多采用串行的方式，對于具有較大的附加數據結構的算法，往往以串行的形式構建它們時就花費了較多的時間.而GPU是細粒度的并行計算設備，在其上執行這種分支較多，具有復雜邏輯判斷的串行算法非常耗時.然而GPU可以快速地以并行的形式創建大規模的附加數據結構，并對其進行并行操作，這大幅節約了算法的執行時間.

N-E模型是一種對約束網絡的細粒度表示方法，它將約束網絡每個點的存在情況記錄在nodes中，將每組值對的相容性記錄在edges中.這種模型內部的元素具有相同的結構，元素之間獨立性強，便于在GPU上生成并進行大規模的并行計算處理.

定理1. 對于給定的標準約束網絡P，其N-E模型表示形式為N_E(P)=(edges,nodes).其中nodes的空間復雜度為O(nd)；每個約束含有的相容和不相容的值對(即edge)為O(dr)，則edges的空間復雜度為O(edr)，對于標準二元約束網絡edges的空間復雜度為O(ed2).

3.2 基于GPU生成N-E模型的并行算法：BuildNE

盡管N-E模型具有較高的復雜度，但得益于GPU在并發處理集合的能力，將約束網絡P轉換為N-E模型的過程較為便捷.本節我們提出生成N-E模型的算法BuildNE算法.

算法2. BuildNE算法.

輸入：P=(X,C)；

輸出：N_E(P).

①ComputeEdgeOffset〈|C|〉 (in:P; out:

ole)；

②oe←Exclusive_Scan(ole)；

③ComputeNodeOffset〈|X|〉 (in:P; out:

oln)；

④on←Exclusive_Scan(oln)；

⑤BuildEdgesLaunch〈|C|〉(in:P,oe; out:

edges)；

⑥BuildNodesLaunch〈|X|〉(in:P,on; out:

nodes).

過程1.ComputeEdgeOffset.

輸入：P；

輸出：ole.

①c←C[thread_id]；

②x0,x1←scp(c)；

③ole[thread_id]←|dom(x0)|×|dom(x1)|.

過程2.ComputeNodeOffset.

輸入：P；

輸出：oln.

①x←X[thread_id]；

②oln[thread_id]←|dom(x)|.

過程3.BuildEdgesLaunch.

輸入：x0,x1,c,start_index；

輸出：edges.

①c←C[thread_id]；

②x0,x1←scp(c)；

③start_index←oe[thread_id]；

④BuildEdges〈|dom(x0)|×|dom(x1)|〉(in:

x0,x1,c,start_index; out:edges).

過程4.BuildEdges.

輸入：P,on；

輸出：edges.

①v0←x0[block_id]；

②v1←x1[thread_id]；

③index←start_index+thread_id+

|dom(x0)|×block_id；

④edges[index]←〈x0,v0,x1,v1,s〉.

過程5.BuildNodesLaunch.

輸入：P,on；

輸出：nodes.

①x←X[thread_id]；

②start_index←on[thread_id]；

③BuildNodes〈|dom(x)|〉 (in:x,

start_index; out:nodes).

過程6.BuildNodes.

輸入：x,start_index；

輸出：nodes.

①v←x[thread_id]；e←1；

②index←start_index+thread_id；

③nodes[index]←〈x,v,e〉.

偽代碼函數“F〈n〉”尖括號括起來的部分表示啟動多個線程執行相同的函數F，n表示啟動的線程數目.

算法2的行①為每個約束啟動1個線程，計算該約束含有edge元組個數存入ole數組中.行②對ole數組進行排他性掃描求oe數組，oe表示對應約束的元組1在edges數組的位置.行⑤啟動過程3，以oe為參數，先啟動|C|個線程，每個線程可以得到1個約束c∈C,scp(c)={x0,x1}，在每個線程中再次啟動|dom(x0)|×|dom(x1)|個線程執行過程4，過程4每個線程從程序運行環境取得v0和v1，從參數中取出x0和x1，并查表取得相容性情況s，填滿edge的5個元組，最后按照oe經計算確定edge元組在edges數組的位置.

算法2的行③④⑥用來計算nodes，所用方式與計算edges的方式類似，這里不予贅述.

4 并行弧相容算法

4.1 AC4GPU

算法2得到了標準二元約束網絡P的N-E模型：N_E(P)=(nodes,edges)，算法3在此基礎上提出基于GPU執行的并行弧相容算法——AC4GPU.

算法3. AC4GPU.

輸入：N_E(P)；

輸出：P滿足AC輸出真，反之輸出假.

loop

①counter[xi,vi,xj,s]←以edges中(x0,v0,x1)和(x1,v1,x0)為鍵歸約edges.s；

② foreach(xi,vi)∈counters.t.counter[xi,vi,xj]=0 in parallel

③nodes[xi,vi].e←0；

④ end foreach

⑤ 同步stream0,stream1；

⑥stream0: foreachedge∈edgesin parallel

⑦ ifnodes[x,v].e=0 then

⑧edges[xi,vi].s←0 andedges[xj,vj].s←0；

⑨ end if

⑩ end foreach

go to loop.

本算法出現的stream語句采用了CUDA流的技術，并發執行算法中可以并行的語句，具有相同名稱的stream語句串行執行，不同名稱的stream之間并行執行，同步語句同步不同的stream.

1) 若P是相容的，因為算法在最后1次循環沒有刪點，那么在最后1次循環計算出的domX數組的每個元素與上次記錄的元素都相同，即domX=domX_pre(第1次循環的domX_pre為各變量初始論域的大小)，這時算法返回true.

2) 若P是不相容的，那么此時domX數組中至少有1個元素為0，即?x∈X，使得domX[x]=0，這時算法返回false.算法啟動后，會一直循環地執行下去，直到P達到滿足上述兩者之一時退出，并給出P的相容性狀況及用N-E表示的約束網絡.

如例1所示的約束網絡上執行算法3的過程如表3～14所示：

Table 3 counter: Obtained from Reduce Table 2 by 3 Keys in Line ①

Table 4 nodes: Modified According to Table 3 in Line ②

Table 5 edges: Modified According to Table 4 in Line ⑥, Parallel with Table 6

Table 6 domX: Modified According to Table 4 in Line , Parallel with Table 5

Table 7 counter: Modified According to Table 5 in 2nd Loop Line ①

Table 8 nodes: Modified According to Table 7 in 2nd Loop Line ②

Table 9 edges: Modified According to Table 8 in 2nd Loop Line ⑥, Parallel with Table 10

Table 10 domX: Modified According to Table 8 in 2nd Loop Line , Parallel with Table 9

Compared with Table 6, some items have been modified. It does not meet the exit condition. The loop will continue.

Table 11 counter Modified According to Table 9 in 3rd Loop Line ①

Table 12 nodes: Modified According to Table 11 in 3rd Loop Line ②

Table 13 edges: Modified According to Table 12 in 3rd Loop Line ⑥, Parallel with Table 14

Table 14 domX: Modified According to Table 12 in 3rd Loop Line , Parallel with Table 13

Compared with Table 10, all items indomXhave not been modified. It means the networkPis AC and the program exits.

定理2. 對于給定的標準二元約束網絡P，對應的N-E模型表示為N_E(P)=(nodes,edges)，在 該 模型上執行AC4GPU算法的步驟復雜度為O(ndloged2)，空間復雜度為O(ed).

該算法附加數據結構為counter，domX與domX_pre，其中counter的空間復雜度O(ed)，domX與domX_pre復雜度均為O(n)，故整個算法的空間復雜度為O(ed).

證畢.

4.2 AC4GPU+

從第5節實驗數據我們可以看出，AC4GPU算法得益于GPU執行大規模并發算法的優秀表現，在單次弧相容算法中取得了較好的效果.但通過分析算法可以發現，每次循環都要執行的按鍵歸約使算法性能大打折扣.以算法3行①為例，該行代碼意在N-E模型上求1次counter數組，相當于原AC4算法中的初始化操作，然而AC4算法只執行1次初始化操作，相比較而言AC4GPU算法則需在每次循環時都執行1次，這大幅增加了執行時間.

AC4GPU+將AC4GPU拆分成初始化和傳播2個階段，通過對N-E模型的擴展，為每個edge元組添加字段m0,m1，用來記錄該edge歸約到counter數組的位置，記為N_E+(P)=(nodes,edges+).算法引入一種特殊的鎖機制——原子操作[16]，通過m0,m1字段，直接計算counter值，從而將原算法循環最耗時的按鍵歸約移到循環之外構成初始化階段；傳播階段則采用循環的方式每次刪掉當前中網絡中所有不相容的點，該階段求domX的按鍵歸約改為復雜度為O(1)的原子自減操作，在網絡相容性判斷的主要判斷分支采用O(1)的判斷語句，最終使得整個傳播過程中絕大多數循環的復雜度為O(1).

算法4. AC4GPU+.

輸入：N_E+(P)；

輸出：P滿足AC輸出真，反之輸出假.

①counter[xi,vi,xj,s]← 以(x0,v0,x1)和(x1,v1,x0)為鍵歸約edges.s；

② 初始化domX；

③ if {ctr|ctr∈counter[xi,vi,xj]=0}≠?

then；

④ return true；

⑤ end if

⑥ foreachcounter[xi,vi,xj]=0 in parallel

⑦nodes[xi,vi].e=0；

⑧domX[xi]←atomic_add(domX[xi],-1)；

⑨ ifdomX[xi]=0 then

⑩ return false；

loop

parallel

gotoloop.

算法4的初始化階段與算法3的第1次循環相似，對edges按鍵歸約求counter，修改node的e值，并且驗證整個網絡相容性.

Fig. 3 A special-case of delete node圖3 刪點的特殊情況

圖3中的特殊情況是scp(ci j)={xi,xj},scp(cj k)={xj,xk},vi=dom(xi),vj∈dom(xj),vk=dom(xk)，現在刪除(xi,vi)和(xk,vk)2個點，且這2個點是(xj,vj)分別在ci j與cj k的最后支持，此時counter[xj,vj,xi]=0且counter[xj,vj,xk]=0，在用原子操作計算domX(xj)會自減2次，這意味著(xj,vj)被刪除了2次，與實際情況不相符，說明在算法執行的過程中domX的值可能是不對的.由于nodes保存的每個點存在的信息是正確的，這時需要對nodes進行1次按鍵歸約，重新修改domX數組的值，如果這個時候發現domX[xi]仍然為0，則表示沒有出現如上這類特殊情況，變量xi的論域確實已被刪空，網絡不相容循環退出.需要注意的是上述這類情況，同樣適用于同時刪除多個點的情況.

使用這種驗證方式是因為在傳播過程的每次循環中，出現這類特殊的刪點情況且將論域刪空的情況并不多見，所以只需要在出現論域刪空的情況下，再次以復雜度較高的方式驗證，而其他時間，則以O(1)的方式進行驗證，這樣達到了節省時間的目的.

同AC4算法相比， AC4GPU+算法初始化階段求counter的過程是通過對edges數組按鍵歸約完成的；整個edges數組由于記錄了約束網絡的關系，可以看作是AC4算法中的S數組；nodes數組則可以看成是AC4算法的待刪除點的隊列Q.

同AC4GPU算法相比，AC4GPU+算法將原算法單一的循環，拆分成初始化和傳播2個階段，整個算法僅在初始化時執行1次最為耗時的按鍵歸約操作，在傳播階段每次循環都刪除所有不相容的點，而在該階段對網絡的相容性驗證，絕大多數情況采用復雜度為O(1)的判斷分支，僅在出現變量值域刪空的情況下，采用步驟復雜度較高的分支，驗證是否產生特殊情況.

定理3. 對于給定的標準二元約束網絡P，對應的N-E模型表示為N_E(P)=(nodes,edges)，在 該 模型上執行AC4GPU+算法的步驟復雜度為O(ndlognd)、空間復雜度為O(ed).

綜上，2個階段步驟復雜度的和為O(loged2+ndlognd)，在標準二元約束網絡中n≥2,d≥1,1≤e≤n2，即loged2≤logn2d2=2×lognd≤nd×lognd，即loged2≤nd×lognd，則理論上整個算法的步驟復雜度為O(ndlognd)，該算法較之AC4GPU在理論上有顯著提升.

算法4附加數據結構為counter，domX，如4.1節所證，counter的空間復雜度O(ed)，domX空間復雜度為O(n)，故整個算法的空間復雜度為O(ed).

證畢.

需要注意的是：我們理論上按照最壞的情況計算了算法復雜度，通過比較忽略了初始化部分的步驟復雜度，但并不表示在所有測試用例中初始化的時間都可以忽略.這里ndlognd的系數nd指的是刪除點所用的循環次數，循環次數的取值范圍是[0,nd]，若傳播階段沒有需要被刪除的點，則該階段耗費時間為0.所以本算法在具體測試用例的執行中，若循環次數小于特定值，則可能出現初始化階段所耗費的時間大于傳播階段，這與循環次數以及循環內部執行的具體的判斷分支有關，因問題而異，這里不能給出準確閾值.

5 實驗數據

5.1 實驗程序設計

由于AC4GPU和AC4GPU+是基于AC4的并行化改進，我們將3個算法進行對比實驗，為了獲得實驗流程各個細節的數據，設計了如下的實驗程序：

如圖4所示，對比程序將benchmark文件讀入內存；再將數據存入顯存中；分別在CPU中建立約束網絡，在GPU中建立N-E模型；最后分別在CPU執行串行的AC4算法，在GPU中執行并行的AC4GPU和AC4GPU+算法.實驗數據記錄各模塊的執行時間，下面會對這些數據進行分析.

本文采用目前國際通用的測試約束滿足問題的標準庫用例，所有測試文件可以從文獻[17]在線下載.

實驗環境：CPU為Intel Core i5-4590 3.33 GHz， 4 GB RAM；GPU為NVIDIA GTX750[18]，512個CUDA核心，GPU Clock 1215 MHz，1 GB RAM.

Fig. 4 The experimental process design圖4 實驗流程設計

5.2 N-E模型測試實驗

為了測試各用例在建立N-E模型的情況，表15給出了常見測試用例的文件大小、傳輸時間、建模時間、edge個數、node個數.

Table 15 Size and Build Model Time of Benchmarks

從實驗數據可以看出，測試用例在內存到顯存的傳輸時間不可忽略，主要由于現有編程模型對GPU的限制，GPU不能直接從文件中讀取數據到顯存，必須通過 CPU先將數據讀入內存，再將數據從內存復制到顯存中，內存的吞吐量遠不及核心計算的吞吐量，受限于多核系統的結構，內存數據傳輸延遲成為了制約并行算法加速的一大難題[19].與之相比，建模時間較短且與測試用例的問題規模有關，即與N-E模型的大小有關，然而如果接下來繼續在GPU上執行弧相容算法或求解算法，則數據拷貝僅此1次且是單向的，此后數據結構被維持在顯存中.

在GPU中對約束網絡進行N-E建模后，GPU中執行并行弧相容算法AC4GPU和AC4GPU+，它們與AC4算法執行時間的比較如表16、表17所示.

Table 16 Experimental Comparison of Composed Problems

Composed問題(表16中簡寫為cpd)是一組具有特殊結構的問題，由于許多實例在執行第1次弧相容時就需要刪點，這里用來測試并行弧相容算法，表16中給出了部分Composed用例的算法性能對比情況，其中循環次數是指在執行算法3、算法4循環的次數；每次循環中刪去點的個數用“”隔開，例如611表示在算法3第1次循環(對應算法4初始化階段)刪去6個點，而后每次循環刪去1個點.從實驗數據可以看出，AC4算法的執行時間隨問題規模的增大而增大，主要原因是由于初始化占用了大量時間.而AC4GPU和AC4GPU+的執行時間與問題規模的關系不明顯，卻與循環次數有關.

通過cpd-25-1-25-1 和cpd-75-1-2-4執行時間進行對比我們可以看出，其AC4GPU算法的執行時間與用例本身的迭代次數成正比，即執行算法的循環迭代的越多，耗費的時間越多，AC4GPU每次增加一次迭代大約增加1ms的執行時間；AC4GPU+對循環進行了改進，執行時間明顯減少.在刪除節點較多的情況下，由于AC4GPU和AC4GPU+算法是在每次循環中并行刪除多個點，而AC4算法則是在每次循環刪除1個點，但由于并行算法的內存讀取速度限制，其性能有否提升還有待研究.

表17給出了3個算法在一些隨機問題上的性能比較，由于這些問題在執行弧相容算法時不刪除點，所以循環次數為1.由實驗數據可以看出，串行算法執行時間隨問題規模的增大而增大；并行算法的執行時間取決于GPU的并行度.若測試用例的問題規模測大于GPU可啟動的線程數，GPU需要分多次啟動線程進行計算，由此可知，對于那些GPU啟動相同次數進行處理的測試用例，它們的執行時間相近，故并行算法的執行時間隨啟動次數增加而增加.

Table 17 The Experimental Comparison of Random Problem

本實驗采用的GPU為GTX750，該GPU具有4個流多處理器(SM)，每個SM上有128個流處理器(SP)，共有512個CUDA核心(SP)，每個SM最多可以啟動2 048個線程，總共可以啟動8 192個線程，如果輸入數組長度過大，則會將原數組分割成多個部分進行計算.由于edge具有對稱性，即〈x0,v0,x1,v1,s〉等價于〈x1,v1,x0,v0,s〉，本實驗程序采用BuildNE建立1個edge時同時生成了2個元組，即對于問題tightness0.8，實際上其edges擁有1 318 400個元素，底層歸約算法將原數組分割成了多個部分計算，故并行算法的執行時間相對較長.同理，對于問題tightness0.9，底層算法將其分割為更多部分進行計算，所以耗時更長.

綜上，并行弧相容算法用并行的手段節約了AC4在初始化階段的巨大花費，并在每次循環中并行地刪除所有不相容的點，最終經實驗結果證明，該算法較AC4算法有所提高.

6 結論與展望

針對基于GPU的并行計算在解決某些問題的優點，本文對結合并行計算的基本算法對弧相容問題進行了研究，提出了一種便于GPU并行計算約束網絡的表示模型——N-E模型；給出了在GPU中生成該模型的算法——BuildNE；并將經典的AC4算法結合并行計算的基本思想，提出一種細粒度的并行弧相容算法——AC4GPU及其改進AC4GPU+.進而通過實驗給出了常見測試用例以N-E模型表示的大小與的創建時間，分析并行弧相容算法較于傳統弧相容算法的特點和優勢.

利用算法并行刪除多個點的優點，對于一些需要快速刪點的求解問題，并行算法一定會發揮很大作用.未來的工作是在此基礎上對并行化約束滿足問題進行進一步的研究，包括對并行化分支搜索、單弧相容算法、約束求解算法等基于弧相容算法的高級算法進行并行化的擴展；利用N-E數組元素間的相對獨立性，研究該模型在約束劃分、辨別特殊的約束網絡的模式上的可能應用.

[1]Freuder E C, Mackworth A K. Constraint Satisfaction: An Emerging Paradigm[M]. Amsterdam: Elsevier, 2006: 13-27

[2]Bessière C. Constraint Propagation[M]. Amsterdam: Elsevier, 2006: 29-84

[3]Mohr R, Henderson T C. Arc and path consistency revisited[J]. Artificial Intelligence, 1986, 28(2): 225-233

[4]Lecoutre C. Generic GAC Algorithms[M]. London: ISTE Ltd, 2009:185-236

[5]Nguyen T, Deville Y. A distributed arc-consistency algorithm[J]. Science of Computer Programming, 1970, 30(12): 227-250

[6]Yokoo M, Durfee E H, Ishida T, et al. The distributed constraint satisfaction problem: Formalization and algorithms[J]. IEEE Trans on Knowledge & Data Engineering, 1998, 10(5): 673-685

[7]Ringwelski G. An arc-consistency algorithm for dynamic and distributed constraint satisfaction problems[J]. Artificial Intelligence Review, 2005, 24(34): 431-454

[8]Campeotto F, Palù A D, Dovier A, et al. Exploring the use of GPUs in constraint solving[C]Practical Aspects of Declarative Languages. Berlin: Springer, 2014: 152-167

[9]Gent I P, Jefferson C, Miguel I, et al. A preliminary review of literature on parallel constraint solving[C]Proc of PMCS 2011 Workshop on Parallel Methods for Constraint Solving. Berlin: Springer, 2011: 499-504

[10]Christophe L. Constraint Networks[M]. London: ISTE Ltd, 2009: 39-92

[11]Wilt N. Reduction[G]CUDA Handbook: A Compr-ehensive Guide to GPU Programming. Reading, MA: Addison-Wesley, 2013: 365-384

[12]NVIDIA. CUDA C Programming Guide, version 7.0[OL]. [2015-06-25]. https:developer.nvidia.comcuda-zone

[13]Wilt N. Scan[G]CUDA Handbook: A Comprehensive Guide to GPU Programming. Reading, MA: Addison-Wesley, 2013: 173-204

[14]Nicholas W. Stream and event[G]CUDA Handbook: A Comprehensive Guide to GPU Programming. Reading, MA: Addison-Wesley, 2013: 385-420

[15]NVIDIA. Thrust Quick Start Guide. version 7.0[OL]. [2015-06-25]. http:docs.nvidia.comcudapdfThrust_Quick_Start_Guide.pdf

[16]Butcher P. Seven Concurrency Models in Seven Weeks: When Threads Unravel[M]. Raleigh, NC: Pragmatic Bookshelf, 2014

[17]Lecoutre C. Benchmarks of constraint satisfaction problem[OL]. [2012-04-16]. http:www.cril.univ-artois.fr～lecoutrebenchmarks.html

[18]NVIDIA. Geforce GTX 750 specifications[OL]. [2015-06-25]. http:www.geforce.cnhardwaredesktop-gpusgeforce-gtx-750

[19]Sun Xianhe, Chen Yong. Reevaluating Amdahl’s law in the multicore era[J]. Journal of Parallel & Distributed Computing, 2010, 70(2): 183-188

Li Zhe, born in 1990. MSc. His main research interest is constraint programming (leezear@live.cn).

Li Zhanshan, born in 1966. PhD, professor of Jilin University. Member of CCF. His main research interests include constraint solving, model-based diagnosis, intelligent planning and scheduling.

Li Ying, born in 1965. PhD, professor, PhD supervisor. Her main research interests include large data processing, geological 3D visualization modeling, 3D geological image processing, geological model data processing, etc (l_ying@jlu.edu.cn).

A Constraint Network Model and Parallel Arc Consistency Algorithms Based on GPU

Li Zhe, Li Zhanshan, and Li Ying

(CollegeofComputerScienceandTechnology,JilinUniversity,Changchun130012) (KeyLaboratoryofSymbolicComputationandKnowledgeEngineering(JilinUniversity),MinistryofEducation,Changchun130012)

Constraint satisfaction problem is a popular paradigm to deal with combinatorial problems in artificial intelligence. Arc consistency (AC) is one of basic solution compression algorithms of constraint satisfaction problem, which is also a core algorithm of many excellent advanced algorithms. When constraints are considered independently, AC corresponds to the strongest form of local reasoning. An effective underlying AC can improve the efficiency of reducing the search space. Recent years, GPU has been used for constituting many super computers, which solve many problems in parallel. Based on GPU’s computation, this paper proposes a constraint networks presentation model N-E and its parallel generation algorithm BuildNE. According to fine-grained arc consistency AC4, a parallel edition AC4GPUand its improved algorithm—AC4GPU+, are proposed. The two parallel algorithms extend arc consistency to GPU. Experimental results verify the feasibility of these new algorithms. Compared with AC4, the parallel versions have made the 10% to 50% acceleration in some smaller instances, and obtained 1 to 2 orders of magnitude in some bigger instances. They provide a core algorithm to other constraint satisfaction problem solving in parallel for further study.

artificial intelligence; constraint satisfaction problem (CSP); arc consistency (AC); GPU; compute unified device architecture (CUDA)

2015-10-13；

2016-08-08

國家自然科學基金項目(61272208,61373052)；吉林省自然科學基金項目(20140101200JC) This work was supported by the National Natural Science Foundation of China (61272208, 61373052) and the Natural Science Foundation of Jilin Province of China (20140101200JC).

李占山(zslizsli@163.com)

TP18