• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于國產(chǎn)異構(gòu)計算平臺的快速SVD算法及其在海洋資料同化的應(yīng)用

    2024-02-27 09:53:38李維釗王偉
    關(guān)鍵詞:數(shù)據(jù)結(jié)構(gòu)進程資料

    李維釗,王偉

    易海陸圓(山東)數(shù)字技術(shù)有限公司,山東 青島 266237

    引 言

    數(shù)值預報屬于微分方程的初值問題,初始場的正確性對于預報結(jié)果具有決定性作用。利用數(shù)值模式對大氣、海洋等進行數(shù)值預報,已經(jīng)成為大氣、海洋環(huán)境預測預報的主要方式。數(shù)值預報的資料同化,實質(zhì)上是利用不同來源的觀測資料不斷修正模式的預報,從而迫使數(shù)值預報產(chǎn)生更加接近實況的模式環(huán)境動力狀態(tài)的過程,它為下一時刻的數(shù)值預報提供更加準確的初始場,從而提高數(shù)值預報的準確性。無論在理論研究還是工程實踐中,資料同化已經(jīng)成為數(shù)值預報的關(guān)鍵研究內(nèi)容。

    變分法和濾波法是現(xiàn)代資料同化中最有代表性的技術(shù)方法。1959年Cressman在該方法的基礎(chǔ)上采用迭代求解方法形成了實際可以操作的逐步迭代法;1963 年Gardin 提出了最優(yōu)插值法,建立了資料同化的理論基礎(chǔ)——統(tǒng)計估計理論;20世紀60年代后期,隨著觀測資料的時空范圍和種類的增加,出現(xiàn)了變分法的資料同化方法,1969年Thompson等更是提出一套較為完善的四維變分理論和工程方法,但是對于強非線性問題性能較差;1960 年,數(shù)學家Kalman 提出了卡爾曼濾波方法;1994年Evensen等首次提出的集合卡爾曼濾波方法,可以應(yīng)用于強非線性系統(tǒng),經(jīng)過大量的理論研究和應(yīng)用實踐,證明集合卡爾曼濾波以及由此發(fā)展出來的集合最優(yōu)插值具有很好的數(shù)據(jù)資料同化效果和不同場景的適應(yīng)能力,得到廣泛的認可。

    集合卡爾曼濾波方法解決了誤差協(xié)方差矩陣預報中遇到的問題,在實際應(yīng)用中,由于大量不同來源、不同應(yīng)用目標、不同空間解析度和對不同物理分量預測預報的需求,導致資料同化的數(shù)據(jù)量和計算量巨大,算法的計算效率成為制約科研和應(yīng)用的一個重要因素?,F(xiàn)代資料同化算法基本上是通過數(shù)值計算的方式迭代求解,特別是矩陣求逆算法,在變分法和集合卡爾曼濾波算法中,成為資料同化的關(guān)鍵技術(shù)之一,也是計算量的主要來源。本文從以下兩個方面提高海量資料同化的效率問題:

    (1)面向大規(guī)模計算環(huán)境的并行策略及實現(xiàn)方式;

    (2)面向異構(gòu)計算平臺的使用奇異值分解(SVD)的矩陣求逆批量快速算法。

    本研究使用C/C++完整地實現(xiàn)了包括SVD在內(nèi)的海洋資料同化軟件。資料同化的實現(xiàn)程序大多采用FORTRAN語言實現(xiàn),而某些國產(chǎn)自主的計算平臺在支持FORTRAN 語言的生態(tài)方面存在不足,與此形成鮮明對照的是工程人員選擇更具備生態(tài)優(yōu)勢的C/C++語言實現(xiàn)程序功能。本研究同時也證明了,只要合理設(shè)計,借助于現(xiàn)代編譯系統(tǒng)的優(yōu)化功能,在同樣的應(yīng)用場景下,C/C++可以達到FORTRAN語言相同的計算效率。

    1 基于集合卡爾曼濾波的資料同化及實現(xiàn)流程

    基于文獻[1]的集合卡爾曼濾波及集合最優(yōu)插值,國內(nèi)外出現(xiàn)大量的基于FORTRAN語言和線性代數(shù)庫的算法在高性能計算環(huán)境或計算卡下的實現(xiàn),并進行了大量的優(yōu)化工作,具體見文獻[2-6],盡管如此,資料同化算法仍然是非常耗時的,尤其是對于國產(chǎn)異構(gòu)計算平臺,提升基礎(chǔ)算法的計算效率是應(yīng)用推廣的關(guān)鍵因素之一。

    參考文獻[2],集合最優(yōu)插值算法的表述如下:

    其中:φa是分析場數(shù)據(jù),φ是背景場數(shù)據(jù),α是區(qū)間( 0,1] 的隨機因子;A是靜態(tài)樣本數(shù)據(jù),A′是擾動預測矩陣,由A中心化而來,A′A′T是背景誤差協(xié)方差的估計;H是測量算子,代表真實模型與預測模型之間的映射;γ是擾動測量誤差向量,γγT是觀測誤差協(xié)方差;d′=d-Hφ是觀測增益,其中d是擾動測量向量。經(jīng)過一系列的推導,參考文獻[2],可以得到如下的公式:

    其中U、Λ 滿足Λ=ΣΣT,對進行SVD分解,得到U向量和奇異值Σ,結(jié)合矩陣的性質(zhì)等,得到如下步驟計算:

    靜態(tài)樣本的選取參考文獻[7],采用文獻[8]中提到的局部集合最優(yōu)插值方法進行資料同化。本文基于國產(chǎn)異構(gòu)計算平臺,采用CPU和類GPU卡協(xié)同方式實現(xiàn)海量資料同化的并行計算。

    2 資料同化工程實現(xiàn)的并行策略

    2.1 資料同化計算任務(wù)分解

    以資料同化為例,一般的大規(guī)模計算應(yīng)用需要完成如下幾項工作:

    2.1.1 數(shù)據(jù)讀取

    從存儲系統(tǒng)將數(shù)據(jù)讀入內(nèi)存。主要的數(shù)據(jù)包括靜態(tài)樣本數(shù)據(jù)、背景場數(shù)據(jù)和觀測場數(shù)據(jù),本研究中靜態(tài)樣本數(shù)據(jù)量在數(shù)百GB,通過數(shù)個IO代理進程從系統(tǒng)的存儲設(shè)備中將數(shù)據(jù)讀入CPU內(nèi)存,再將數(shù)據(jù)分發(fā)給資料同化的計算進程。數(shù)據(jù)在文件中的數(shù)據(jù)結(jié)構(gòu)對IO的性能有顯著的影響,原則上,應(yīng)該優(yōu)化數(shù)據(jù)在文件中的存儲方式,使得打開/關(guān)閉文件的操作數(shù)量盡量少,并且文件操作的并發(fā)數(shù)量要與存儲系統(tǒng)的配置相匹配。

    2.1.2 數(shù)據(jù)分發(fā)

    采用并行計算的方式,資料同化軟件需要將靜態(tài)樣本、背景場和觀測場的數(shù)據(jù)通過通信網(wǎng)絡(luò)分發(fā)到計算進程。

    2.1.3 資料同化計算

    資料同化計算部分主要包括三個部分:SVD之前的數(shù)據(jù)插值等前處理;SVD計算任務(wù);SVD 后處理。SVD 前、后的處理主要為矩陣的乘法運算。同時,SVD 的主要計算量來源是矩陣的二對角化。單純用CPU 或計算卡進行SVD,非常耗時,因此,在異構(gòu)計算平臺上,發(fā)揮CPU 處理存在復雜邏輯的計算,結(jié)合計算卡的數(shù)值計算性能,實現(xiàn)高效SVD的工程實現(xiàn)算法,是提高資料同化的計算效率的關(guān)鍵技術(shù)。SVD算法在量子計算、人工智能、機器視覺、信號降噪等領(lǐng)域應(yīng)用廣泛,實現(xiàn)高效SVD 的工程實現(xiàn)算法具有重要的實用價值。

    2.1.4 數(shù)據(jù)回寫

    計算進程計算完成后,需要將結(jié)果數(shù)據(jù)回寫到分析場數(shù)據(jù)文件中。

    已有大量的文獻對資料同化工程實現(xiàn)的并行策略進行了研究,并取得不同程度的性能提升。與文獻[3]和文獻[2]相似,本研究資料同化的并行策略可以總結(jié)為計算任務(wù)和流水線并行策略,其核心為任務(wù)的編排。

    2.2 任務(wù)編排

    計算任務(wù)的并行度,取決于任務(wù)粒度、通信方式和負載平衡三個方面因素。因此,資料同化任務(wù)的分解和調(diào)度是選擇并行策略的基礎(chǔ)。

    海洋資料同化的一般流程如圖1所示。

    圖1 海洋資料同化流程圖Fig.1 Ocean data assimilation flowchart

    按照海洋資料同化流程,結(jié)合國產(chǎn)異構(gòu)計算平臺的環(huán)境配置,實現(xiàn)并行化設(shè)計,海洋資料同化的任務(wù)模塊劃分如下:

    2.2.1 數(shù)據(jù)訪問任務(wù)模塊

    同化軟件啟動Nio個IO 代理進程實現(xiàn)海洋資料同化的數(shù)據(jù)文件的輸入和輸出。數(shù)據(jù)訪問的效率取決于存儲系統(tǒng)的并行特性和文件的數(shù)據(jù)結(jié)構(gòu)。其中文件的數(shù)據(jù)結(jié)構(gòu)是關(guān)鍵。在HYCOM 數(shù)值模式生成的數(shù)據(jù)中,共有120 個包括垂向36層網(wǎng)格數(shù)據(jù)的文件。資料同化程序的代理進程讀取相應(yīng)的數(shù)據(jù),并將數(shù)據(jù)按照計算進程的數(shù)量及序號,平均分配到相應(yīng)的計算進程執(zhí)行同化計算。原則上,120份垂直36層的靜態(tài)樣本數(shù)據(jù)可以組織到1或多個文件中,工程實現(xiàn)時,文件的數(shù)量NF和每個文件中的數(shù)據(jù)結(jié)構(gòu)對數(shù)據(jù)IO 的效率影響明顯。以溫度樣本為例,原程序文件組織形式為NF=120 個文件,每個文件中的溫度的數(shù)據(jù)結(jié)構(gòu)如圖2所示。

    圖2 靜態(tài)樣本數(shù)據(jù)原始組織方式示意圖Fig.2 Schematic diagram of the original organization of static sample data

    由于每個計算進程分配的數(shù)據(jù)為36層數(shù)據(jù)中相同區(qū)域的網(wǎng)格點數(shù)據(jù),如圖3所示(Np為計算進程的數(shù)量)。

    圖3 靜態(tài)樣本數(shù)據(jù)在計算進程中的分配示意圖Fig.3 Schematic diagram of static sample data allocation in the calculation process

    其中,總的文件數(shù)量為NF=120,每個文件中都有圖3 所示的一份數(shù)據(jù)。每層數(shù)據(jù)的讀取均遵循“打開文件—讀取數(shù)據(jù)—關(guān)閉文件”的文件讀取操作,每個進程執(zhí)行文件讀取操作的數(shù)量K=NF×36,文件讀取操作的總數(shù)量為KTotal=Nio×NF×36,當Nio=8,NF=120 時,KTotal=34,560 次文件操作。由于存在多個代理進程讀取數(shù)據(jù),120 個文件的讀寫存在大量的沖突,超出存儲系統(tǒng)的并發(fā)性能,造成文件讀取效率低下。數(shù)據(jù)IO 優(yōu)化的目標有兩個:一是通過優(yōu)化數(shù)據(jù)結(jié)構(gòu),減少文件操作的次數(shù),增加單次文件連續(xù)讀操作的數(shù)據(jù)量;二是將文件合理放置在存儲系統(tǒng)不同物理存儲單元上,實現(xiàn)數(shù)據(jù)訪問的高并發(fā)性,如文獻[2],在基于Lustre 的文件系統(tǒng)上,通過將文件放置在不同OST上,取得明顯的性能優(yōu)化。本文考慮到存儲文件系統(tǒng)的多樣性和復雜性,通過優(yōu)化文件的數(shù)據(jù)結(jié)構(gòu)實現(xiàn)IO性能的提升。根據(jù)局部集合最優(yōu)插值算法,網(wǎng)格點的同化計算順序不影響最終結(jié)果,因此,設(shè)定1個最大數(shù)量的IO進程數(shù)Nio-max=32,在資料同化進行之前,將靜態(tài)樣本數(shù)據(jù)的文件存儲格式按照如圖4組織(文件數(shù)量NF=120)。

    圖4 靜態(tài)樣本數(shù)據(jù)優(yōu)化組織方式示意圖Fig.4 Schematic diagram of static sample data optimization organization method

    由圖4可知,數(shù)據(jù)優(yōu)化的方式就是將原來按照層序號連續(xù)放置的一維方式,改變?yōu)榘凑誟層、進程]的二維方式,適應(yīng)計算進程對數(shù)據(jù)結(jié)構(gòu)要求,遵循同樣的文件操作方式,文件讀操作的最大數(shù)量變?yōu)椋∟io=8 時)KTotal=Nio×NF=8×120=960,不同的文件代理進程讀取不同的文件,結(jié)合數(shù)據(jù)分發(fā)任務(wù)模塊的具體特點,避免了進程間文件訪問的沖突,實驗證明,優(yōu)化措施帶來顯著的文件性能IO 提升,并且具有一定的普遍適用性。表1 是Nio=4,8,12,Nf=120 時,兩種讀寫方式的運行時間對比。優(yōu)化后運行效率提高到原來的10倍(Nio=8 時性能最優(yōu),不同存儲系統(tǒng)存在差異)。

    表1 不同讀寫策略讀寫效率對比Table 1 Comparison of reading and writing efficiency of different reading and writing strategies

    2.2.2 數(shù)據(jù)分發(fā)任務(wù)模塊

    通過文件讀取操作的優(yōu)化,文件中的數(shù)據(jù)結(jié)構(gòu)得到優(yōu)化,IO代理進程在讀取數(shù)據(jù)后,不需要進一步組織數(shù)據(jù),就可以把讀取的數(shù)據(jù)直接分發(fā)給計算進程,計算進程接收到數(shù)據(jù)后,啟動執(zhí)行同化計算的任務(wù)模塊,按照網(wǎng)格點實施同化計算,文獻[2]中對采用MPI 相關(guān)通信機制進行了詳細說明。研究證明,通過對文件數(shù)據(jù)結(jié)構(gòu)的優(yōu)化,結(jié)合MPI的通信機制,充分利用分布式計算中心通信網(wǎng)絡(luò)的特性,可以實現(xiàn)高效的數(shù)據(jù)讀取和分發(fā)。這不是本研究的重點,請參考文獻[2]。

    2.2.3 計算任務(wù)模塊

    采用局部集合最優(yōu)插值的資料同化算法,HYCOM 數(shù)值模式中,地球按照經(jīng)緯度劃分為4,500×3,200 的網(wǎng)格,對周圍一定范圍內(nèi)的觀測數(shù)據(jù),按照網(wǎng)格點逐一進行同化計算,即公式(2)和公式(3)-(7)。因此,資料同化具有良好的計算并行性,將每個網(wǎng)格點的計算作為一個基本的計算任務(wù),劃分為M=4,500×3,200=14,400,000 個基本計算單元,連同數(shù)據(jù),將M 個計算單元按照一定的規(guī)則(如文獻[2]的負載均衡)分配到Np個計算進程上,由Np個計算進程并行完成。文獻[2]中對于全球網(wǎng)格的非海洋網(wǎng)格點進行剔除后,重新進行計算負載的分配,計算效率獲得顯著提升。文獻[2]中通過計算負載均衡整個同化效率得到提升,而有效計算核時并沒有改變,同化計算本身的效率并沒有提高,本文第3部分,提出了一種適合在國產(chǎn)異構(gòu)計算平臺顯著提升同化計算效率(核心是SVD 算法的實現(xiàn))的方法。

    2.3 任務(wù)調(diào)度管理策略

    整個資料同化軟件的操作可以分為:數(shù)據(jù)讀入、數(shù)據(jù)分發(fā)、同化計算和數(shù)據(jù)回寫四個部分。四個部分順序執(zhí)行,而采用流水線的設(shè)計思路,在數(shù)據(jù)讀入運行開始后,當數(shù)據(jù)讀入到一定數(shù)量,啟動數(shù)據(jù)分發(fā),而計算進程在接收到必要的數(shù)據(jù)后,啟動同化計算,當計算進程完成計算后,一次性將結(jié)果數(shù)據(jù)回寫到代理進程,由代理進程將最終的數(shù)據(jù)寫入文件。這樣實現(xiàn)四個部分的操作在時間上的重疊,如圖5所示。

    圖5 任務(wù)模塊流水線操作示意圖Fig.5 Schematic diagram of task module pipeline operation

    3 異構(gòu)計算平臺高效SVD算法的工程實現(xiàn)

    3.1 SVD優(yōu)化策略

    本文的第2部分和文獻[2]中,通過對資料同化的任務(wù)編排、流水線操作和負載均衡進行優(yōu)化,可以有效提高包括資料同化在內(nèi)的并行計算應(yīng)用的計算效率,而資料同化的基本計算單元包括SVD算法在內(nèi)的矩陣運算,計算量大,現(xiàn)有算法并行化程度低,對其進行優(yōu)化可以從基礎(chǔ)算法方面提高計算速度。

    SVD 算法通常采用基于QR 迭代和Jacobi旋轉(zhuǎn)的方法實現(xiàn),文獻[4]、[5]分別對基于Jacobi旋轉(zhuǎn)和QR迭代的方法,單獨使用計算卡或通過CPU-GPU 協(xié)同并行計算實現(xiàn)SVD 算法做了研究。本文通過研究QR迭代的SVD求解算法,提出CPU-GPU 協(xié)同并行計算實現(xiàn)SVD 算法的詳細步驟和驗證結(jié)果。

    對于批量矩陣的SVD算法,優(yōu)化策略如下:

    (1)根據(jù)計算卡雙精度浮點處理器內(nèi)核的數(shù)量、計算卡內(nèi)存數(shù)量和矩陣的大小,確定SVD算法每個批次的矩陣數(shù)量Nm,該數(shù)量通過對計算卡的測試獲得。

    (2)GPU承擔矩陣的二對角化計算任務(wù),并負責Givens變換系數(shù)計算奇異向量U、V矩陣的任務(wù)。

    (3)CPU 承擔QR 迭代,保存每一步Givens變換的系數(shù),用于計算卡計算U、V矩陣。

    (4)合理規(guī)劃和使用計算卡內(nèi)存和主機內(nèi)存,設(shè)計CPU和計算卡的計算任務(wù)調(diào)度機制,實現(xiàn)高效的調(diào)度算法。

    3.2 基于國產(chǎn)異構(gòu)計算平臺的QR迭代SVD算法及實現(xiàn)方法

    基于QR 迭代的方法是求解二對角矩陣奇異值的最快速的算法。矩陣的SVD 求解,可以分為三步:第一步是通過Householder變換,將矩陣轉(zhuǎn)換為二對角矩陣;第二步采用QR 迭代,得到矩陣的奇異值,保存Givens變換的系數(shù);第三步使用Householder 變換時得到U、V,然后再使用QR迭代過程中Givens變換的系數(shù),迭代計算奇異向量U、V矩陣。

    3.2.1 SVD算法原理

    SVD的定義如下:

    一個m×n矩陣,定義A=UΣVT,其中U是一個m×m矩陣,U是一個n×n矩陣,Σ 是一個m×n矩陣,其除主對角線上元素外,其他元素為0,并且主對角線上元素σi>0。

    在資料同化應(yīng)用中,由于m=n,即待分解的矩陣為方陣,且為對稱方陣,因此,U=V。在資料同化的SVD操作,見公式(1)-公式(7)。

    3.2.2 SVD解算的實現(xiàn)流程和數(shù)據(jù)結(jié)構(gòu)

    從QR 迭代法求解矩陣的奇異值和奇異向量的一般算法出發(fā),根據(jù)國產(chǎn)異構(gòu)計算平臺的特點,本文設(shè)計了SVD 解算的實現(xiàn)流程和數(shù)據(jù)結(jié)構(gòu),具體流程如圖6。

    圖6 異構(gòu)計算平臺SVD算法流程Fig.6 SVD algorithm process for heterogeneous computing platforms

    說明:

    (1)CPU 程序負責申請主機內(nèi)存和計算卡的存儲,CPU內(nèi)存使用鎖頁內(nèi)存,從而提高了內(nèi)存訪問的速度。

    (2)從流程圖中可以看到,在SVD 解算過程中,CPU 和計算卡均需要等待,設(shè)CPU 執(zhí)行時間為tcpu,計算卡執(zhí)行時間為tgpu,4 次數(shù)據(jù)交換時間為tdata,1 次批量SVD 的矩陣數(shù)量為n,所以執(zhí)行1 個矩陣SVD 總的運行時間為t=

    (3)在SVD流程中,CPU和GPU均存在等待時間,且等待時間基本上等于對方的執(zhí)行時間。

    為實現(xiàn)高效的SVD 解算,在主機鎖頁內(nèi)存中開辟2份緩沖區(qū),分別存放n個矩陣的數(shù)據(jù)及計算結(jié)果數(shù)據(jù),按照如圖7 的方式執(zhí)行SVD 算法流程。

    通過合理配置批量SVD 矩陣的數(shù)量,使得CPU 執(zhí)行時間、GPU 執(zhí)行時間和數(shù)據(jù)交換時間達到平衡,可以顯著提升批量SVD的計算效率。

    3.2.3 矩陣二對角化(Householder變換)

    Householder 變換最初由A.C Aitken 在1932年提出,是一種線性變換,其變換矩陣被稱作Householder矩陣,其定義具體如下:

    設(shè)ω∈Rn,‖ω‖2=1,定義H=I-2ωωT(H∈Rn×n)。

    Householder變換可以將矩陣的某些元素置零而矩陣的范數(shù)不變,正是利用此性質(zhì),將待分解的矩陣先變?yōu)槎蔷仃嚕粗鲗蔷€和上對角線之外的元素都為零。如下所示3×3 二對角矩陣:

    二對角變換算法如下:

    (1)設(shè)矩陣An×n是n 行n 列的輸入矩陣。如下:

    (2)x=[a11,a21,…,an1],計算v=x±‖x‖2e1,其中e1為除第一個元素為1,其他元素是0的單位向量。

    (3)計算ω=v/‖v‖2。

    (4)計算H=I-2ωωT。

    (5)計算A′n×n=HAn×n,A′n×n為:

    (6)x=[a′12,a′13,…,a′1n],計算v=x±‖x‖2e1。

    (7)計算ω=v/‖v‖2。

    (8)計算H=I-2ωωT。

    (9)計算A″n×n=HA′n×n,A″n×n為:

    (10)將A″n×n從a″22開始,組成新的矩陣An-1×n-1,重復(2)-(9),直到矩陣的大小變?yōu)?×2。

    由于二對角采用計算卡計算,并且二對角過程中的ω組成U 和V 向量矩陣需要保存,按照如下原則進行計算,提升計算速度:

    (1)利用矩陣運算的性質(zhì)減少乘法次數(shù)。具體如下:

    并且在計算ω=v/‖v‖2時從而有:

    可以明顯減少乘法次數(shù),提升計算速度。

    (2)批量計算。以‖v‖2的計算為例,具體如下:

    如果按照常規(guī)計算,設(shè)計一個循環(huán),對v的每個元素進行平方后再求和,每次循環(huán)做一次乘法和一次累加,這樣的計算并沒有充分發(fā)揮GPU 的內(nèi)核緩存性能;優(yōu)化的算法是每個循環(huán)進行k次的乘法和k次的累加,k一般取4,8,16。對于某國產(chǎn)異構(gòu)計算平臺的計算卡,當k=8 時,計算速度最快。在二對角化過程中,對類似的運算(主要是矩陣乘法等運算)進行如是處理,可以明顯提升計算速度。

    3.2.4 QR迭代

    二對角化計算完后,得到二對角矩陣Bn×n和正交向量Un×n、Vn×n,由于二對角化是在異構(gòu)計算平臺的計算卡上運行的,在進行QR 迭代前,要將數(shù)據(jù)搬移到CPU的內(nèi)存中,見圖8,而正交向量Un×n、Vn×n將在QR 迭代后,使用保存的Givens變換的系數(shù)進行計算,得到奇異向量。

    圖8 Givens變換矩陣元素的數(shù)據(jù)結(jié)構(gòu)Fig.8 The data structure of givens transformation matrix elements

    二對角矩陣Bn×n如下所示(以n=3為例):

    對Bn×n沿主對角線和上對角線元素順次進行Givens 變換(具體見相關(guān)文獻),使得次對角線的元素為零,并保存Givens 變換的系數(shù)c和s,不斷迭代,直到將次對角線元素全部變?yōu)榱悖▽嶋H計算時設(shè)置一個閾值,當小于該閾值時就停止迭代)。

    由于QR 迭代存在大量的分支,并且QR 迭代是串行進行的,因此,不適合使用計算卡并行實現(xiàn),該部分計算由更適合進行邏輯判斷的CPU 進行。計算完成后,對得到矩陣的主對角線元素進行如下處理:

    (1)如果主對角線元素大于零,不變;如果小于零,取絕對值,同時,將對應(yīng)的U和V矩陣對應(yīng)向量的元素符號取反。

    (2)對主對角線按照由大到小排序,同時調(diào)整U和V矩陣對應(yīng)向量的順序,以保持一致。

    以上兩步操作在完成奇異向量計算后進行。

    3.2.5 奇異向量的計算

    QR迭代完成后,將保存的Givens變換的系數(shù)搬移到計算卡內(nèi)存,對二對角化后得到的U和V,與Givens變換矩陣相乘,最終得到奇異向量U和V。由于QR 迭代產(chǎn)生的Givens 變換矩陣必須按照相同的順序作用到二對角化后得到的U和V,設(shè)計如圖8數(shù)據(jù)結(jié)構(gòu):

    其中,Nj:第j次迭代的矩陣數(shù)量;SR和ER:進行矩陣運算時,Givens變換矩陣非零部分是2×2 的大小,與U和V對應(yīng)的行列的子矩陣進行運算,SR和ER分別是起始行和結(jié)束行的位置。

    在計算卡中計算奇異矩陣時,通過遍歷該數(shù)據(jù)結(jié)構(gòu),完成奇異向量的高效計算。

    3.3 編程實現(xiàn)語言

    包括資料同化在內(nèi)的大多數(shù)的科學計算算法均由Fortran 語言實現(xiàn),F(xiàn)ortran 的語言特性是接近數(shù)學公式描述,具有很高的執(zhí)行效率,可以直接對矩陣和復數(shù)進行運算。盡管Fortran具備上述優(yōu)點,但是相對于C/C++語言來說,開發(fā)環(huán)境、運行環(huán)境和生態(tài)方面,差距巨大,尤其是對于非科學計算領(lǐng)域的工程人員,不掌握Fortran編程開發(fā)能力。加之某些異構(gòu)計算平臺不支持Fortran開發(fā)環(huán)境,而均支持C/C++開發(fā)和運行環(huán)境。

    鑒于上述原因,本文所涉及的資料同化軟件,全部使用C/C++開發(fā),對于SVD解算核心算法,使用C語言開發(fā),并采用模塊化設(shè)計思路,矩陣乘法、矩陣二對角化、QR分解、向量計算等均是單獨的模塊,封裝成動態(tài)或靜態(tài)鏈接庫,可以通過設(shè)置不同的參數(shù)靈活調(diào)用。經(jīng)實際驗證,采用C語言設(shè)計實現(xiàn)的SVD算法在相同的CPU上的運行速度,與LAPACK庫相比,運行速度相當,而對于采用CPU和計算卡協(xié)同計算相比,具有顯著的性能提升,具體見3.4實測結(jié)果。

    整個程序采用C/C++設(shè)計,可以運行在Linux操作系統(tǒng)和MPI并行環(huán)境下,也可以運行在Windows 操作系統(tǒng)下,Windows 版本提供方便的圖形化操作界面,便于配置和交互。

    3.4 實測結(jié)果

    實際測試的環(huán)境為某國產(chǎn)異構(gòu)計算平臺,每個計算節(jié)點1 顆CPU,2 張國產(chǎn)加速卡。分為SVD核心算法測試和整個同化軟件的測試。

    類GPU環(huán)境與資源配置如下:

    (1)CPU 節(jié)點配置1 塊32 核心的CPU,主頻2.0 GHz 處理器雙精度浮點理論值0.256 TFLOPS。

    (2)異構(gòu)節(jié)點,在CPU節(jié)點基礎(chǔ)上,配置4塊16G 內(nèi)存的類GPU 加速卡,單卡雙精度浮點理論值6.96TFLOPS。

    3.4.1 SVD核心算法測試

    測試方法分為兩種,一種是采用隨機數(shù)測試,具體如下:

    (1)對比測試的庫為rocm的dgesvd_batched函數(shù)。

    (2)采用8,192個100階方陣,矩陣元素為隨機數(shù)。

    測試結(jié)果如表2。

    表2 SVD核心算法對比測試數(shù)據(jù)(方法1)Table 2 SVD core algorithm comparison test data(method 1)

    從實際測試數(shù)據(jù)看,自主開發(fā)的批量SVD算法的計算速度是ROCM 庫的9.51 倍,計算速度顯著提升。

    另一種測試方法是在同化程序優(yōu)化前后,對SVD 解算部分的單獨進行計時,進而測試SVD在國產(chǎn)異構(gòu)計算平臺同化程序中的性能提升,由于采用并行計算,為確保測算準確性,抽取優(yōu)化前后運行時解算矩陣數(shù)量最多的進程的運行時間。具體如表3。

    表3 SVD核心算法對比測試數(shù)據(jù)(方法2)Table 3 SVD core algorithm comparison test data(method 2)

    同化程序中矩陣的大小均為100×100,雖然測試中矩陣元素的數(shù)據(jù)不同,由于矩陣數(shù)量較多,因此,測試數(shù)據(jù)的有效性可以保證。從測試數(shù)據(jù)可知,計算效率提升到23.33 倍,提升效果顯著。

    方法1 和方法2 效率的不同主要原因是兩種測試的數(shù)據(jù)不同,方法1 的數(shù)據(jù)為隨機數(shù)據(jù),迭代次數(shù)多,而方法2 的數(shù)據(jù)在進行二對角化前,主對角線的元素絕對值遠遠大于其他元素,導致對角化后QR的迭代次數(shù)較少。

    3.4.2 同化軟件測試

    整個同化軟件在采用上述的各項優(yōu)化措施后,經(jīng)測試計算性能得到顯著提升。

    運行環(huán)境如下:

    硬件測試環(huán)境:CPU為20顆21-core Processor,本地硬盤為240GB M.2 Nvme SSD,內(nèi)存為8*16GB DDR4 2666MHz,存儲為Parastor300S。

    軟件測試環(huán)境:操作系統(tǒng):CentOS7.6,國產(chǎn)異構(gòu)加速卡環(huán)境:DTK-21.10.1庫,MPI:intelmpi/2020.1.217,編譯器:mpiifort,hipcc。

    測試算例如下:

    基于國產(chǎn)異構(gòu)計算架構(gòu)的全球高分辨率海洋資料同化高效軟件,共分為三個算例,分別是對三個不同時間的海表溫度的海洋資料同化。算例1、2、3分別是對2020年9月29日、2013年9月22日、2012年9月9日的海表溫度進行資料同化。每個算例包括優(yōu)化前程序與優(yōu)化后程序,優(yōu)化前程序是原始的同化程序,優(yōu)化后程序是使用自研batchedsvd 算法以及多種策略優(yōu)化的同化程序。

    具體測試數(shù)據(jù)采用3 個算例平均,具體如表4。

    表4 海洋資料同化優(yōu)化前后性能比較Table 4 Ocean data assimilation performance comparison

    從測試數(shù)據(jù)看,程序運行速度提高到原來的2.5 倍以上,IO 運行速度提升到20 倍左右,計算部分的加速比顯著提升。

    3.5 本文的貢獻和下一步的工作

    本文所涉及的同化軟件及SVD算法方面在大規(guī)模計算應(yīng)用方面的貢獻如下:

    (1)使用C/C++完整實現(xiàn)同化軟件,同化軟件便于可以集成或被集成到其他軟件系統(tǒng),對于豐富大規(guī)模計算尤其是國產(chǎn)異構(gòu)平臺的生態(tài)體系、促進同化軟件的發(fā)展意義重大。

    (2)使用C 語言基于國產(chǎn)異構(gòu)計算平臺,完整實現(xiàn)了SVD 算法,包括矩陣乘法、轉(zhuǎn)置、Householder 變換、Givens 變換、矩陣二對角化、QR迭代,在批量矩陣SVD速度方面取得顯著提升,對于科學計算、人工智能、信號降噪、視頻推薦等領(lǐng)域的工程應(yīng)用意義重大,具體見文獻[9-11]。

    (3)以資料同化應(yīng)用為例,對于數(shù)據(jù)和計算雙密集型的應(yīng)用在國產(chǎn)異構(gòu)平臺的設(shè)計思路、優(yōu)化方法做了一般性的探索,并形成了較為成熟的應(yīng)用開發(fā)經(jīng)驗和指導原則,對于同類型的大規(guī)模計算應(yīng)用開發(fā)具有一定的借鑒意義。

    下一步的工作:

    (1)同化軟件中,除SVD算法外,插值、距離計算等計算任務(wù)的計算量比較大,并且適合使用計算卡,配合篩選算法,實現(xiàn)CPU和計算卡的協(xié)同計算方式,進一步提升計算速度。

    (2)開展以現(xiàn)有SVD解算為基礎(chǔ),面向國產(chǎn)異構(gòu)計算平臺,使用C/C++語言,開發(fā)和豐富跨平臺的線性代數(shù)的相關(guān)算法庫。

    4 總 結(jié)

    本文主要針對國產(chǎn)異構(gòu)計算平臺環(huán)境下的數(shù)據(jù)和計算雙密集型的同化算法的高效實現(xiàn)給出相關(guān)策略、方法和創(chuàng)新,可對類似的應(yīng)用提供借鑒,特別是采用C/C++實現(xiàn)的以SVD 解算為核心的自主開發(fā)線性代數(shù)基本計算庫,豐富了國產(chǎn)異構(gòu)計算平臺的通用計算卡底層生態(tài),擴展了國產(chǎn)異構(gòu)計算平臺的應(yīng)用領(lǐng)域。

    利益沖突聲明

    所有作者聲明不存在利益沖突關(guān)系。

    猜你喜歡
    數(shù)據(jù)結(jié)構(gòu)進程資料
    Party Time
    PAIRS & TWOS
    JUST A THOUGHT
    債券市場對外開放的進程與展望
    中國外匯(2019年20期)2019-11-25 09:54:58
    “翻轉(zhuǎn)課堂”教學模式的探討——以《數(shù)據(jù)結(jié)構(gòu)》課程教學為例
    高職高專數(shù)據(jù)結(jié)構(gòu)教學改革探討
    中國市場(2016年45期)2016-05-17 05:15:48
    TRIZ理論在“數(shù)據(jù)結(jié)構(gòu)”多媒體教學中的應(yīng)用
    社會進程中的新聞學探尋
    民主與科學(2014年3期)2014-02-28 11:23:03
    《數(shù)據(jù)結(jié)構(gòu)》教學方法創(chuàng)新探討
    河南科技(2014年5期)2014-02-27 14:08:57
    我國高等教育改革進程與反思
    五月伊人婷婷丁香| 成人黄色视频免费在线看| 亚洲最大成人手机在线| 亚洲经典国产精华液单| 日韩欧美精品免费久久| 国产探花极品一区二区| 性色avwww在线观看| 国产精品久久久久久久电影| 欧美日韩国产mv在线观看视频 | 性色av一级| 亚洲性久久影院| 3wmmmm亚洲av在线观看| 国产成人aa在线观看| 少妇高潮的动态图| 亚洲高清免费不卡视频| 18+在线观看网站| 黄色视频在线播放观看不卡| 99久久精品国产国产毛片| 嫩草影院新地址| 午夜老司机福利剧场| 国产精品人妻久久久影院| 大香蕉久久网| 草草在线视频免费看| 欧美+日韩+精品| 成人二区视频| 99九九线精品视频在线观看视频| 一个人观看的视频www高清免费观看| 成人黄色视频免费在线看| 成人综合一区亚洲| 视频中文字幕在线观看| 国产老妇女一区| 国产一区二区亚洲精品在线观看| 老司机影院成人| 国产精品不卡视频一区二区| 国产精品一区www在线观看| 精品酒店卫生间| 亚洲一区二区三区欧美精品 | 国产一区二区亚洲精品在线观看| 亚洲三级黄色毛片| 69av精品久久久久久| 午夜日本视频在线| 午夜福利高清视频| 亚洲在久久综合| 日韩欧美一区视频在线观看 | 国产亚洲av片在线观看秒播厂| videossex国产| 青青草视频在线视频观看| 久久久成人免费电影| 伦理电影大哥的女人| 麻豆乱淫一区二区| 插阴视频在线观看视频| 亚洲美女视频黄频| 熟女电影av网| 精品少妇黑人巨大在线播放| 国产乱来视频区| 视频区图区小说| 欧美日韩一区二区视频在线观看视频在线 | 亚洲精品日韩在线中文字幕| 男人添女人高潮全过程视频| 国产亚洲一区二区精品| 色视频www国产| 亚洲成人一二三区av| 国产精品麻豆人妻色哟哟久久| 久久久久国产精品人妻一区二区| 在线天堂最新版资源| 国模一区二区三区四区视频| 国产精品三级大全| av国产精品久久久久影院| 中国美白少妇内射xxxbb| 乱码一卡2卡4卡精品| 久久久久久久亚洲中文字幕| 国产日韩欧美亚洲二区| 久久精品国产亚洲av天美| 欧美极品一区二区三区四区| 午夜精品国产一区二区电影 | 国产日韩欧美在线精品| 国产免费又黄又爽又色| 国产在线男女| 欧美潮喷喷水| 成人免费观看视频高清| 99久久中文字幕三级久久日本| 嫩草影院入口| 最近最新中文字幕大全电影3| 少妇的逼水好多| 欧美人与善性xxx| 国产白丝娇喘喷水9色精品| 国产精品久久久久久久电影| 街头女战士在线观看网站| 亚洲欧美日韩无卡精品| 2018国产大陆天天弄谢| 青青草视频在线视频观看| 久久97久久精品| 欧美成人a在线观看| 亚洲国产最新在线播放| 少妇人妻精品综合一区二区| 在线亚洲精品国产二区图片欧美 | 久久久欧美国产精品| 熟女av电影| 国产一区二区三区综合在线观看 | 国产在线一区二区三区精| 国产亚洲av片在线观看秒播厂| 大又大粗又爽又黄少妇毛片口| 男女无遮挡免费网站观看| 国产在线男女| 国产老妇女一区| 两个人的视频大全免费| 婷婷色综合大香蕉| 黄色怎么调成土黄色| 又爽又黄无遮挡网站| 国产毛片在线视频| 深夜a级毛片| 美女xxoo啪啪120秒动态图| 天堂俺去俺来也www色官网| 国产伦理片在线播放av一区| 国产一区二区在线观看日韩| 男插女下体视频免费在线播放| 久久久久久久久久久丰满| 国产男女超爽视频在线观看| 2018国产大陆天天弄谢| 韩国av在线不卡| 日韩电影二区| .国产精品久久| 交换朋友夫妻互换小说| freevideosex欧美| 亚洲一级一片aⅴ在线观看| 色哟哟·www| 我要看日韩黄色一级片| 啦啦啦中文免费视频观看日本| 青春草亚洲视频在线观看| 亚洲精品亚洲一区二区| 亚洲欧美成人综合另类久久久| 欧美少妇被猛烈插入视频| 日韩欧美精品免费久久| 国产69精品久久久久777片| a级一级毛片免费在线观看| 高清午夜精品一区二区三区| 精品国产乱码久久久久久小说| 国产亚洲最大av| 一级毛片我不卡| 中文字幕免费在线视频6| 免费大片黄手机在线观看| 亚洲第一区二区三区不卡| 国产成人91sexporn| 久久97久久精品| 国产中年淑女户外野战色| 亚洲国产精品国产精品| 80岁老熟妇乱子伦牲交| 国产成人一区二区在线| 插逼视频在线观看| 街头女战士在线观看网站| 高清av免费在线| videossex国产| 亚洲美女搞黄在线观看| 在线精品无人区一区二区三 | 边亲边吃奶的免费视频| 亚洲综合色惰| 一级毛片电影观看| 日韩免费高清中文字幕av| 欧美成人午夜免费资源| 久久鲁丝午夜福利片| 欧美日韩视频精品一区| 久久久精品欧美日韩精品| 精品一区二区三区视频在线| 国产成人aa在线观看| 国产精品久久久久久久久免| 热re99久久精品国产66热6| 国产精品一二三区在线看| 秋霞在线观看毛片| 亚洲av在线观看美女高潮| 亚洲欧洲日产国产| 国精品久久久久久国模美| 嫩草影院入口| 亚洲人与动物交配视频| 国产日韩欧美在线精品| 人人妻人人澡人人爽人人夜夜| 制服丝袜香蕉在线| a级一级毛片免费在线观看| 女人久久www免费人成看片| 不卡视频在线观看欧美| 国产精品伦人一区二区| 日韩成人伦理影院| 天堂中文最新版在线下载 | 久久热精品热| 制服丝袜香蕉在线| 交换朋友夫妻互换小说| 欧美区成人在线视频| 成人亚洲欧美一区二区av| av网站免费在线观看视频| 久久久久久九九精品二区国产| 26uuu在线亚洲综合色| 国产综合精华液| 国产美女午夜福利| 欧美极品一区二区三区四区| 成年女人在线观看亚洲视频 | 激情五月婷婷亚洲| 一本久久精品| 亚洲精品中文字幕在线视频 | 亚洲国产精品专区欧美| 国产成人免费无遮挡视频| 在线观看免费高清a一片| 久久久久久伊人网av| 久久精品综合一区二区三区| 中文乱码字字幕精品一区二区三区| 亚洲av成人精品一区久久| 高清午夜精品一区二区三区| 我的女老师完整版在线观看| 插阴视频在线观看视频| 精品人妻一区二区三区麻豆| 在线天堂最新版资源| 成人国产av品久久久| 成人高潮视频无遮挡免费网站| 精华霜和精华液先用哪个| 国产黄片视频在线免费观看| 九九爱精品视频在线观看| 插逼视频在线观看| 午夜激情久久久久久久| 日本一二三区视频观看| av黄色大香蕉| 国产毛片a区久久久久| av在线观看视频网站免费| 国产精品福利在线免费观看| 亚洲色图av天堂| 神马国产精品三级电影在线观看| 婷婷色综合大香蕉| 人妻夜夜爽99麻豆av| 777米奇影视久久| 毛片女人毛片| 蜜桃久久精品国产亚洲av| 人妻一区二区av| 国产午夜福利久久久久久| 久久久久久九九精品二区国产| xxx大片免费视频| 男人爽女人下面视频在线观看| 国产成人精品一,二区| 亚洲经典国产精华液单| 99热网站在线观看| 亚洲av免费高清在线观看| 又黄又爽又刺激的免费视频.| 日韩国内少妇激情av| .国产精品久久| 自拍偷自拍亚洲精品老妇| 有码 亚洲区| 禁无遮挡网站| 涩涩av久久男人的天堂| 高清午夜精品一区二区三区| 亚洲av成人精品一二三区| 久久影院123| 久久久久精品性色| 免费少妇av软件| 在线看a的网站| av又黄又爽大尺度在线免费看| 国产精品熟女久久久久浪| 欧美潮喷喷水| 街头女战士在线观看网站| 国产成人a∨麻豆精品| 老师上课跳d突然被开到最大视频| 亚洲欧美精品自产自拍| 日本猛色少妇xxxxx猛交久久| 欧美日韩国产mv在线观看视频 | 女人被狂操c到高潮| 国产伦精品一区二区三区四那| 成人亚洲欧美一区二区av| 五月伊人婷婷丁香| 青青草视频在线视频观看| 99热全是精品| 亚洲精品一二三| 我的老师免费观看完整版| 国产伦精品一区二区三区视频9| 国产成人福利小说| 男人爽女人下面视频在线观看| 国产欧美日韩精品一区二区| 亚洲欧美一区二区三区黑人 | 亚洲,一卡二卡三卡| 亚洲在久久综合| 亚洲熟女精品中文字幕| 18禁在线播放成人免费| 国产精品嫩草影院av在线观看| 婷婷色av中文字幕| 搞女人的毛片| 三级国产精品欧美在线观看| 国产伦精品一区二区三区视频9| 91久久精品国产一区二区三区| av女优亚洲男人天堂| 精品视频人人做人人爽| 麻豆久久精品国产亚洲av| 国产成人精品婷婷| 国内精品宾馆在线| 亚洲欧美日韩卡通动漫| 亚洲av免费高清在线观看| 又爽又黄a免费视频| av天堂中文字幕网| 一级毛片黄色毛片免费观看视频| 汤姆久久久久久久影院中文字幕| 久久国产乱子免费精品| 七月丁香在线播放| 最近手机中文字幕大全| 18禁在线无遮挡免费观看视频| 久久人人爽人人片av| 欧美+日韩+精品| 黄片wwwwww| 欧美日韩视频高清一区二区三区二| 韩国av在线不卡| 亚洲精品一二三| 久久国产乱子免费精品| 在线观看国产h片| 婷婷色av中文字幕| 熟女av电影| 97人妻精品一区二区三区麻豆| 欧美日韩视频高清一区二区三区二| 狂野欧美白嫩少妇大欣赏| 亚洲在线观看片| 全区人妻精品视频| 色吧在线观看| 亚洲国产精品成人久久小说| 午夜免费男女啪啪视频观看| 午夜福利网站1000一区二区三区| 综合色av麻豆| 国产精品久久久久久精品电影小说 | 午夜激情久久久久久久| 国产中年淑女户外野战色| kizo精华| 亚洲精品日韩av片在线观看| 欧美成人a在线观看| 成人欧美大片| 3wmmmm亚洲av在线观看| 亚洲国产最新在线播放| 亚洲无线观看免费| 99热网站在线观看| 在线免费观看不下载黄p国产| 国产老妇女一区| 久久ye,这里只有精品| 国产欧美日韩一区二区三区在线 | 九九久久精品国产亚洲av麻豆| 亚洲va在线va天堂va国产| 99热这里只有是精品50| 欧美日韩视频精品一区| 亚洲va在线va天堂va国产| 国产精品国产av在线观看| 91aial.com中文字幕在线观看| 久久韩国三级中文字幕| 国产精品久久久久久精品古装| 99精国产麻豆久久婷婷| 国产国拍精品亚洲av在线观看| 最近手机中文字幕大全| 欧美3d第一页| 国产精品不卡视频一区二区| 午夜福利在线在线| 麻豆国产97在线/欧美| 日本一二三区视频观看| av专区在线播放| 久久久久久九九精品二区国产| 久久久久久久午夜电影| 51国产日韩欧美| 午夜视频国产福利| 中文天堂在线官网| 亚洲一级一片aⅴ在线观看| 老司机影院成人| 亚洲av中文字字幕乱码综合| 欧美日韩在线观看h| 国产黄频视频在线观看| 啦啦啦在线观看免费高清www| 男人爽女人下面视频在线观看| h日本视频在线播放| 七月丁香在线播放| 国产乱来视频区| 精品视频人人做人人爽| 久久久午夜欧美精品| 久久韩国三级中文字幕| 亚洲丝袜综合中文字幕| 亚洲精品成人久久久久久| 毛片一级片免费看久久久久| 国产在线一区二区三区精| 亚洲精品,欧美精品| 国产精品熟女久久久久浪| 亚洲国产精品999| 亚洲伊人久久精品综合| 国内精品宾馆在线| 另类亚洲欧美激情| 国产精品成人在线| 性色avwww在线观看| 赤兔流量卡办理| 日韩欧美精品v在线| 亚洲精品久久午夜乱码| 精品久久久噜噜| 亚洲自拍偷在线| 国语对白做爰xxxⅹ性视频网站| 国产 一区 欧美 日韩| 如何舔出高潮| 夜夜爽夜夜爽视频| 国内少妇人妻偷人精品xxx网站| 国模一区二区三区四区视频| 啦啦啦中文免费视频观看日本| 人妻 亚洲 视频| 免费电影在线观看免费观看| 成人鲁丝片一二三区免费| 国产一区二区三区综合在线观看 | 亚洲精品乱码久久久v下载方式| 纵有疾风起免费观看全集完整版| 综合色av麻豆| 久久久久久久国产电影| 色吧在线观看| 精品久久久精品久久久| 日本一本二区三区精品| 国产爱豆传媒在线观看| 国产成人a∨麻豆精品| av专区在线播放| 久久精品久久久久久噜噜老黄| 国产精品99久久99久久久不卡 | 午夜视频国产福利| 哪个播放器可以免费观看大片| 婷婷色综合www| 在线精品无人区一区二区三 | 亚州av有码| 建设人人有责人人尽责人人享有的 | 大片免费播放器 马上看| 可以在线观看毛片的网站| 国产一区二区在线观看日韩| 亚洲av免费高清在线观看| 九草在线视频观看| 午夜激情久久久久久久| 亚洲av不卡在线观看| 亚洲色图av天堂| 国国产精品蜜臀av免费| 中文字幕亚洲精品专区| 久久热精品热| 日韩强制内射视频| 尤物成人国产欧美一区二区三区| 男人舔奶头视频| 丰满少妇做爰视频| 熟女av电影| 一本色道久久久久久精品综合| 久久久久网色| 黄色配什么色好看| 精品久久久久久久久av| 国产爽快片一区二区三区| 日韩欧美 国产精品| 好男人视频免费观看在线| 另类亚洲欧美激情| 国产成人免费观看mmmm| 狂野欧美白嫩少妇大欣赏| 成人亚洲精品一区在线观看 | 亚洲精品aⅴ在线观看| 亚洲精品久久午夜乱码| 国产欧美日韩精品一区二区| 国产真实伦视频高清在线观看| av在线播放精品| tube8黄色片| 看非洲黑人一级黄片| 亚洲欧美中文字幕日韩二区| 久久亚洲国产成人精品v| 久久久色成人| 女人久久www免费人成看片| 美女国产视频在线观看| 99久国产av精品国产电影| 成人综合一区亚洲| 看黄色毛片网站| 欧美日韩视频精品一区| 在线免费观看不下载黄p国产| 肉色欧美久久久久久久蜜桃 | 男女无遮挡免费网站观看| 亚洲伊人久久精品综合| 又爽又黄无遮挡网站| 交换朋友夫妻互换小说| 国产免费一区二区三区四区乱码| 国产探花极品一区二区| 免费不卡的大黄色大毛片视频在线观看| 成人二区视频| 只有这里有精品99| 99热6这里只有精品| 免费观看av网站的网址| 日本爱情动作片www.在线观看| 色5月婷婷丁香| 一级a做视频免费观看| 男人添女人高潮全过程视频| 18禁裸乳无遮挡动漫免费视频 | 亚洲一级一片aⅴ在线观看| 新久久久久国产一级毛片| 亚洲国产av新网站| 国产精品蜜桃在线观看| 街头女战士在线观看网站| 超碰97精品在线观看| 免费看a级黄色片| 午夜老司机福利剧场| 高清午夜精品一区二区三区| 少妇猛男粗大的猛烈进出视频 | 国产在视频线精品| 欧美日韩一区二区视频在线观看视频在线 | 午夜福利视频精品| 久久久久久久久久久免费av| 午夜亚洲福利在线播放| 中国美白少妇内射xxxbb| 一边亲一边摸免费视频| 97在线视频观看| 久久精品久久久久久久性| 麻豆久久精品国产亚洲av| 大片免费播放器 马上看| 黄色配什么色好看| 一级二级三级毛片免费看| 日韩欧美一区视频在线观看 | 女的被弄到高潮叫床怎么办| 国产乱人视频| 天美传媒精品一区二区| 亚洲三级黄色毛片| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 欧美一级a爱片免费观看看| 少妇人妻久久综合中文| 成人鲁丝片一二三区免费| 国产免费福利视频在线观看| 性色avwww在线观看| 看十八女毛片水多多多| 在线观看美女被高潮喷水网站| 日韩欧美一区视频在线观看 | 秋霞伦理黄片| 久久久久久伊人网av| 国产精品国产三级国产专区5o| 麻豆久久精品国产亚洲av| 国模一区二区三区四区视频| 黄色日韩在线| av在线播放精品| 亚洲在久久综合| 一级毛片久久久久久久久女| 国产精品人妻久久久久久| 精品久久久久久久久av| 国产欧美日韩精品一区二区| 久久人人爽人人爽人人片va| 国产成人a区在线观看| 菩萨蛮人人尽说江南好唐韦庄| 欧美极品一区二区三区四区| 久久精品国产自在天天线| 日本猛色少妇xxxxx猛交久久| 亚洲av电影在线观看一区二区三区 | 国产欧美日韩精品一区二区| 草草在线视频免费看| 少妇熟女欧美另类| 欧美激情久久久久久爽电影| 亚洲精品乱久久久久久| 国产中年淑女户外野战色| 成人免费观看视频高清| 亚洲精品乱久久久久久| 麻豆乱淫一区二区| 不卡视频在线观看欧美| 国产精品一及| 99久国产av精品国产电影| 一级毛片aaaaaa免费看小| 91aial.com中文字幕在线观看| 18禁裸乳无遮挡免费网站照片| 九九爱精品视频在线观看| 国产老妇女一区| 国产一区二区亚洲精品在线观看| 欧美日韩精品成人综合77777| 日日摸夜夜添夜夜爱| 成人漫画全彩无遮挡| 国产乱来视频区| 色综合色国产| 亚洲国产精品专区欧美| 在现免费观看毛片| 一级二级三级毛片免费看| 在线看a的网站| 在线观看美女被高潮喷水网站| 久久这里有精品视频免费| 日韩不卡一区二区三区视频在线| 国产v大片淫在线免费观看| 欧美激情久久久久久爽电影| 校园人妻丝袜中文字幕| 女的被弄到高潮叫床怎么办| 日韩欧美 国产精品| 18禁在线无遮挡免费观看视频| 亚洲人成网站在线观看播放| 精品久久久久久电影网| 免费黄网站久久成人精品| 精品熟女少妇av免费看| 美女cb高潮喷水在线观看| 日韩在线高清观看一区二区三区| 嫩草影院精品99| 国产91av在线免费观看| 99热国产这里只有精品6| 人妻少妇偷人精品九色| 国产爱豆传媒在线观看| av卡一久久| 成人亚洲欧美一区二区av| 午夜福利视频精品| 国产毛片在线视频| 搡老乐熟女国产| 精品人妻一区二区三区麻豆| 国产伦精品一区二区三区视频9| 欧美3d第一页| 哪个播放器可以免费观看大片| 在线 av 中文字幕| 久久久久久久久久人人人人人人| av卡一久久| 人人妻人人爽人人添夜夜欢视频 | 国产精品偷伦视频观看了| a级毛色黄片| 国产男女超爽视频在线观看| 国模一区二区三区四区视频| av福利片在线观看| 国产av国产精品国产| 亚洲国产日韩一区二区| 狂野欧美激情性xxxx在线观看| 亚洲欧美成人精品一区二区| 成人国产麻豆网| 老女人水多毛片| 国产精品一及| 18禁动态无遮挡网站| 亚洲精品成人久久久久久| 男人和女人高潮做爰伦理| 精品熟女少妇av免费看| 国产淫片久久久久久久久| 黄色日韩在线| 精品久久久精品久久久| 五月开心婷婷网| 久久午夜福利片| 亚洲国产高清在线一区二区三| 少妇猛男粗大的猛烈进出视频 | 日韩伦理黄色片|