色哟哟视频在线观看-色哟哟视频在线-色哟哟欧美15最新在线-色哟哟免费在线观看-国产l精品国产亚洲区在线观看-国产l精品国产亚洲区久久

0
  • 聊天消息
  • 系統(tǒng)消息
  • 評論與回復
登錄后你可以
  • 下載海量資料
  • 學習在線課程
  • 觀看技術視頻
  • 寫文章/發(fā)帖/加入社區(qū)
會員中心
創(chuàng)作中心

完善資料讓更多小伙伴認識你,還能領取20積分哦,立即完善>

3天內不再提示

解析OneFlow Element-Wise算子實現方法

jf_pmFSk4VX ? 來源:GiantPandaCV ? 作者:GiantPandaCV ? 2022-12-12 10:54 ? 次閱讀

0x0. 前言

由于CUDA水平太菜,所以一直沒寫過這方面的筆記。現在日常的工作中已經不能離開寫CUDA代碼,所以準備學習ZZK隨緣做一做CUDA的筆記記錄一下學習到的知識和技巧。這篇文章記錄的是閱讀OneFlow的Element-Wise系列CUDA算子實現方案學習到的技巧,希望可以幫助到一起入門CUDA的小伙伴們。Elemet-Wise算子指的是針對輸入Tensor進行逐元素操作,比如ReLU就是針對輸入Tensor的每個值進行判斷是否大于0,大于0的話輸出就是輸入否則就是0。用CUDA來表達最簡單的寫法就是:

__global__voidrelu_kernel(float*input,float*output){
int32_tidx=blockIdx.x*blockDim.x+threadIdx.x;
output[idx]=input[idx]>>(src,dst);

cudaDeviceSynchronize();
cudaFree(src);
cudaFree(dst);
return0;
}

雖然這種寫法非常簡單明了,但卻存在明顯的性能問題。所以這篇文章將基于OneFlow開源的Element-Wise CUDA算子方案來解釋如何寫一個高性能的Element-Wise CUDA算子。

0x1. 性能

以GELU激活函數為例子,分別測試 dtype = float32,不同shape下的前向耗時以及帶寬利用率(NVIDIA A100-PCIE-40GB)。性能情況如下圖所示:

9f2cb390-7987-11ed-8abf-dac502259ad0.png

在這里插入圖片描述

9f2cb390-7987-11ed-8abf-dac502259ad0.png

在這里插入圖片描述

可以看到對于 GeLU 來說,無論是性能還是帶寬 OneFlow 的實現都是更優(yōu)的,接下來我們就來了解一下為什么 OneFlow 的 Element-Wise 算子性能可以做到更優(yōu)。

0x2. 用法

OneFlow在 elementwise.cuh 文件中分別針對一元,二元,三元運算的 Element-Wise 操作實現了模板函數。在包含這個頭文件之后我們可以使用 cuda::Unary/Binary/Ternary 這幾個模板函數來針對我們自己定義的 Element-Wise 操作進行計算。注意,這里說的一元,二元,三元代表的是這個 Element-Wise 操作有幾個輸入 Tensor。

我們舉個例子,假設我們要做的 Element-Wise 操作是逐點乘法,也即有 2 個輸入Tensor x 和 y,然后 x 和 y的形狀和數據類型都是一致的。那么我們可以定義一個模板類:

template
structMultiplyFunctor{
OF_DEVICE_FUNCToperator()(Tx,Ty)const{
returnx*y;
}
};

這里 OF_DEVICE_FUNC 表示我們定義的這個函數既可以運行在 CPU 又可以運行在 GPU 上,它的定義是:

#ifdefined(__CUDACC__)
#defineOF_DEVICE_FUNCTION__device____host____forceinline__
#else
#defineOF_DEVICE_FUNCTIONinline
#endif

然后我們就可以使用 cuda::Binary 這個模板函數來完成這個二元的 Element-Wise 算子了。示例代碼如下:

constuser_op::Tensor*x=ctx->Tensor4ArgNameAndIndex("x",0);
constuser_op::Tensor*y=ctx->Tensor4ArgNameAndIndex("y",0);
user_op::Tensor*out=ctx->Tensor4ArgNameAndIndex("out",0);
constint64_telem_cnt=x->shape().elem_cnt();
OF_CUDA_CHECK(cuda::Binary(MultiplyFunctor(),elem_cnt,out->mut_dptr(),
x->dptr(),
y->dptr(),
ctx->device_ctx()->cuda_stream()));

這里的 x, y, out 分別代表這個 Element-Wise 操作的輸入輸出 Tensor,然后 element_cnt 表示 Tensor 的元素個數,輸出張量的數據首地址 out->mut_dptr(), 輸入張量的數據首地址 x->dptr() && y->dptr() ,最后一個參數則是當前 Kernel 運行的 cuda Stream對象。

0x3. 原理&&代碼實現解析

我個人認為這里有幾個要點,分別是一個線程處理多個數據,向量化數據訪問提升帶寬,設置合理的Block數量(GridSize)和線程數量(BlockSize)以及在合適的地方進行循環(huán)展開(unrool)以及一些編程上的技巧。

0x3.1 給 Element-Wise 操作設置合理的 GridSize 和 BlockSize

下面這段代碼展示了 OneFlow 針對 Element-Wise 算子是如何設置 GridSize 和 BlockSize 的。對應的源碼地址為:https://github.com/Oneflow-Inc/oneflow/blob/master/oneflow/core/cuda/elementwise.cuh#L30-L52 。

constexprintkBlockSize=256;
constexprintkNumWaves=32;

inlinecudaError_tGetNumBlocks(int64_tn,int*num_blocks){
intdev;
{
cudaError_terr=cudaGetDevice(&dev);
if(err!=cudaSuccess){returnerr;}
}
intsm_count;
{
cudaError_terr=cudaDeviceGetAttribute(&sm_count,cudaDevAttrMultiProcessorCount,dev);
if(err!=cudaSuccess){returnerr;}
}
inttpm;
{
cudaError_terr=cudaDeviceGetAttribute(&tpm,cudaDevAttrMaxThreadsPerMultiProcessor,dev);
if(err!=cudaSuccess){returnerr;}
}
*num_blocks=std::max(1,std::min((n+kBlockSize-1)/kBlockSize,
sm_count*tpm/kBlockSize*kNumWaves));
returncudaSuccess;
}

這個地方 BlockSize 直接被設置為了 256 ,對應 constexpr int kBlockSize = 256; 這行代碼,也就是說每個 Block 有 256 個線程。為什么是 256 ?大家不妨讀一下俊丞大佬這篇經典的 給CUDA Kernel設置合適的 GridSize 和 Block Size 的文章 。文章中通過對 SM 的資源分析確定在主流的GPU上將 BlockSize 設置為 128 或者 256 是比較合適,在這里直接設置為了 256 。

確定了 BlockSize 之后需要確定 Kernel 啟動線程塊的數量,我一直覺得上述文章中對這一段的分析是尤其精彩的,這里再截圖展示一下:

9f4990fa-7987-11ed-8abf-dac502259ad0.png

選自OneFlow CUDA Kernel 中 grid_size 和 block_size 應該怎么設置 一文

根據這里的分析,對于 Element-Wise 操作要設置合適的 GridSize 不僅需要考慮元素的數量還要考慮由于 SM 硬件本身帶來的限制。如下公式所述:

*num_blocks=std::max(1,std::min((n+kBlockSize-1)/kBlockSize,
sm_count*tpm/kBlockSize*kNumWaves));

這里的 (n + kBlockSize - 1) / kBlockSize 就是根據 Element-Wise 操作的元素個數來計算需要啟動多少個線程塊,比如在文章開頭的例子中有 = 個元素,那么就一共需要 個線程塊。然后這里以GTX 3080Ti為例,它的SM個數也就是sm_count=80,每個SM最多調度的線程數tpm=1536,那么sm_count * tpm / kBlockSize * kNumWaves = 80 * 1536 / 256 * 32 = 15360,所以在這個例子中我們最終設置的線程塊個數為 588 個。

通過上述講解和分析我們已經確定了啟動 Element-Wise CUDA Kernel 的 GridSize 和 BlockSize。

0x3.2 向量化數據訪問提升帶寬

對于大多數 Element-Wise 算子來說,一般它們的計算量不會太大,所以它們的瓶頸一般在GPU的帶寬上。在 NVIDIA 的性能優(yōu)化博客 https://developer.nvidia.com/blog/cuda-pro-tip-increase-performance-with-vectorized-memory-access/ 中提到,對于很多 CUDA 核函數我們都可以通過向量化數據訪問的方式來提升帶寬受限的 Kernel 的性能,特別是對于架構比較新的 GPU 向量化數據訪問的效果會更加明顯。

在 OneFlow 的 Element-Wise 系列算子中,為了更好的進行向量化的數據訪問,俊丞設計了如下的 Pack 數據結構(代碼位置:https://github.com/Oneflow-Inc/oneflow/blob/master/oneflow/core/cuda/elementwise.cuh#L54-L70):

template
structGetPackType{
usingtype=typenamestd::aligned_storage::type;
};

template
usingPackType=typenameGetPackType::type;

template
unionPack{
static_assert(sizeof(PackType)==sizeof(T)*pack_size,"");
__device__Pack(){
//donothing
}
PackTypestorage;
Telem[pack_size];
};

對GetPackType理解有誤請看知乎的修改后正確版本用了 std::aligned_storage 先聲明了一個內存對齊的數據類型 type ,注意這個 type 的內存長度為 pack_size * sizeof(T) 。然后這里的 T 是我們需要進行 Pack 的數據類型,而 pack_size 則表示我們需要 Pack 的元素個數。接下來我們看到 Pack 聯合體中聲明了 storage 和 elem 兩個數組,它們公用同一段對齊的內存。然后 Pack 聯合體的入口有一個檢查: static_assert(sizeof(PackType) == sizeof(T) * pack_size, ""); 這是用來判斷我們之前聲明的 type 的內存長度是否符合預期。

接下來我們從 https://github.com/Oneflow-Inc/oneflow/blob/master/oneflow/core/cuda/elementwise.cuh#L155-L194 這里可以看到這個 Pack 聯合體主要是用在 Kernel 啟動之前判斷 Element-Wise 操作的輸入輸出 Tensor 對應的數據指針地址是否滿足內存對齊的條件,如果不滿足則這個 Element-Wise 操作無法執(zhí)行數據 Pack 。對應下圖2個畫紅色框的地方。

9f77468a-7987-11ed-8abf-dac502259ad0.png

接下來,OneFlow 定義了真正要執(zhí)行數據 Pack 的數據結構 Packed 并且定義了計算 PackSize 的工具函數。代碼位置為:https://github.com/Oneflow-Inc/oneflow/blob/master/oneflow/core/cuda/elementwise.cuh#L72-L95 。

template
structalignas(sizeof(T)*pack_size)Packed{
__device__Packed(){
//donothing
}
union{
Telem[pack_size];
};
};

constexprintkMaxPackBytes=128/8;
constexprintkMaxPackSize=8;

constexprintMin(inta,intb){returna
constexprintPackSize(){
returnMin(kMaxPackBytes/sizeof(T),kMaxPackSize);
}

template
constexprintPackSize(){
returnMin(PackSize(),PackSize());
}

這里需要注意的是對于 CUDA 來說,最多支持 128 個 bit 的訪問粒度,也就是說 PackSize 的大小不能超過 128 個bit。然后對于各種數據類型來說,Half 數據類型的 bit 數是最少的即 16,所以一次性可以支持 Pack 8個half類型的數據,4個float32的數據,以此類推。所以這里的定義的 kMaxPackSize 表示 128/16=8 ,然后 kMaxPackBytes 則表示最大可以 Pack 的 byte 數 。

請注意區(qū)分 bit 和 byte 。

接下來 https://github.com/Oneflow-Inc/oneflow/blob/master/oneflow/core/cuda/elementwise.cuh#L97-L144 則是真正的為 Element-Wise 操作完成數據 Pack 并執(zhí)行計算。

首先來看這段充滿技巧的代碼:

9f848cbe-7987-11ed-8abf-dac502259ad0.png

在這里插入圖片描述

首先這里定義了一個 HasApply2 類用來判斷是否可以支持一次性Pack 2個 char/int8/half2 類型的元素,這個地方是一個針對 int8/half2/char 數據類型的特殊處理,某些 Element-Wise 算子 Kernel 確實需要支持這種數據類型的計算。也就是說對于 half2 的話,在一個內存訪問粒度里我們其實是可以 Pack 128 / 8 = 16個的。然后用了C++模板元編程的 std::enable_if 來控制針對 half2 類型的特殊 Pack 處理,也就是上圖代碼中的兩個 ApplyPack 函數。可以看到對于 half2 類型的 Element-Wise 操作我們需要給對應的 Functor 定義一個 Apply2 函數,比如對于 Cast 操作的 Functor 定義如下:

template
structCastFunctor{
__device__Tooperator()(Fromfrom)const{returnstatic_cast(from);}
};

template
structCastFunctor::value>::type>{
__device__Tooperator()(halffrom)const{returnstatic_cast(static_cast(from));}

__device__voidApply2(To*to,consthalf*from)const{
constfloat2f2=__half22float2(*reinterpret_cast(from));
to[0]=static_cast(f2.x);
to[1]=static_cast(f2.y);
}
};

0x3.3 啟動 Kernel

我們接下來看一下 Element-Wise 的 Kernel 實現:https://github.com/Oneflow-Inc/oneflow/blob/master/oneflow/core/cuda/elementwise.cuh#L133-L144 。

9f98a0b4-7987-11ed-8abf-dac502259ad0.png

在這里插入圖片描述

在 Kernel 中我們發(fā)現每一個線程實際上處理了多個 Pack 后的數據,也即:for (int64_t i = global_tid; i < n_pack; i += blockDim.x * gridDim.x) 。初學者看到這個循環(huán)也許會比較疑惑,為什么它的步幅是 blockDim.x * gridDim.x ?? 這個 blockDim.x * gridDim.x 表示的是 CUDA 線程網格中的線程總數。假設線程網格中有 1280 個線程,線程 0 將計算元素 0、1280、2560 等。通過使用步幅等于網格大小的循環(huán),確保了 warp 中的所有尋址都是單位步幅,可以獲得最大的內存合并。想了解更多細節(jié)可以查看:https://zhuanlan.zhihu.com/p/571320529 。

除此之外,使用這種技巧的還有個好處就是如果對于 Kernel 中存在每個線程都包含一個公共的操作,那么線程數的增多,也代表著這部分的開銷變大。這個時候我們減少線程的數量并循環(huán)進行處理的話那么這個公共操作的開銷就會更低。

最后,在循環(huán)之外,我們還需要根據傳入的 n_tail 參數,看一下還有沒有因為沒有被 pack_size 整除的剩余元素,如果有的話就單獨調用 functor 進行處理。

0x3.4 unroll

實際上就是代碼中的 #pragma unroll ,這個宏會對我們的 for 循環(huán)做循環(huán)展開,讓更多的指令可以并行執(zhí)行。但容易想到,只有處理的數據沒有前后依賴關系的時候我們可以做。對于大多數的 ElementWise 算子來說一般是滿足這個條件的。

0x3.5 Kernel Launch的細節(jié)

在 https://github.com/Oneflow-Inc/oneflow/blob/master/oneflow/core/cuda/elementwise.cuh#L166-L181 這個位置 OneFlow 展示了 Element-Wise Kernel 的啟動細節(jié),我們簡單注釋一下:

template
cudaError_tLaunchKernel(FactoryTfactory,int64_tn,R*r,constIN*...in,cudaStream_tstream){
constint64_tn_pack=n/pack_size;//根據元素個數和pack_size,計算pack數目,比如1026/4=256。
constint64_ttail_offset=n_pack*pack_size;//如果存在不被整除的情況,我們計算使用pack的偏移量:256*4;
constint64_tn_tail=n-tail_offset;////元素數目-偏移量=剩下的元素個數->1026-1024=2
intnum_blocks;
{
cudaError_terr=GetNumBlocks(n_pack,&num_blocks);//計算線程塊數目
if(err!=cudaSuccess){returnerr;}
}
ApplyGeneric<<>>(
factory,n_pack,reinterpret_cast*>(r),
(reinterpret_cast*>(in))...,n_tail,r+tail_offset,
(in+tail_offset)...);
returncudaPeekAtLastError();
}

0x4. 總結

以上就是我對 OneFlow Element-Wise 系列 CUDA 算子實現的解析,后續(xù)有空會持續(xù)更新學習到的新知識。

審核編輯:郭婷

聲明:本文內容及配圖由入駐作者撰寫或者入駐合作網站授權轉載。文章觀點僅代表作者本人,不代表電子發(fā)燒友網立場。文章及其配圖僅供工程師學習之用,如有內容侵權或者其他違規(guī)問題,請聯系本站處理。 舉報投訴
  • 代碼
    +關注

    關注

    30

    文章

    4802

    瀏覽量

    68738
  • CUDA
    +關注

    關注

    0

    文章

    121

    瀏覽量

    13642

原文標題:【BBuf 的CUDA筆記】一,解析OneFlow Element-Wise 算子實現

文章出處:【微信號:GiantPandaCV,微信公眾號:GiantPandaCV】歡迎添加關注!文章轉載請注明出處。

收藏 人收藏

    評論

    相關推薦

    AUTOSAR通信協議解析 如何實現AUTOSAR通信

    通信協議棧是一個復雜的系統(tǒng),它涵蓋了多種通信方式和模塊,以實現車內ECU之間的高效、可靠的數據交換。以下是對AUTOSAR通信協議的解析實現AUTOSAR通信的方法: 一、AUTOS
    的頭像 發(fā)表于 12-17 14:54 ?725次閱讀

    PLC數據采集模塊的編程方法解析

    PLC數據采集模塊的編程方法主要依賴于所使用的PLC品牌和型號,以及具體的應用場景和需求。以下是對PLC數據采集模塊編程方法的一般性解析: 一、PLC數據采集模塊概述 PLC數據采集模塊(也稱為
    的頭像 發(fā)表于 11-26 13:53 ?293次閱讀

    用ezdsp5535板子實現麥克風輸入語音信息,耳機輸出語音,為什么編譯的時候通不過?

    最近在用ezdsp5535板子實現麥克風輸入語音信息,耳機輸出語音。嘗試使用TI官方的例程aic3204實現,修復各種bug之后,現在還是不好用。 例程aic3204中有兩個函數
    發(fā)表于 10-28 07:40

    電源常用ic腳位解析方法 7腳電源芯片怎么看型號

    電源常用IC腳位解析方法 電源常用IC(集成電路)的腳位解析方法主要依賴于對IC引腳功能的理解,以及參考相關的技術手冊或數據手冊。以下是一些通用的
    的頭像 發(fā)表于 10-07 17:10 ?2284次閱讀

    基于 DSP5509 進行數字圖像處理中 Sobel 算子邊緣檢測的硬件連接電路圖

    以下是基于 DSP5509 進行數字圖像處理中 Sobel 算子邊緣檢測的硬件設計方案: 一、總體架構 圖像采集:使用合適的圖像傳感器,如 CMOS 傳感器,通過相應的接口(如 SPI、I2C 等
    發(fā)表于 09-25 15:25

    摩爾線程攜手智源研究院完成基于Triton的大模型算子庫適配

    里,即成功完成了近60個算子的功能驗證,精度符合交付標準,并實現對Bert-large模型的全面支持。FlagGems算子庫在摩爾線程MUSA架構上展現出了接近手寫算子的計算性能,且性
    的頭像 發(fā)表于 08-02 11:06 ?904次閱讀

    微創(chuàng)軟件推出AI大模型應用平臺WISE

    微創(chuàng)軟件在“2024微創(chuàng)人工智能戰(zhàn)略發(fā)布會”上,正式推出了企業(yè)級AI大模型應用平臺WISE。該平臺以其獨特的技術架構和卓越性能,為企業(yè)開發(fā)AI應用提供了全新的解決方案。
    的頭像 發(fā)表于 05-31 11:31 ?885次閱讀

    Arm發(fā)布全新終端計算子系統(tǒng),加速AI體驗與產品上市

    全球領先的半導體知識產權(IP)提供商Arm控股有限公司(納斯達克股票代碼:ARM)今日正式推出全新的Arm終端計算子系統(tǒng)(CSS),以推動人工智能(AI)體驗的前沿發(fā)展,并助力芯片合作伙伴在構建基于Arm架構的解決方案時實現更高效、更快速的流程,從而加速產品上市。
    的頭像 發(fā)表于 05-30 14:23 ?593次閱讀

    Arm宣布推出終端計算子系統(tǒng)(CSS),提供領先的人工智能體驗

    Arm 控股有限公司(納斯達克股票代碼:ARM,以下簡稱“Arm”)今日宣布推出 Arm? 終端計算子系統(tǒng) (CSS),以提供領先的人工智能 (AI) 體驗,助力芯片合作伙伴更輕松、快速地構建基于 Arm 架構的解決方案,并加速其產品上市進程。
    的頭像 發(fā)表于 05-30 14:11 ?1274次閱讀
    Arm宣布推出終端計<b class='flag-5'>算子</b>系統(tǒng)(CSS),提供領先的人工智能體驗

    微創(chuàng)軟件正式發(fā)布AI大模型應用平臺WISE

    上海2024年5月28日?/美通社/ -- 5月20日,微創(chuàng)軟件召開“2024微創(chuàng)人工智能戰(zhàn)略發(fā)布會”,并正式推出企業(yè)級AI大模型應用平臺WISE(Wicresoft Intelligence
    的頭像 發(fā)表于 05-28 17:18 ?616次閱讀
    微創(chuàng)軟件正式發(fā)布AI大模型應用平臺<b class='flag-5'>WISE</b>

    Arm新Arm Neoverse計算子系統(tǒng)(CSS):Arm Neoverse CSS V3和Arm Neoverse CSS N3

    Arm宣布了兩款新的Arm Neoverse計算子系統(tǒng)(CSS),它們基于“迄今為止最好的一代Neoverse技術”。是什么讓這些新產品在擁擠的計算技術領域脫穎而出? Arm的兩個新Arm
    的頭像 發(fā)表于 04-24 17:53 ?1102次閱讀
    Arm新Arm Neoverse計<b class='flag-5'>算子</b>系統(tǒng)(CSS):Arm Neoverse CSS V3和Arm Neoverse CSS N3

    stm32h750既要實現主機,也要實現從機功能,要怎么實現呢?

    h750的板子實現OTG功能,作為主機識別接入的u盤;作為從機接入電腦時,模擬成一個u盤。 在cubemx中USB_OTG_FS的mode中選擇OTG/Dual_Role_Device,在
    發(fā)表于 03-19 06:46

    基于TPU-MLIR:詳解EinSum的完整處理過程!

    、Reduce。EinSum支持任意多的輸入,只要計算中只包含點乘(element-wise)、廣播(broadcast)、歸約求和(reductionsum)都可以使
    的頭像 發(fā)表于 02-19 13:08 ?722次閱讀
    基于TPU-MLIR:詳解EinSum的完整處理過程!

    三相異步電動機調速的方法有哪些?四種常用方法解析

    三相異步電動機調速的方法有哪些?四種常用方法解析? 三相異步電動機調速的方法有很多種,其中較為常用的包括電壓調制、變頻調速、轉差調速和自耦調速等。下面將對這四種常用
    的頭像 發(fā)表于 02-01 16:24 ?8124次閱讀

    詳細解析二相電機反轉的改變方法

    詳細解析二相電機反轉的改變方法? 二相電機反轉是指通過改變電機的工作方式和接線方式來改變電機的旋轉方向。以下是對二相電機反轉的改變方法的詳細解析。 首先,要了解二相電機的工作原理。二相
    的頭像 發(fā)表于 01-23 14:45 ?2783次閱讀
    主站蜘蛛池模板: 卫生间被教官做好爽HH视频| 99久女女精品视频在线观看| 亚洲AV无码乱码A片无码蜜桃| 国产在线午夜| 精品午夜视频| 女的把腿张开男的往里面插| 日韩亚洲欧美中文高清在线| 亚洲精品乱码久久久久久直播 | 快播电影网址| 日本护士性生活| 亚洲精品久久区二区三区蜜桃臀| 最近免费中文MV在线字幕 | 亚洲男同tv| ankha成人| 国产人妻777人伦精品HD| 老男人粗大猛| 我年轻漂亮的继坶2中字在线播放| 伊人久99久女女视频精品免| 超碰99热在线精品视频| 精品国产九九| 日韩在线av免费视久久| 亚洲精品一区国产欧美| gogogo视频在线观看| 狠狠躁日日躁人人爽| 日本美国群交P片内射捆绑| 一天不停的插BB十几次| 成在线人免费| 久久夜色噜噜噜亚洲AV0000| 收集最新中文国产中文字幕| 18video性欧美19sex高清| 中国女人精69xxxxxx视频| 中文字幕乱偷无码AV蜜桃| 国产GV天堂亚洲国产GV刚刚碰| 久久久午夜精品福利内容| 失禁 调教 刺激 哭喊男男| 2020久久精品永久免费| 国产日韩精品SUV| 欧美午夜理伦三级在线观看| 亚洲无吗视频| 国产精品AV无码免费播放| 女教师杨雪的性荡生活|