Func input = BoundaryConditions::repeat_edge(in);
Func refer = BoundaryConditions::repeat_edge(ref);
// the position in motion vector output tableVar bid_x, bid_y;// variable for motion vectorVar mv_x, mv_y;RDom rblk(0, BLOCK_SIZE, 0, BLOCK_SIZE)
// position in input frame
Expr blk_x = (bid_x * BLOCK_SIZE) + rblk.x;Expr blk_y = (bid_y * BLOCK_SIZE) + rblk.y;
// position in reference frameExpr ref_x = blk_x + mv_x;
Expr ref_y = blk_y + mv_y;
Expr in_val = i16(input(blk_x, blk_y);Expr ref_val = i16(refer(ref_x, ref_y));Func sad;sad(mv_x, mv_y, bid_x, bid_y) = sum(abs(in_val - ref_val));
// setup search window
RDom search_win(-range/2, range, -range/2, range);// argmin will find the one in search_window which get minimal sad.Tuple me = argmin(sad(search_win.x, search_win.y, bid_x, bid_y));// output the motion estimation result to MV tableVar ch;out_mv(bid_x, bid_y, ch) = me[ch];
// AVX-512 versionvoid __foo_avx512(void){ ... }// AVX versionvoid __foo_avx2(void){ ... }// SSE4 versionvoid __foo_sse4(void){ ... }// C versionvoid __foo(void){ ... }void foo(void){if(_avx512_available())__foo_avx512();else if(_avx2_available())__foo_avx2();else if(_sse4_available())__foo_sse4();else__foo();}
// AVX-512 version__attribute__ ((target ("avx512f")))void foo(void){ ... }// AVX2 version__attribute__ ((target ("avx2")))void foo(void){ ... }// SSE4 version__attribute__ ((target ("sse4.2")))void foo(void){ ... }// C version__attribute__ ((target ("default")))void foo(void){ ... }
// C version__attribute__ ((target_clones ("avx512f", "avx2", "sse4.2", "default")))void foo(void){ ... }
#define TSIZE 8接著以 gcc 編譯:
#if defined (__clang__)
typedef float vfloat __attribute__((ext_vector_type(TSIZE)));
#else
typedef float vfloat __attribute__ ((vector_size (TSIZE*4)));
#endif
__attribute__ ((target_clones ("avx512f", "avx2", "avx", "sse4.2", "default")))
void gemm_vec(float *a, int sa, float *b, int sb, float *c, int sc)
{
vfloat vb[TSIZE];
for(int y = 0; y < TSIZE; y++){
vb[y] = *((vfloat*)(b + sb*y));
}
for(int y = 0; y < TSIZE; y++){
vfloat vc = *((vfloat*)(c + sc*y));
vfloat va = *((vfloat*)(a + sa*y));
for(int x = 0; x < TSIZE; x++){
vc += va[x] * vb[x];
}
*((vfloat*)(c + sc*y)) = vc;
}
}
$ gcc -O3 -shared mm.c -o libgemm.so
$ objdump -t libgemm.so
0000000000001720 l F .text 0000000000000c96 gemm_vec.default.4
00000000000023c0 l F .text 00000000000005df gemm_vec.avx512f.0
00000000000029a0 l F .text 00000000000005df gemm_vec.avx2.1
0000000000002f80 l F .text 000000000000060c gemm_vec.avx.2
0000000000003590 l F .text 0000000000000c95 gemm_vec.sse4_2.3
u32x8 c = {
a[0] > b[0] ? a[0] : b[0],
a[1] > b[1] ? a[1] : b[1],
a[2] > b[2] ? a[2] : b[2],
a[3] > b[3] ? a[3] : b[3],
a[4] > b[4] ? a[4] : b[4],
a[5] > b[5] ? a[5] : b[5],
a[6] > b[6] ? a[6] : b[6],
a[7] > b[7] ? a[7] : b[7],
};
u32x8 a, b, c;
...
c = a > b ? a : b;
for(int i = 0; i < 8; i++)
c[i] = a[i] > b[i] ? a[i] : b[i];
u32x8 d = (u32x8)(a > b);
u32x8 c = (d & a) | (~d & b);
u32x8 c = _mm256_max_epu32(a, b);
u32x4 a, b, c;
...
c.s02 = a.s02; // or c.xz = a.xz;
c.s13 = b.s02; // or c.yw = b.xz;
u32x4 a, b, c, vout;
...
vout = fa*a + fb*b + c;
for(int i = 0; i < 4; i++){
assert(out[ofs+i] == vout[i]);
}
C = A x B
import spatial.dsl._
@spatial object MatMult_outer extends SpatialApp {
type X = FixPt[TRUE,_16,_16]
def main(args: Array[String]): Unit = {
// Get sizes for matrices from command line
val m = args(0).to[Int]
val n = args(1).to[Int]
val p = args(2).to[Int]
// Generate data for input matrices A and B, and initialize C
val a = (0::m, 0::p){(i,j) => ((i + j * p) % 8).to[X] }
val b = (0::p, 0::n){(i,j) => ((i + j * n) % 8).to[X] }
val c_init = (0::m, 0::n){(_,_) => 0.to[X] }
// Communicate dimensions to FPGA with ArgIns
val M = ArgIn[Int]
val N = ArgIn[Int]
val P = ArgIn[Int]
setArg(M,m)
setArg(N,n)
setArg(P,p)
// Create pointers to matrices
val A = DRAM[X](M, P)
val B = DRAM[X](P, N)
val C = DRAM[X](M, N)
// Set up parallelizations
val op = 1
val mp = 1
val ip = 1
val bm = 16
val bn = 64
val bp = 64
setMem(A, a)
setMem(B, b)
setMem(C, c_init)
Accel {
// Tile by output regions in C
Foreach(M by bm par op, N by bn par op) { (i,j) =>
val tileC = SRAM[X](bm, bn)
// Prefetch C tile
tileC load C(i::i+bm, j::j+bn par ip)
// Accumulate on top of C tile over all tiles in P dimension
MemFold(tileC)(P by bp) { k =>
val tileA = SRAM[X](bm, bp)
val tileB = SRAM[X](bp, bn)
val accum = SRAM[X](bm, bn)
// Load A and B tiles
Parallel {
tileA load A(i::i+bm, k::k+bp)
tileB load B(k::k+bp, j::j+bn)
}
// Perform matrix multiply on tile
MemReduce(accum)(bp by 1 par mp){ kk =>
val tileC_partial = SRAM[X](bm,bn)
Foreach(bm by 1, bn by 1 par ip){ (ii,jj) =>
tileC_partial(ii,jj) = tileA(ii,kk) * tileB(kk,jj)
}
tileC_partial
}{_+_}
}{_+_}
// Store C tile to DRAM
C(i::i+bm, j::j+bn par ip) store tileC
}
}
// Fetch result on host
val result = getMatrix(C)
// Compute correct answer
val gold = (0::m, 0::n){(i,j) =>
Array.tabulate(p){k => a(i,k) * b(k,j)}.reduce{_+_}
}
// Show results
println(r"expected cksum: ${gold.map(a => a).reduce{_+_}}")
println(r"result cksum: ${result.map(a => a).reduce{_+_}}")
printMatrix(gold, "Gold: ")
printMatrix(result, "Result: ")
assert(gold == result)
}
}
a, b, 與 c_init 等在 Host 端的資料結構, 使用了用以產生一個 Matrix的語法. 精確地說是 (0::m, 0::p){(i,j) => … }. 這能用以產生高達 5D 的資料結構. 為了產生 c_init, 以底線來標示 iterators, 這是由於他們並沒有明確地被使用在此函數, 這在 Scala 中這是相當普遍的作法.C 有被初始化為0. 因此正確的作法是在既有存在於 C 中的數值之上做疊加. 在此應用的最外層 loop 為一個 3-stage coarse-grain 的 pipeline. 而 prefetch 行為第一個 stage.k = 0) 覆蓋掉 prefetch 至 tileC 中的數值, 因此這裡選擇使用 MemFold controller. 並且建立了一個 local accum sram, 用來存放兩個來自A 與 B 的 tile 乘法的部份結果. 整個 MemFold 組成了最外部 pipeline 的第二個 stage.kk 自 tileA 選擇了一個欄位以及自 tileB 的一列. 對於 MemReduce 與 MemFold controllers 必須針對其 map 函數傳回一個 SRAM, 所以為此目的建構了 tileC_partial .tileA 的所有列以及 tileB 的所有欄. 這有效進行在由 kk 選出的兩向量中的外積運算(窮舉的成對乘法).MemReduce 回傳此 SRAM, 所以此特定 MemReduce 為外部 MemFold map 函數的輸出. 這也是為何有著兩個顯得奇特多餘的 {_+_} lambdas.tileC 傳輸回 DRAM 的存放動作. 下圖動畫示範如何以 triple-buffered 記憶體的方式思考 tileC.
getMatrix 與 printMatrix 函式來自 DRAM 取得資料並顯示這些資料結構的內容. 同樣地, 能夠使用 getTensor3, getTensor4, 與 getTensor5, 作為搭配 print 印出高維度的資料結構.import spatial.dsl._
@spatial object MatMult_outer extends SpatialApp {
type X = FixPt[TRUE,_16,_16]
def main(args: Array[String]): Unit = {
val m = args(0).to[Int]
val n = args(1).to[Int]
val p = args(2).to[Int]
val a = (0::m, 0::p){(i,j) => ((i + j * p) % 8).to[X] }
val b = (0::p, 0::n){(i,j) => ((i + j * n) % 8).to[X] }
val c_init = (0::m, 0::n){(_,_) => 0.to[X] }
val M = ArgIn[Int]
val N = ArgIn[Int]
val P = ArgIn[Int]
setArg(M,m)
setArg(N,n)
setArg(P,p)
val A = DRAM[X](M, P)
val B = DRAM[X](P, N)
val C = DRAM[X](M, N)
// *** Set mp and ip > 1
val op = 1
val mp = 2
val ip = 4
val bm = 16
val bn = 64
val bp = 64
setMem(A, a)
setMem(B, b)
setMem(C, c_init)
Accel {
Foreach(M by bm par op, N by bn par op) { (i,j) =>
val tileC = SRAM[X](bm, bn)
tileC load C(i::i+bm, j::j+bn par ip)
MemFold(tileC)(P by bp) { k =>
val tileA = SRAM[X](bm, bp)
val tileB = SRAM[X](bp, bn)
val accum = SRAM[X](bm, bn)
Parallel {
tileA load A(i::i+bm, k::k+bp)
// *** Parallelize writer by 8
tileB load B(k::k+bp, j::j+bn par 8)
}
MemReduce(accum)(bp by 1 par mp){ kk =>
val tileC_partial = SRAM[X](bm,bn)
Foreach(bm by 1, bn by 1 par ip){ (ii,jj) =>
tileC_partial(ii,jj) = tileA(ii,kk) * tileB(kk,jj)
}
tileC_partial
}{_+_}
}{_+_}
C(i::i+bm, j::j+bn par ip) store tileC
}
}
val result = getMatrix(C)
val gold = (0::m, 0::n){(i,j) =>
Array.tabulate(p){k => a(i,k) * b(k,j)}.reduce{_+_}
}
println(r"expected cksum: ${gold.map(a => a).reduce{_+_}}")
println(r"result cksum: ${result.map(a => a).reduce{_+_}}")
printMatrix(gold, "Gold: ")
printMatrix(result, "Result: ")
assert(gold == result)
}
}
tileB 作為示範記憶體. Banking 的計算方式基於 “Theory and algorithm for generalized memory partitioning in high-level synthesis” (Wang et al), 並有著次要的修改以及在 pattern 搜尋上的改進.mp = 2 與 ip = 4. 若檢查上面使用到此平行的應用, 將發現 mp 平行化外部的 pipe 而 ip 平行化內部的 pipe. 下圖動畫顯示一個對於 tileB 在此平行化(僅考慮唯讀)下可能的 banking 方案. 請注意這並非唯一可行的方案, 其他分層與平面的方案是存在的.
tileB 的 writer. 這表示上面的 banking 方案不再可行. 需要靜態對 SRAM 作 bank 切分來保證在記憶體存取中沒有任何 bank 衝突. 在這樣的情況中, compiler 可能會選擇如同下圖動畫顯示的平坦的 banking 方案. 請注意這動畫僅顯示以 mp = 2 與 ip = 1 所設置的 reader, 而其他不同 reader 能夠嘗試調整這些參數, 藉此來看對於記憶體會發生什麼改變.
Function Multi-Versioning (FMV) 過往的 CPU 發展歷程中, x86 平台由於因應各種應用需求的提出, 而陸陸續續加入了不同的指令集, 此外也可能因為針對市場做等級區隔, 支援的數量與種類也不等. 在 Linux 平台上這些 CPU 資訊可以透過...
