diff --git a/README.md b/README.md index 8ab07fdf6c..fa6da2a9d9 100644 --- a/README.md +++ b/README.md @@ -1799,7 +1799,7 @@ Early results are shown [here](results/README.md) Authored and maintained by Zheming Jin (https://github.com/zjin-lcf) ## Acknowledgement -Abhishek Bagusetty, Andrew Barker, Andrey Alekseenko, Anton Gorshkov, Beau Johnston, Bernhard Esslinger, Bert de Jong, Chengjian Liu, Chris Knight, Daine McNiven, David Oro, David Sanchez, Douglas Franz, Edson Borin, Gabriell Araujo, Georgi Mirazchiyski, Henry Linjamäki, Henry Gabb, Hugh Delaney, Ian Karlin, Istvan Reguly, Jack Kirk, Jason Lau, Jeff Hammond, Jenny Chen, Jeongnim Kim, Jianxin Qiu, Jakub Chlanda, Jiya Su, John Tramm, Ju Zheng, Junchao Zhang, Kali Uday Balleda, Kinman Lei, Luke Drummond, Manish Motwani, Mark Stephenson, Martin Burtscher, Matthias Noack, Matthew Davis, Michael Kruse, Michel Migdal, Mike Franusich, Mike Giles, Mikhail Dvorskiy, Mohammed Alser, Mohammed ,Tarek Ibn Ziad, Muhammad Haseeb, Muaaz Awan, Nevin Liber, Nicholas Miller, Oscar Ludwig, Pavel Samolysov, Pedro Valero Lara, Piotr Różański, Rahulkumar Gayatri, Samyak Gangwal, Shaoyi Peng, Steffen Larsen, Rafal Bielski, Robert Harrison, Robin Kobus, Rod Burns, Rodrigo Vimieiro, Romanov Vlad, Tadej Ciglarič, Thomas Applencourt, Tiago Carneiro, Tiago Cogumbreiro, Timmie Smith, Tobias Baumann, Usman Roshan, Wayne Joubert, Ye Luo, Yongbin Gu, Zhe Chen +Abhishek Bagusetty, Andrew Barker, Andrey Alekseenko, Anton Gorshkov, Beau Johnston, Bernhard Esslinger, Bert de Jong, Chengjian Liu, Chris Knight, Daine McNiven, David Oro, David Sanchez, Douglas Franz, Edson Borin, Gabriell Araujo, Georgi Mirazchiyski, Henry Linjamäki, Henry Gabb, Hugh Delaney, Ian Karlin, Istvan Reguly, Jack Kirk, Jason Lau, Jeff Hammond, Jenny Chen, Jeongnim Kim, Jianxin Qiu, Jakub Chlanda, Jiya Su, John Tramm, Ju Zheng, Junchao Zhang, Kali Uday Balleda, Kinman Lei, Luke Drummond, Manish Motwani, Mark Stephenson, Martin Burtscher, Matthias Noack, Matthew Davis, Michael Kruse, Michel Migdal, Mike Franusich, Mike Giles, Mikhail Dvorskiy, Mohammed Alser, Mohammed ,Tarek Ibn Ziad, Muhammad Haseeb, Muaaz Awan, Nevin Liber, Nicholas Miller, Oscar Ludwig, Pavel Samolysov, Pedro Valero Lara, Piotr Różański, Rahulkumar Gayatri, Samyak Gangwal, Shaoyi Peng, Steffen Larsen, Rafal Bielski, Robert Harrison, Robin Kobus, Rod Burns, Rodrigo Vimieiro, Romanov Vlad, Tadej Ciglarič, Thomas Applencourt, Tiago Carneiro, Tiago Cogumbreiro, Timmie Smith, Tobias Baumann, Usman Roshan, Wayne Joubert, Ye Luo, Yongbin Gu, Zhe Chen, Zhihui Du Codeplay® and Intel® for their contributions to the oneAPI ecosystem diff --git a/benchmarks.yaml b/benchmarks.yaml index 04fdd73199..5a6e6da9d3 100644 --- a/benchmarks.yaml +++ b/benchmarks.yaml @@ -63,7 +63,7 @@ adv: models: [cuda, hip, omp, sycl] test: regex: '(?:elapsed time=)([0-9.+-e]+)(?: )' - args: ["16", "16", "16"] + args: ["7", "15", "16"] timeout: 300 aes: diff --git a/src/background-subtract-cuda/main.cu b/src/background-subtract-cuda/main.cu index 18987cbcb4..b3a32b2971 100644 --- a/src/background-subtract-cuda/main.cu +++ b/src/background-subtract-cuda/main.cu @@ -8,6 +8,7 @@ #define BLOCK_SIZE 256 +__launch_bounds__(BLOCK_SIZE) __global__ void findMovingPixels( const size_t imgSize, const unsigned char *__restrict__ Img, @@ -26,6 +27,7 @@ __global__ void findMovingPixels( } // alpha = 0.92 +__launch_bounds__(BLOCK_SIZE) __global__ void updateBackground( const size_t imgSize, const unsigned char *__restrict__ Img, @@ -38,6 +40,7 @@ __global__ void updateBackground( } // alpha = 0.92, c = 3 +__launch_bounds__(BLOCK_SIZE) __global__ void updateThreshold( const size_t imgSize, const unsigned char *__restrict__ Img, @@ -56,6 +59,7 @@ __global__ void updateThreshold( // // merge three kernels into a single kernel // +__launch_bounds__(BLOCK_SIZE) __global__ void merge( const size_t imgSize, const unsigned char *__restrict__ Img, diff --git a/src/background-subtract-hip/main.cu b/src/background-subtract-hip/main.cu index 586f8c2019..d987dc9913 100644 --- a/src/background-subtract-hip/main.cu +++ b/src/background-subtract-hip/main.cu @@ -8,6 +8,7 @@ #define BLOCK_SIZE 256 +__launch_bounds__(BLOCK_SIZE) __global__ void findMovingPixels( const size_t imgSize, const unsigned char *__restrict__ Img, @@ -26,6 +27,7 @@ __global__ void findMovingPixels( } // alpha = 0.92 +__launch_bounds__(BLOCK_SIZE) __global__ void updateBackground( const size_t imgSize, const unsigned char *__restrict__ Img, @@ -38,6 +40,7 @@ __global__ void updateBackground( } // alpha = 0.92, c = 3 +__launch_bounds__(BLOCK_SIZE) __global__ void updateThreshold( const size_t imgSize, const unsigned char *__restrict__ Img, @@ -56,6 +59,7 @@ __global__ void updateThreshold( // // merge three kernels into a single kernel // +__launch_bounds__(BLOCK_SIZE) __global__ void merge( const size_t imgSize, const unsigned char *__restrict__ Img, diff --git a/src/background-subtract-sycl/main.cpp b/src/background-subtract-sycl/main.cpp index 5a4f6e96e8..57164c8ef5 100644 --- a/src/background-subtract-sycl/main.cpp +++ b/src/background-subtract-sycl/main.cpp @@ -158,7 +158,7 @@ int main(int argc, char* argv[]) { auto start = std::chrono::steady_clock::now(); q.submit([&] (sycl::handler &cgh) { cgh.parallel_for( - sycl::nd_range<1>(gws, lws), [=] (sycl::nd_item<1> item) { + sycl::nd_range<1>(gws, lws), [[sycl::reqd_work_group_size(BLOCK_SIZE)]] [=] (sycl::nd_item<1> item) { merge ( item, imgSize, d_Img, d_Img1, d_Img2, d_Tn, d_Bn ); }); }).wait(); @@ -169,19 +169,19 @@ int main(int argc, char* argv[]) { auto start = std::chrono::steady_clock::now(); q.submit([&] (sycl::handler &cgh) { cgh.parallel_for( - sycl::nd_range<1>(gws, lws), [=] (sycl::nd_item<1> item) { + sycl::nd_range<1>(gws, lws), [[sycl::reqd_work_group_size(BLOCK_SIZE)]] [=] (sycl::nd_item<1> item) { findMovingPixels ( item, imgSize, d_Img, d_Img1, d_Img2, d_Tn, d_Mp ); }); }); q.submit([&] (sycl::handler &cgh) { cgh.parallel_for( - sycl::nd_range<1>(gws, lws), [=] (sycl::nd_item<1> item) { + sycl::nd_range<1>(gws, lws), [[sycl::reqd_work_group_size(BLOCK_SIZE)]] [=] (sycl::nd_item<1> item) { updateBackground ( item, imgSize, d_Img, d_Mp, d_Bn ); }); }); q.submit([&] (sycl::handler &cgh) { cgh.parallel_for( - sycl::nd_range<1>(gws, lws), [=] (sycl::nd_item<1> item) { + sycl::nd_range<1>(gws, lws), [[sycl::reqd_work_group_size(BLOCK_SIZE)]] [=] (sycl::nd_item<1> item) { updateThreshold ( item, imgSize, d_Img, d_Mp, d_Bn, d_Tn ); }); }); diff --git a/src/bitonic-sort-cuda/main.cu b/src/bitonic-sort-cuda/main.cu index 413732b76e..ad5be1d5ca 100644 --- a/src/bitonic-sort-cuda/main.cu +++ b/src/bitonic-sort-cuda/main.cu @@ -43,6 +43,7 @@ #define BLOCK_SIZE 256 +__launch_bounds__(BLOCK_SIZE) __global__ void bitonic_sort (const int seq_len, const int two_power, int *a) { diff --git a/src/bitonic-sort-hip/main.cu b/src/bitonic-sort-hip/main.cu index 6f97dcbe80..b9d4489be6 100644 --- a/src/bitonic-sort-hip/main.cu +++ b/src/bitonic-sort-hip/main.cu @@ -43,6 +43,7 @@ #define BLOCK_SIZE 256 +__launch_bounds__(BLOCK_SIZE) __global__ void bitonic_sort (const int seq_len, const int two_power, int *a) { diff --git a/src/bitonic-sort-sycl/main.cpp b/src/bitonic-sort-sycl/main.cpp index 8c78e7bdc8..4827c045bb 100644 --- a/src/bitonic-sort-sycl/main.cpp +++ b/src/bitonic-sort-sycl/main.cpp @@ -76,7 +76,7 @@ void ParallelBitonicSort(int input[], int n) { // Offload the work to kernel. q.submit([&](sycl::handler &h) { h.parallel_for(sycl::nd_range<1>( - sycl::range<1>(size), sycl::range<1>(BLOCK_SIZE)), [=](sycl::nd_item<1> item) { + sycl::range<1>(size), sycl::range<1>(BLOCK_SIZE)), [[sycl::reqd_work_group_size(BLOCK_SIZE)]] [=](sycl::nd_item<1> item) { int i = item.get_global_id(0); // Assign the bitonic sequence number. diff --git a/src/channelSum-cuda/main.cu b/src/channelSum-cuda/main.cu index db54e0d2c4..d09f5c01d9 100644 --- a/src/channelSum-cuda/main.cu +++ b/src/channelSum-cuda/main.cu @@ -49,6 +49,7 @@ __device__ void BlockReduce(T &input1, T &input2) { template __global__ +__launch_bounds__(kBlockDimX * kBlockDimY) void ChannelSumNCHW( const int N, const int C, @@ -79,6 +80,7 @@ void ChannelSumNCHW( template __global__ +__launch_bounds__(NUM_THREADS) void ChannelSumNHWC( const int N, const int C, diff --git a/src/channelSum-hip/main.cu b/src/channelSum-hip/main.cu index 33a8a8786c..29ad0270f6 100644 --- a/src/channelSum-hip/main.cu +++ b/src/channelSum-hip/main.cu @@ -49,6 +49,7 @@ __device__ void BlockReduce(T &input1, T &input2) { template __global__ +__launch_bounds__(kBlockDimX * kBlockDimY) void ChannelSumNCHW( const int N, const int C, @@ -78,6 +79,7 @@ void ChannelSumNCHW( template __global__ +__launch_bounds__(NUM_THREADS) void ChannelSumNHWC( const int N, const int C, diff --git a/src/channelSum-sycl/main.cpp b/src/channelSum-sycl/main.cpp index 51952daec9..fe7ca29ccd 100644 --- a/src/channelSum-sycl/main.cpp +++ b/src/channelSum-sycl/main.cpp @@ -19,7 +19,7 @@ typedef int scalar_t; sycl::range<2> lws (128, 1); \ q.submit([&](sycl::handler &cgh) { \ cgh.parallel_for( \ - sycl::nd_range<2>(gws, lws), [=] (sycl::nd_item<2> item) { \ + sycl::nd_range<2>(gws, lws), [=] [[sycl::reqd_work_group_size(128, 1)]] (sycl::nd_item<2> item) { \ Func(item, __VA_ARGS__); \ }); \ }); \ @@ -28,7 +28,7 @@ typedef int scalar_t; sycl::range<2> lws (64, 2); \ q.submit([&](sycl::handler &cgh) { \ cgh.parallel_for( \ - sycl::nd_range<2>(gws, lws), [=] (sycl::nd_item<2> item) { \ + sycl::nd_range<2>(gws, lws), [=] [[sycl::reqd_work_group_size(64, 2)]] (sycl::nd_item<2> item) { \ Func(item, __VA_ARGS__); \ }); \ }); \ @@ -37,7 +37,7 @@ typedef int scalar_t; sycl::range<2> lws (32, 4); \ q.submit([&](sycl::handler &cgh) { \ cgh.parallel_for( \ - sycl::nd_range<2>(gws, lws), [=] (sycl::nd_item<2> item) { \ + sycl::nd_range<2>(gws, lws), [=] [[sycl::reqd_work_group_size(32, 4)]] (sycl::nd_item<2> item) { \ Func(item, __VA_ARGS__); \ }); \ }); \ @@ -46,7 +46,7 @@ typedef int scalar_t; sycl::range<2> lws (16, 8); \ q.submit([&](sycl::handler &cgh) { \ cgh.parallel_for( \ - sycl::nd_range<2>(gws, lws), [=] (sycl::nd_item<2> item) { \ + sycl::nd_range<2>(gws, lws), [=] [[sycl::reqd_work_group_size(16, 8)]] (sycl::nd_item<2> item) { \ Func(item, __VA_ARGS__); \ }); \ }); \ diff --git a/src/dpid-cuda/Makefile b/src/dpid-cuda/Makefile index 4152565799..c8c1ab4e52 100644 --- a/src/dpid-cuda/Makefile +++ b/src/dpid-cuda/Makefile @@ -37,7 +37,7 @@ endif # Optimization Flags ifeq ($(OPTIMIZE),yes) - CFLAGS += -O3 + CFLAGS += -O3 -use_fast_math endif #=============================================================================== # Targets to Build @@ -46,7 +46,10 @@ endif $(program): $(obj) Makefile $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) -%.o: %.cu shared.h Makefile +main.o: main.cu shared.h reference.h Makefile + $(CC) $(CFLAGS) -c $< -o $@ + +kernels.o: kernels.cu shared.h Makefile $(CC) $(CFLAGS) -c $< -o $@ clean: diff --git a/src/dpid-cuda/main.cu b/src/dpid-cuda/main.cu index 1556fb4802..121ca099d1 100644 --- a/src/dpid-cuda/main.cu +++ b/src/dpid-cuda/main.cu @@ -8,6 +8,7 @@ #include #include #include "shared.h" +#include "reference.h" double LCG_random_double(uint64_t * seed) { @@ -70,14 +71,35 @@ int main(int argc, char** argv) { // run downsampling on a device run(p, hInput, hOutput); - int x = 0, y = 0, z = 0; + // per-pixel comparison for validation + uchar3 *hRef = (uchar3*) malloc(sizeof(uchar3) * p.oWidth * p.oHeight); + reference(p.oWidth, p.oHeight, p.iWidth, p.iHeight, + p.pWidth, p.pHeight, p.lambda, + hInput, hRef); + uint32_t mismatches = 0; + int maxDiff = 0; for (uint32_t i = 0; i < p.oWidth * p.oHeight; i++) { - x += hOutput[i].x; - y += hOutput[i].y; - z += hOutput[i].z; + int dx = std::abs((int)hOutput[i].x - (int)hRef[i].x); + int dy = std::abs((int)hOutput[i].y - (int)hRef[i].y); + int dz = std::abs((int)hOutput[i].z - (int)hRef[i].z); + if (dx > 0 || dy > 0 || dz > 0) mismatches++; + if (dx > maxDiff) maxDiff = dx; + if (dy > maxDiff) maxDiff = dy; + if (dz > maxDiff) maxDiff = dz; + } + uint32_t total = p.oWidth * p.oHeight; + if (mismatches == 0) { + printf("Verification PASS: GPU matches CPU reference exactly.\n"); + } else { + printf("Verification: %u / %u pixels differ (%.2f%%), max channel diff = %d\n", + mismatches, total, 100.0 * mismatches / total, maxDiff); + if (maxDiff <= 1) + printf(" -> within rounding tolerance (maxDiff<=1, likely OK)\n"); + else + printf(" -> WARNING: differences exceed rounding tolerance\n"); } - printf("Checksums %d %d %d\n", x, y, z); + free(hRef); free(hInput); free(hOutput); return 0; diff --git a/src/dpid-cuda/reference.h b/src/dpid-cuda/reference.h new file mode 100644 index 0000000000..0a619bc286 --- /dev/null +++ b/src/dpid-cuda/reference.h @@ -0,0 +1,179 @@ +// Pass 1 (guidance): weighted average color per output patch +// Pass 2 (downsampling): L2-distance reweighted downsampling with lambda + +#ifndef DPID_REF_TYPES_DEFINED +#define DPID_REF_TYPES_DEFINED + +struct ref_uchar3 { unsigned char x, y, z; }; +struct ref_float4 { float x, y, z, w; }; + +static inline ref_uchar3 make_ref_uchar3(unsigned char x, unsigned char y, + unsigned char z) { + ref_uchar3 v; v.x = x; v.y = y; v.z = z; return v; +} +static inline ref_float4 make_ref_float4(float x, float y, float z, float w) { + ref_float4 v; v.x = x; v.y = y; v.z = z; v.w = w; return v; +} + +#endif // DPID_REF_TYPES_DEFINED + +static inline void ref_normalize(ref_float4& v) { + v.x /= v.w; v.y /= v.w; v.z /= v.w; v.w = 1.f; +} + +static inline void ref_add(ref_float4& out, unsigned char cx, unsigned char cy, + unsigned char cz, float f) { + out.x += cx * f; out.y += cy * f; out.z += cz * f; out.w += f; +} + +static inline float ref_contribution(float sx, float ex, float sy, float ey, + float f, uint32_t x, uint32_t y) { + if (x < sx) f *= 1.f - (sx - x); + if ((x+1.f) > ex) f *= 1.f - ((x+1.f) - ex); + if (y < sy) f *= 1.f - (sy - y); + if ((y+1.f) > ey) f *= 1.f - ((y+1.f) - ey); + return f; +} + +static inline float ref_lambda(float lambda_val, float dist) { + if (lambda_val == 0.f) return 1.f; + if (lambda_val == 1.f) return dist; + return powf(dist, lambda_val); +} + +static inline float ref_l2dist(const ref_float4& avg, + unsigned char cx, unsigned char cy, + unsigned char cz) { + float dx = avg.x - cx, dy = avg.y - cy, dz = avg.z - cz; + return std::sqrt(dx*dx + dy*dy + dz*dz) / 441.6729559f; +} + +static void ref_guidance(const unsigned char* input_bytes, // uchar3 array + unsigned char* patches_bytes, // uchar3 array + uint32_t oWidth, uint32_t oHeight, + uint32_t iWidth, + float pWidth, float pHeight) { + for (uint32_t py = 0; py < oHeight; py++) { + for (uint32_t px = 0; px < oWidth; px++) { + float sx = std::fmax( px * pWidth, 0.f); + float ex = std::fmin((px+1) * pWidth, (float)iWidth); + float sy = std::fmax( py * pHeight, 0.f); + float ey = std::fmin((py+1) * pHeight, (float)(oHeight * pHeight)); + + uint32_t sxr = (uint32_t)std::floor(sx), syr = (uint32_t)std::floor(sy); + uint32_t exr = (uint32_t)std::ceil(ex), eyr = (uint32_t)std::ceil(ey); + uint32_t xCount = exr - sxr, yCount = eyr - syr; + + ref_float4 color = make_ref_float4(0.f, 0.f, 0.f, 0.f); + + for (uint32_t i = 0; i < xCount * yCount; i++) { + uint32_t x = sxr + (i % xCount), y = syr + (i / xCount); + const unsigned char* pix = input_bytes + (x + y * iWidth) * 3; + float f = ref_contribution(sx, ex, sy, ey, 1.f, x, y); + ref_add(color, pix[0], pix[1], pix[2], f); + } + + ref_normalize(color); + unsigned char* out = patches_bytes + (px + py * oWidth) * 3; + out[0] = (unsigned char)color.x; + out[1] = (unsigned char)color.y; + out[2] = (unsigned char)color.z; + } + } +} + +static ref_float4 ref_calc_average(const unsigned char* patches_bytes, + uint32_t oWidth, uint32_t oHeight, + uint32_t px, uint32_t py) { + const float corner = 1.f, edge = 2.f, center = 4.f; + ref_float4 avg = make_ref_float4(0.f, 0.f, 0.f, 0.f); + + auto add_patch = [&](uint32_t qx, uint32_t qy, float w) { + const unsigned char* p = patches_bytes + (qx + qy * oWidth) * 3; + ref_add(avg, p[0], p[1], p[2], w); + }; + + if (py > 0) { + if (px > 0) add_patch(px-1, py-1, corner); + add_patch(px, py-1, edge); + if (px+1 < oWidth) add_patch(px+1, py-1, corner); + } + if (px > 0) add_patch(px-1, py, edge); + add_patch(px, py, center); + if (px+1 < oWidth) add_patch(px+1, py, edge); + if (py+1 < oHeight) { + if (px > 0) add_patch(px-1, py+1, corner); + add_patch(px, py+1, edge); + if (px+1 < oWidth) add_patch(px+1, py+1, corner); + } + + ref_normalize(avg); + return avg; +} + +static void ref_downsampling(const unsigned char* input_bytes, + const unsigned char* patches_bytes, + unsigned char* output_bytes, + uint32_t oWidth, uint32_t oHeight, + uint32_t iWidth, + float pWidth, float pHeight, + float lambda_val) { + for (uint32_t py = 0; py < oHeight; py++) { + for (uint32_t px = 0; px < oWidth; px++) { + float sx = std::fmax( px * pWidth, 0.f); + float ex = std::fmin((px+1) * pWidth, (float)iWidth); + float sy = std::fmax( py * pHeight, 0.f); + float ey = std::fmin((py+1) * pHeight, (float)(oHeight * pHeight)); + + uint32_t sxr = (uint32_t)std::floor(sx), syr = (uint32_t)std::floor(sy); + uint32_t exr = (uint32_t)std::ceil(ex), eyr = (uint32_t)std::ceil(ey); + uint32_t xCount = exr - sxr, yCount = eyr - syr; + + ref_float4 avg = ref_calc_average(patches_bytes, oWidth, oHeight, px, py); + ref_float4 color = make_ref_float4(0.f, 0.f, 0.f, 0.f); + + for (uint32_t i = 0; i < xCount * yCount; i++) { + uint32_t x = sxr + (i % xCount), y = syr + (i / xCount); + const unsigned char* pix = input_bytes + (x + y * iWidth) * 3; + float f = ref_l2dist(avg, pix[0], pix[1], pix[2]); + f = ref_lambda(lambda_val, f); + f = ref_contribution(sx, ex, sy, ey, f, x, y); + ref_add(color, pix[0], pix[1], pix[2], f); + } + + unsigned char* out = output_bytes + (px + py * oWidth) * 3; + if (color.w == 0.f) { + out[0] = (unsigned char)avg.x; + out[1] = (unsigned char)avg.y; + out[2] = (unsigned char)avg.z; + } else { + ref_normalize(color); + out[0] = (unsigned char)color.x; + out[1] = (unsigned char)color.y; + out[2] = (unsigned char)color.z; + } + } + } +} + +// hInput and hOutput are arrays of uchar3 (3 bytes per pixel, packed). +// Call as: reference(p, hInput, hOutput) +// where hInput/hOutput are cast to void* from uchar3* in the caller. + +void reference(uint32_t oWidth, uint32_t oHeight, + uint32_t iWidth, uint32_t iHeight, + float pWidth, float pHeight, + float lambda_val, + const void* hInput, void* hOutput) { + unsigned char* patches = + (unsigned char*) malloc(oWidth * oHeight * 3); + + ref_guidance((const unsigned char*)hInput, patches, + oWidth, oHeight, iWidth, pWidth, pHeight); + + ref_downsampling((const unsigned char*)hInput, patches, + (unsigned char*)hOutput, + oWidth, oHeight, iWidth, pWidth, pHeight, lambda_val); + + free(patches); +} diff --git a/src/dpid-hip/Makefile b/src/dpid-hip/Makefile index 4eca2a79e4..ba255a5918 100644 --- a/src/dpid-hip/Makefile +++ b/src/dpid-hip/Makefile @@ -23,7 +23,7 @@ obj = $(source:.cu=.o) #=============================================================================== # Standard Flags -CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../dpid-cuda # Linker Flags LDFLAGS = @@ -36,7 +36,7 @@ endif # Optimization Flags ifeq ($(OPTIMIZE),yes) - CFLAGS += -O3 + CFLAGS += -O3 -ffast-math endif #=============================================================================== # Targets to Build @@ -45,7 +45,10 @@ endif $(program): $(obj) Makefile $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) -%.o: %.cu shared.h Makefile +main.o: main.cu shared.h ../dpid-cuda/reference.h Makefile + $(CC) $(CFLAGS) -c $< -o $@ + +kernels.o: kernels.cu shared.h Makefile $(CC) $(CFLAGS) -c $< -o $@ clean: diff --git a/src/dpid-hip/main.cu b/src/dpid-hip/main.cu index 1556fb4802..121ca099d1 100644 --- a/src/dpid-hip/main.cu +++ b/src/dpid-hip/main.cu @@ -8,6 +8,7 @@ #include #include #include "shared.h" +#include "reference.h" double LCG_random_double(uint64_t * seed) { @@ -70,14 +71,35 @@ int main(int argc, char** argv) { // run downsampling on a device run(p, hInput, hOutput); - int x = 0, y = 0, z = 0; + // per-pixel comparison for validation + uchar3 *hRef = (uchar3*) malloc(sizeof(uchar3) * p.oWidth * p.oHeight); + reference(p.oWidth, p.oHeight, p.iWidth, p.iHeight, + p.pWidth, p.pHeight, p.lambda, + hInput, hRef); + uint32_t mismatches = 0; + int maxDiff = 0; for (uint32_t i = 0; i < p.oWidth * p.oHeight; i++) { - x += hOutput[i].x; - y += hOutput[i].y; - z += hOutput[i].z; + int dx = std::abs((int)hOutput[i].x - (int)hRef[i].x); + int dy = std::abs((int)hOutput[i].y - (int)hRef[i].y); + int dz = std::abs((int)hOutput[i].z - (int)hRef[i].z); + if (dx > 0 || dy > 0 || dz > 0) mismatches++; + if (dx > maxDiff) maxDiff = dx; + if (dy > maxDiff) maxDiff = dy; + if (dz > maxDiff) maxDiff = dz; + } + uint32_t total = p.oWidth * p.oHeight; + if (mismatches == 0) { + printf("Verification PASS: GPU matches CPU reference exactly.\n"); + } else { + printf("Verification: %u / %u pixels differ (%.2f%%), max channel diff = %d\n", + mismatches, total, 100.0 * mismatches / total, maxDiff); + if (maxDiff <= 1) + printf(" -> within rounding tolerance (maxDiff<=1, likely OK)\n"); + else + printf(" -> WARNING: differences exceed rounding tolerance\n"); } - printf("Checksums %d %d %d\n", x, y, z); + free(hRef); free(hInput); free(hOutput); return 0; diff --git a/src/dpid-sycl/Makefile b/src/dpid-sycl/Makefile index ac27c77963..a996fd9f17 100644 --- a/src/dpid-sycl/Makefile +++ b/src/dpid-sycl/Makefile @@ -30,8 +30,8 @@ obj = $(source:.cpp=.o) #=============================================================================== # Standard Flags -CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../include -fsycl \ - --gcc-toolchain=$(GCC_TOOLCHAIN) +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../dpid-cuda \ + -fsycl --gcc-toolchain=$(GCC_TOOLCHAIN) ifeq ($(VENDOR), AdaptiveCpp) CFLAGS_TMP := $(CFLAGS) @@ -59,7 +59,7 @@ endif # Optimization Flags ifeq ($(OPTIMIZE),yes) - CFLAGS += -O3 + CFLAGS += -O3 -ffast-math endif ifeq ($(GPU),yes) @@ -72,7 +72,10 @@ endif $(program): $(obj) $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) -%.o: %.cpp shared.h +main.o: main.cpp shared.h ../dpid-cuda/reference.h + $(CC) $(CFLAGS) -c $< -o $@ + +kernels.o: kernels.cpp shared.h $(CC) $(CFLAGS) -c $< -o $@ clean: diff --git a/src/dpid-sycl/kernels.cpp b/src/dpid-sycl/kernels.cpp index 3f093ac4e6..75f7463e6e 100644 --- a/src/dpid-sycl/kernels.cpp +++ b/src/dpid-sycl/kernels.cpp @@ -19,17 +19,17 @@ inline void normalize(float4& var) { - var.x() /= var.w(); - var.y() /= var.w(); - var.z() /= var.w(); - var.w() = 1.f; + var[0] /= var[3]; + var[1] /= var[3]; + var[2] /= var[3]; + var[3] = 1.f; } void add(float4& output, const uchar3& color, const float factor) { - output.x() += (float)color.x() * factor; - output.y() += (float)color.y() * factor; - output.z() += (float)color.z() * factor; - output.w() += factor; + output[0] += (float)color[0] * factor; + output[1] += (float)color[1] * factor; + output[2] += (float)color[2] * factor; + output[3] += factor; } inline @@ -76,10 +76,10 @@ float4 __shfl_down(sycl::nd_item<2> &item, const float4 var, const uint32_t srcL { float4 output; auto sg = item.get_sub_group(); - output.x() = sycl::shift_group_left(sg, var.x(), srcLane); - output.y() = sycl::shift_group_left(sg, var.y(), srcLane); - output.z() = sycl::shift_group_left(sg, var.z(), srcLane); - output.w() = sycl::shift_group_left(sg, var.w(), srcLane); + output[0] = sycl::shift_group_left(sg, var[0], srcLane); + output[1] = sycl::shift_group_left(sg, var[1], srcLane); + output[2] = sycl::shift_group_left(sg, var[2], srcLane); + output[3] = sycl::shift_group_left(sg, var[3], srcLane); return output; } @@ -94,9 +94,9 @@ void reduce(sycl::nd_item<2> &item, float4& value) { inline float distance(const float4& avg, const uchar3& color) { - const float x = avg.x() - color.x(); - const float y = avg.y() - color.y(); - const float z = avg.z() - color.z(); + const float x = avg[0] - color[0]; + const float y = avg[1] - color[1]; + const float z = avg[2] - color[2]; return sycl::sqrt(x * x + y * y + z * z) / 441.6729559f; // L2-Norm / sqrt(255^2 * 3) } @@ -118,7 +118,7 @@ void kernelGuidance(sycl::nd_item<2> &item, float f = contribution(l, 1.f, x, y); const uchar3& pixel = input[x + y * p.iWidth]; - color += {pixel.x() * f, pixel.y() * f, pixel.z() * f, f}; + color += {pixel[0] * f, pixel[1] * f, pixel[2] * f, f}; } // reduce warps @@ -127,7 +127,7 @@ void kernelGuidance(sycl::nd_item<2> &item, // store results if((TX % 32) == 0) { normalize(color); - patches[PX + PY * p.oWidth] = uchar3(color.x(), color.y(), color.z()); + patches[PX + PY * p.oWidth] = uchar3(color[0], color[1], color[2]); } } @@ -213,11 +213,11 @@ void kernelDownsampling(sycl::nd_item<2> &item, if(WTHREAD == 0) { uchar3& ref = output[PX + PY * p.oWidth]; - if(color.w() == 0.0f) - ref = uchar3((unsigned char)avg.x(), (unsigned char)avg.y(), (unsigned char)avg.z()); + if(color[3] == 0.0f) + ref = uchar3((unsigned char)avg[0], (unsigned char)avg[1], (unsigned char)avg[2]); else { normalize(color); - ref = uchar3((unsigned char)color.x(), (unsigned char)color.y(), (unsigned char)color.z()); + ref = uchar3((unsigned char)color[0], (unsigned char)color[1], (unsigned char)color[2]); } } } @@ -249,14 +249,14 @@ void run(const Params& p, const uchar3* hInput, uchar3* hOutput) { for (uint32_t i = 0; i < p.repeat; i++) { q.submit([&] (sycl::handler &cgh) { cgh.parallel_for(sycl::nd_range<2>(gws, lws), - [=] (sycl::nd_item<2> item) [[intel::reqd_sub_group_size(32)]] { + [=] (sycl::nd_item<2> item) { kernelGuidance (item, dInput, dGuidance, p); }); }); q.submit([&] (sycl::handler &cgh) { cgh.parallel_for(sycl::nd_range<2>(gws, lws), - [=] (sycl::nd_item<2> item) [[intel::reqd_sub_group_size(32)]] { + [=] (sycl::nd_item<2> item) { kernelDownsampling (item, dInput, dGuidance, p, dOutput); }); }); diff --git a/src/dpid-sycl/main.cpp b/src/dpid-sycl/main.cpp index 772659b102..990d255b32 100644 --- a/src/dpid-sycl/main.cpp +++ b/src/dpid-sycl/main.cpp @@ -8,10 +8,7 @@ #include #include #include "shared.h" - -//------------------------------------------------------------------- -// SHARED -//------------------------------------------------------------------- +#include "reference.h" double LCG_random_double(uint64_t * seed) { @@ -24,7 +21,6 @@ double LCG_random_double(uint64_t * seed) void run(const Params& i, const uchar3* hInput, uchar3* hOutput); -//------------------------------------------------------------------- int main(int argc, char** argv) { // generate a random image for testing if(argc != 5) { @@ -66,23 +62,44 @@ int main(int argc, char** argv) { uchar3 *hOutput = (uchar3*) malloc (sizeof(uchar3) * p.oWidth * p.oHeight); uint64_t seed = 123; - for (int i = 0; i < p.iWidth * p.iHeight; i++) { - hInput[i].x() = (unsigned char)(256*LCG_random_double(&seed)); - hInput[i].y() = (unsigned char)(256*LCG_random_double(&seed)); - hInput[i].z() = (unsigned char)(256*LCG_random_double(&seed)); + for (uint32_t i = 0; i < p.iWidth * p.iHeight; i++) { + hInput[i][0] = (unsigned char)(256*LCG_random_double(&seed)); + hInput[i][1] = (unsigned char)(256*LCG_random_double(&seed)); + hInput[i][2] = (unsigned char)(256*LCG_random_double(&seed)); } // run downsampling on a device run(p, hInput, hOutput); - int x = 0, y = 0, z = 0; - for (int i = 0; i < p.oWidth * p.oHeight; i++) { - x += hOutput[i].x(); - y += hOutput[i].y(); - z += hOutput[i].z(); + // per-pixel comparison for validation + uchar3 *hRef = (uchar3*) malloc(sizeof(uchar3) * p.oWidth * p.oHeight); + reference(p.oWidth, p.oHeight, p.iWidth, p.iHeight, + p.pWidth, p.pHeight, p.lambda, + hInput, hRef); + uint32_t mismatches = 0; + int maxDiff = 0; + for (uint32_t i = 0; i < p.oWidth * p.oHeight; i++) { + int dx = std::abs((int)hOutput[i][0] - (int)hRef[i][0]); + int dy = std::abs((int)hOutput[i][1] - (int)hRef[i][1]); + int dz = std::abs((int)hOutput[i][2] - (int)hRef[i][2]); + if (dx > 0 || dy > 0 || dz > 0) mismatches++; + if (dx > maxDiff) maxDiff = dx; + if (dy > maxDiff) maxDiff = dy; + if (dz > maxDiff) maxDiff = dz; + } + uint32_t total = p.oWidth * p.oHeight; + if (mismatches == 0) { + printf("Verification PASS: GPU matches CPU reference exactly.\n"); + } else { + printf("Verification: %u / %u pixels differ (%.2f%%), max channel diff = %d\n", + mismatches, total, 100.0 * mismatches / total, maxDiff); + if (maxDiff <= 1) + printf(" -> within rounding tolerance (maxDiff<=1, likely OK)\n"); + else + printf(" -> WARNING: differences exceed rounding tolerance\n"); } - printf("Checksums %d %d %d\n", x, y, z); + free(hRef); free(hInput); free(hOutput); return 0; diff --git a/src/f16sp-cuda/kernels.h b/src/f16sp-cuda/kernels.h index 3242fda362..df72925506 100644 --- a/src/f16sp-cuda/kernels.h +++ b/src/f16sp-cuda/kernels.h @@ -1,5 +1,6 @@ #define NUM_OF_BLOCKS (1024 * 1024) #define NUM_OF_THREADS 128 +// __launch_bounds__ tells the compiler the block size for register allocation // fp16 reduction __forceinline__ __device__ @@ -29,6 +30,7 @@ void reduceInShared_native(float2 * const v) } +__launch_bounds__(NUM_OF_THREADS) __global__ void scalarProductKernel_intrinsics( half2 const *__restrict__ const a, @@ -36,12 +38,13 @@ void scalarProductKernel_intrinsics( float *__restrict__ const results, size_t const size) { - const int stride = gridDim.x*blockDim.x; + const size_t stride = (size_t)gridDim.x * blockDim.x; __shared__ half2 shArray[NUM_OF_THREADS]; half2 value = __float2half2_rn(0.f); - for (int i = threadIdx.x + blockDim.x * blockIdx.x; i < size; i+=stride) + for (size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + i < size; i+=stride) { value = __hfma2(a[i], b[i], value); } @@ -59,6 +62,7 @@ void scalarProductKernel_intrinsics( } // compute type is float for high precision +__launch_bounds__(NUM_OF_THREADS) __global__ void scalarProductKernel_native_fp32( half2 const *__restrict__ const a, @@ -66,12 +70,13 @@ void scalarProductKernel_native_fp32( float *__restrict__ const results, size_t const size) { - const int stride = gridDim.x*blockDim.x; + const size_t stride = (size_t)gridDim.x * blockDim.x; __shared__ float2 shArray[NUM_OF_THREADS]; float2 value = {0.f, 0.f}; - for (int i = threadIdx.x + blockDim.x * blockIdx.x; i < size; i+=stride) + for (size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + i < size; i+=stride) { value.x += (float)a[i].x * (float)b[i].x; value.y += (float)a[i].y * (float)b[i].y; @@ -90,6 +95,7 @@ void scalarProductKernel_native_fp32( // CUB reduction with float2 seems unsupported +__launch_bounds__(NUM_OF_THREADS) __global__ void scalarProductKernel_native2_fp32( half2 const *__restrict__ const a, @@ -97,14 +103,15 @@ void scalarProductKernel_native2_fp32( float *__restrict__ const results, size_t const size) { - const int stride = gridDim.x*blockDim.x; + const size_t stride = (size_t)gridDim.x * blockDim.x; typedef cub::BlockReduce BlockReduce; __shared__ typename BlockReduce::TempStorage temp_storage; float2 value = {0.f, 0.f}; - for (int i = threadIdx.x + blockDim.x * blockIdx.x; i < size; i+=stride) + for (size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + i < size; i+=stride) { value.x += (float)a[i].x * (float)b[i].x; value.y += (float)a[i].y * (float)b[i].y; diff --git a/src/f16sp-cuda/main.cu b/src/f16sp-cuda/main.cu index 51a187dd64..a4ab42e3e6 100644 --- a/src/f16sp-cuda/main.cu +++ b/src/f16sp-cuda/main.cu @@ -12,12 +12,28 @@ #include #include #include -#include +#include #include #include +#include +#include #include #include "kernels.h" +const int WARMUP = 1000; + +struct Half2Mul { + const half2* a; + const half2* b; + Half2Mul(const half2* _a, const half2* _b) : a(_a), b(_b) {} + + __host__ __device__ + float operator()(int idx) const { + half2 prod = __hmul2(a[idx], b[idx]); + return __low2float(prod) + __high2float(prod); + } +}; + // The analytical result of dot product is 65504 void generateInput(half2 * a, size_t size) { @@ -71,8 +87,8 @@ int main(int argc, char *argv[]) printf("\nGPU grid size is %d\n", grid); - // warmup - for (int i = 0; i < 1000; i++) + // WARMUP + for (int i = 0; i < WARMUP; i++) scalarProductKernel_intrinsics<<>>(d_a, d_b, d_r, size); cudaDeviceSynchronize(); @@ -91,8 +107,8 @@ int main(int argc, char *argv[]) cudaMemcpy(&r, d_r, result_bytes, cudaMemcpyDeviceToHost); printf("Error rate: %e\n", fabsf(r - 65504.f)/65504.f); - // warmup - for (int i = 0; i < 1000; i++) + // WARMUP + for (int i = 0; i < WARMUP; i++) scalarProductKernel_native_fp32<<>>(d_a, d_b, d_r, size); cudaDeviceSynchronize(); @@ -111,8 +127,8 @@ int main(int argc, char *argv[]) cudaMemcpy(&r, d_r, result_bytes, cudaMemcpyDeviceToHost); printf("Error rate: %e\n", fabsf(r - 65504.f)/65504.f); - // warmup - for (int i = 0; i < 1000; i++) + // WARMUP + for (int i = 0; i < WARMUP; i++) scalarProductKernel_native2_fp32<<>>(d_a, d_b, d_r, size); cudaDeviceSynchronize(); @@ -142,25 +158,51 @@ int main(int argc, char *argv[]) xType = yType = rType = CUDA_R_16F; eType = CUDA_R_32F; - // warmup - for (int i = 0; i < 1000; i++) { - cublasDotEx(h, size*2, (half*)d_a, xType, 1, (half*)d_b, + // WARMUP + for (int i = 0; i < WARMUP; i++) { + cublasDotEx_64(h, size*2, (half*)d_a, xType, 1, (half*)d_b, yType, 1, d_r2, rType, eType); } cudaDeviceSynchronize(); auto start = std::chrono::steady_clock::now(); for (int i = 0; i < repeat; i++) { - cublasDotEx(h, size*2, (half*)d_a, xType, 1, (half*)d_b, + cublasDotEx_64(h, size*2, (half*)d_a, xType, 1, (half*)d_b, yType, 1, d_r2, rType, eType); } cudaDeviceSynchronize(); auto end = std::chrono::steady_clock::now(); auto time = std::chrono::duration_cast(end - start).count(); - printf("Average kernel (cublasDot) execution time %f (us)\n", (time * 1e-3f) / repeat); + printf("Average cublasDotEx_64 execution time %f (us)\n", (time * 1e-3f) / repeat); cudaMemcpy(&r2, d_r2, result2_bytes, cudaMemcpyDeviceToHost); printf("Error rate: %e\n", fabsf((float)r2 - 65504.f)/65504.f); + printf("\n"); + + for (int i = 0; i < WARMUP + repeat; i++) { + if (i == WARMUP) { + cudaDeviceSynchronize(); + start = std::chrono::steady_clock::now(); + } + + thrust::counting_iterator counting(0); + auto iter = thrust::make_transform_iterator(counting, Half2Mul(d_a, d_b)); + + void* d_temp = nullptr; size_t temp_bytes = 0; + cub::DeviceReduce::Sum(d_temp, temp_bytes, iter, d_r, size); + if (temp_bytes != 0) cudaMalloc(&d_temp, temp_bytes); + auto status = cub::DeviceReduce::Sum(d_temp, temp_bytes, iter, d_r, size); + if (status != cudaSuccess) printf("cub::DeviceReduce::Sum failed\n"); + if (temp_bytes != 0) cudaFree(d_temp); + } + cudaDeviceSynchronize(); + end = std::chrono::steady_clock::now(); + time = std::chrono::duration_cast(end - start).count(); + printf("Average cub::DeviceReduce::Sum execution time: %f (us)\n", time * 1e-3f / repeat); + + cudaMemcpy(&r, d_r, result_bytes, cudaMemcpyDeviceToHost); + printf("Error rate: %e\n", fabsf(r - 65504.f)/65504.f); + cudaFree(d_a); cudaFree(d_b); cudaFree(d_r); diff --git a/src/f16sp-hip/kernels.h b/src/f16sp-hip/kernels.h index 9e773abdc3..ce604632a3 100644 --- a/src/f16sp-hip/kernels.h +++ b/src/f16sp-hip/kernels.h @@ -1,5 +1,6 @@ #define NUM_OF_BLOCKS (1024 * 1024) #define NUM_OF_THREADS 128 +// __launch_bounds__ tells the compiler the block size for register allocation // fp16 reduction __forceinline__ __device__ @@ -29,6 +30,7 @@ void reduceInShared_native(float2 * const v) } +__launch_bounds__(NUM_OF_THREADS) __global__ void scalarProductKernel_intrinsics( half2 const *__restrict__ const a, @@ -36,12 +38,13 @@ void scalarProductKernel_intrinsics( float *__restrict__ const results, size_t const size) { - const int stride = gridDim.x*blockDim.x; + const size_t stride = (size_t)gridDim.x * blockDim.x; __shared__ half2 shArray[NUM_OF_THREADS]; half2 value = __float2half2_rn(0.f); - for (int i = threadIdx.x + blockDim.x * blockIdx.x; i < size; i+=stride) + for (size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + i < size; i+=stride) { value = __hfma2(a[i], b[i], value); } @@ -59,6 +62,7 @@ void scalarProductKernel_intrinsics( } // compute type is float for high precision +__launch_bounds__(NUM_OF_THREADS) __global__ void scalarProductKernel_native_fp32( half2 const *__restrict__ const a, @@ -66,12 +70,13 @@ void scalarProductKernel_native_fp32( float *__restrict__ const results, size_t const size) { - const int stride = gridDim.x*blockDim.x; + const size_t stride = (size_t)gridDim.x * blockDim.x; __shared__ float2 shArray[NUM_OF_THREADS]; float2 value = {0.f, 0.f}; - for (int i = threadIdx.x + blockDim.x * blockIdx.x; i < size; i+=stride) + for (size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + i < size; i+=stride) { value.x += (float)a[i].x * (float)b[i].x; value.y += (float)a[i].y * (float)b[i].y; @@ -90,6 +95,7 @@ void scalarProductKernel_native_fp32( // CUB reduction with float2 +__launch_bounds__(NUM_OF_THREADS) __global__ void scalarProductKernel_native2_fp32( half2 const *__restrict__ const a, @@ -97,14 +103,15 @@ void scalarProductKernel_native2_fp32( float *__restrict__ const results, size_t const size) { - const int stride = gridDim.x*blockDim.x; + const size_t stride = (size_t)gridDim.x * blockDim.x; typedef hipcub::BlockReduce BlockReduce; __shared__ typename BlockReduce::TempStorage temp_storage; float2 value = {0.f, 0.f}; - for (int i = threadIdx.x + blockDim.x * blockIdx.x; i < size; i+=stride) + for (size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + i < size; i+=stride) { value.x += (float)a[i].x * (float)b[i].x; value.y += (float)a[i].y * (float)b[i].y; diff --git a/src/f16sp-hip/main.cu b/src/f16sp-hip/main.cu index f7234c123f..e71c4b3b7d 100644 --- a/src/f16sp-hip/main.cu +++ b/src/f16sp-hip/main.cu @@ -18,6 +18,21 @@ #include #include "kernels.h" +const int WARMUP = 1000; + +struct Half2Mul { + const half2* a; + const half2* b; + Half2Mul(const half2* _a, const half2* _b) : a(_a), b(_b) {} + + __host__ __device__ + float operator()(int idx) const { + half2 prod = __hmul2(a[idx], b[idx]); + return __low2float(prod) + __high2float(prod); + } +}; + + // The analytical result of dot product is 65504 void generateInput(half2 * a, size_t size) { @@ -71,8 +86,8 @@ int main(int argc, char *argv[]) printf("\nGPU grid size is %d\n", grid); - // warmup - for (int i = 0; i < 1000; i++) + // WARMUP + for (int i = 0; i < WARMUP; i++) scalarProductKernel_intrinsics<<>>(d_a, d_b, d_r, size); hipDeviceSynchronize(); @@ -91,8 +106,8 @@ int main(int argc, char *argv[]) hipMemcpy(&r, d_r, result_bytes, hipMemcpyDeviceToHost); printf("Error rate: %e\n", fabsf(r - 65504.f)/65504.f); - // warmup - for (int i = 0; i < 1000; i++) + // WARMUP + for (int i = 0; i < WARMUP; i++) scalarProductKernel_native_fp32<<>>(d_a, d_b, d_r, size); hipDeviceSynchronize(); @@ -111,8 +126,8 @@ int main(int argc, char *argv[]) hipMemcpy(&r, d_r, result_bytes, hipMemcpyDeviceToHost); printf("Error rate: %e\n", fabsf(r - 65504.f)/65504.f); - // warmup - for (int i = 0; i < 1000; i++) + // WARMUP + for (int i = 0; i < WARMUP; i++) scalarProductKernel_native2_fp32<<>>(d_a, d_b, d_r, size); hipDeviceSynchronize(); @@ -142,8 +157,8 @@ int main(int argc, char *argv[]) xType = yType = rType = HIP_R_16F; eType = HIP_R_32F; - // warmup - for (int i = 0; i < 1000; i++) { + // WARMUP + for (int i = 0; i < WARMUP; i++) { hipblasDotEx_64(h, size*2, (half*)d_a, xType, 1, (half*)d_b, yType, 1, d_r2, rType, eType); } @@ -157,10 +172,38 @@ int main(int argc, char *argv[]) hipDeviceSynchronize(); auto end = std::chrono::steady_clock::now(); auto time = std::chrono::duration_cast(end - start).count(); - printf("Average kernel (hipblasDot) execution time %f (us)\n", (time * 1e-3f) / repeat); + printf("Average hipblasDotEx_64 execution time %f (us)\n", (time * 1e-3f) / repeat); hipMemcpy(&r2, d_r2, result2_bytes, hipMemcpyDeviceToHost); printf("Error rate: %e\n", fabsf((float)r2 - 65504.f)/65504.f); + printf("\n"); + + for (int i = 0; i < WARMUP + repeat; i++) { + if (i == WARMUP) { + hipDeviceSynchronize(); + start = std::chrono::steady_clock::now(); + } + hipcub::CountingInputIterator counting(0); + hipcub::TransformInputIterator> + iter(counting, Half2Mul(d_a, d_b)); + + void* d_temp = nullptr; size_t temp_bytes = 0; + hipcub::DeviceReduce::Sum(d_temp, temp_bytes, iter, d_r, size); + if (temp_bytes != 0) hipMalloc(&d_temp, temp_bytes); + auto error = hipcub::DeviceReduce::Sum(d_temp, temp_bytes, iter, d_r, size); + if (error != hipSuccess) + printf("hipcub::DeviceReduce::Sum failed\n"); + if (temp_bytes != 0) hipFree(d_temp); + } + hipDeviceSynchronize(); + end = std::chrono::steady_clock::now(); + time = std::chrono::duration_cast(end - start).count(); + printf("Average hipcub::DeviceReduce::Sum execution time: %f (us)\n", time * 1e-3f / repeat); + + hipMemcpy(&r, d_r, result_bytes, hipMemcpyDeviceToHost); + printf("Error rate: %e\n", fabsf(r - 65504.f)/65504.f); + hipFree(d_a); hipFree(d_b); hipFree(d_r); diff --git a/src/f16sp-omp/main.cpp b/src/f16sp-omp/main.cpp index a657435c73..65163e7c50 100644 --- a/src/f16sp-omp/main.cpp +++ b/src/f16sp-omp/main.cpp @@ -36,11 +36,11 @@ void kernel1(const __half* __restrict__ ax, #pragma omp parallel num_threads(NUM_OF_THREADS) { int lid = omp_get_thread_num(); - int stride = grid * NUM_OF_THREADS; - int base = omp_get_team_num() * NUM_OF_THREADS + lid; + size_t stride = (size_t)grid * NUM_OF_THREADS; + size_t base = (size_t)omp_get_team_num() * NUM_OF_THREADS + lid; __half val_x = 0, val_y = 0; - for (size_t i = base; (size_t)i < size; i += stride) { + for (size_t i = base; i < size; i += stride) { val_x += ax[i] * bx[i]; val_y += ay[i] * by[i]; } @@ -82,11 +82,11 @@ void kernel2(const __half* __restrict__ ax, #pragma omp parallel num_threads(NUM_OF_THREADS) { int lid = omp_get_thread_num(); - int stride = grid * NUM_OF_THREADS; - int base = omp_get_team_num() * NUM_OF_THREADS + lid; + size_t stride = (size_t)grid * NUM_OF_THREADS; + size_t base = (size_t)omp_get_team_num() * NUM_OF_THREADS + lid; float val_x = 0.f, val_y = 0.f; - for (int i = base; (size_t)i < size; i += stride) { + for (size_t i = base; i < size; i += stride) { val_x += float(ax[i]) * float(bx[i]); val_y += float(ay[i]) * float(by[i]); } diff --git a/src/f16sp-sycl/kernels.h b/src/f16sp-sycl/kernels.h index 393325172f..2bf54d4019 100644 --- a/src/f16sp-sycl/kernels.h +++ b/src/f16sp-sycl/kernels.h @@ -33,12 +33,12 @@ void scalarProductKernel_native(const sycl::half2 *a, { int lid = item.get_local_id(0); - const int stride = item.get_group_range(0) * item.get_local_range(0); + const size_t stride = item.get_group_range(0) * item.get_local_range(0); sycl::half2 value(0.f, 0.f); shArray[lid] = value; - for (int i = item.get_global_id(0); i < size; i += stride) + for (size_t i = item.get_global_id(0); i < size; i += stride) { value = a[i] * b[i] + value; } @@ -67,12 +67,12 @@ void scalarProductKernel_native_fp32(const sycl::half2 *a, { int lid = item.get_local_id(0); - const int stride = item.get_group_range(0) * item.get_local_range(0); + const size_t stride = item.get_group_range(0) * item.get_local_range(0); sycl::float2 value(0.f, 0.f); shArray[lid] = value; - for (int i = item.get_global_id(0); i < size; i += stride) + for (size_t i = item.get_global_id(0); i < size; i += stride) { value.x() += (float)a[i].x() * (float)b[i].x(); value.y() += (float)a[i].y() * (float)b[i].y(); @@ -100,11 +100,11 @@ void scalarProductKernel_native2_fp32(const sycl::half2 *a, { int lid = item.get_local_id(0); - const int stride = item.get_group_range(0) * item.get_local_range(0); + const size_t stride = item.get_group_range(0) * item.get_local_range(0); sycl::float2 value(0.f, 0.f); - for (int i = item.get_global_id(0); i < size; i += stride) + for (size_t i = item.get_global_id(0); i < size; i += stride) { value.x() += (float)a[i].x() * (float)b[i].x(); value.y() += (float)a[i].y() * (float)b[i].y(); diff --git a/src/f16sp-sycl/main.cpp b/src/f16sp-sycl/main.cpp index 7bc124d57a..e8337d4156 100644 --- a/src/f16sp-sycl/main.cpp +++ b/src/f16sp-sycl/main.cpp @@ -78,7 +78,7 @@ int main(int argc, char *argv[]) auto fp16Fn = [&](sycl::handler &cgh) { sycl::local_accessor sm(sycl::range<1>(NUM_OF_THREADS), cgh); cgh.parallel_for(sycl::nd_range<1>(sycl::range<1>(grid*NUM_OF_THREADS), - sycl::range<1>(NUM_OF_THREADS)), [=](sycl::nd_item<1> item) { + sycl::range<1>(NUM_OF_THREADS)), [[sycl::reqd_work_group_size(NUM_OF_THREADS)]] [=](sycl::nd_item<1> item) { scalarProductKernel_native(d_a, d_b, d_r, sm.get_multi_ptr().get(), size, item); @@ -87,7 +87,7 @@ int main(int argc, char *argv[]) auto fp32Fn = [&](sycl::handler &cgh) { sycl::local_accessor sm(sycl::range<1>(NUM_OF_THREADS), cgh); cgh.parallel_for(sycl::nd_range<1>(sycl::range<1>(grid*NUM_OF_THREADS), - sycl::range<1>(NUM_OF_THREADS)), [=](sycl::nd_item<1> item) { + sycl::range<1>(NUM_OF_THREADS)), [[sycl::reqd_work_group_size(NUM_OF_THREADS)]] [=](sycl::nd_item<1> item) { scalarProductKernel_native_fp32(d_a, d_b, d_r, sm.get_multi_ptr().get(), size, item); @@ -142,7 +142,7 @@ int main(int argc, char *argv[]) cgh.parallel_for( sycl::nd_range<1>(sycl::range<1>(grid*NUM_OF_THREADS), sycl::range<1>(NUM_OF_THREADS)), - [=](sycl::nd_item<1> item) { + [[sycl::reqd_work_group_size(NUM_OF_THREADS)]] [=](sycl::nd_item<1> item) { scalarProductKernel_native2_fp32(d_a, d_b, d_r, size, item); }); }); @@ -157,7 +157,7 @@ int main(int argc, char *argv[]) cgh.parallel_for( sycl::nd_range<1>(sycl::range<1>(grid*NUM_OF_THREADS), sycl::range<1>(NUM_OF_THREADS)), - [=](sycl::nd_item<1> item) { + [[sycl::reqd_work_group_size(NUM_OF_THREADS)]] [=](sycl::nd_item<1> item) { scalarProductKernel_native2_fp32(d_a, d_b, d_r, size, item); }); }); diff --git a/src/mmcsf-cuda/mttkrp_gpu.h b/src/mmcsf-cuda/mttkrp_gpu.h index 8b22602bb8..dd664a4149 100644 --- a/src/mmcsf-cuda/mttkrp_gpu.h +++ b/src/mmcsf-cuda/mttkrp_gpu.h @@ -2,7 +2,11 @@ #define MTTKRP_GPU_H #include "util.h" -#include +#if defined(__HIPCC__) + #include +#else + #include +#endif int MTTKRP_MIHCSR_GPU(TiledTensor *TiledX, Matrix *U, const Options &Opt); diff --git a/src/opticalFlow-cuda/Common/helper_image.h b/src/opticalFlow-cuda/Common/helper_image.h index 9b7edc062c..f35f482cc0 100644 --- a/src/opticalFlow-cuda/Common/helper_image.h +++ b/src/opticalFlow-cuda/Common/helper_image.h @@ -530,7 +530,7 @@ inline bool compareData(const T *reference, const T *data, error_count += !comp; -#if 0 +#ifdef DEBUG if (!comp) { std::cerr << "ERROR, i = " << i << ",\t " diff --git a/src/opticalFlow-cuda/common.h b/src/opticalFlow-cuda/common.h index 028aa80341..2d6597b21e 100644 --- a/src/opticalFlow-cuda/common.h +++ b/src/opticalFlow-cuda/common.h @@ -47,7 +47,7 @@ /////////////////////////////////////////////////////////////////////////////// // Common constants /////////////////////////////////////////////////////////////////////////////// -const int StrideAlignment = 32; +const int StrideAlignment = 256; /////////////////////////////////////////////////////////////////////////////// // Common functions diff --git a/src/opticalFlow-cuda/derivativesKernel.cuh b/src/opticalFlow-cuda/derivativesKernel.cuh index f26c56ee8f..b862db9562 100644 --- a/src/opticalFlow-cuda/derivativesKernel.cuh +++ b/src/opticalFlow-cuda/derivativesKernel.cuh @@ -1,148 +1,70 @@ -/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * * Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * * Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * * Neither the name of NVIDIA CORPORATION nor the names of its - * contributors may be used to endorse or promote products derived - * from this software without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY - * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR - * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY - * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) */ - #include "common.h" -/////////////////////////////////////////////////////////////////////////////// -/// \brief compute image derivatives -/// -/// CUDA kernel, relies heavily on texture unit -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[out] Ix x derivative -/// \param[out] Iy y derivative -/// \param[out] Iz temporal derivative -/////////////////////////////////////////////////////////////////////////////// +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(t,x,y) tex2D_sw(t,x,y) +#define TEX_CREATE(v,s,w,h,st) Tex2DDesc v{s,w,h,st} +#define TEX_DESTROY(v) + __global__ void ComputeDerivativesKernel(int width, int height, int stride, float *Ix, float *Iy, float *Iz, - cudaTextureObject_t texSource, - cudaTextureObject_t texTarget) { - const int ix = threadIdx.x + blockIdx.x * blockDim.x; - const int iy = threadIdx.y + blockIdx.y * blockDim.y; - + TEX_TYPE texSource, TEX_TYPE texTarget) { + const int ix = threadIdx.x + blockIdx.x * blockDim.x; + const int iy = threadIdx.y + blockIdx.y * blockDim.y; const int pos = ix + iy * stride; - if (ix >= width || iy >= height) return; - float dx = 1.0f / (float)width; - float dy = 1.0f / (float)height; - - float x = ((float)ix + 0.5f) * dx; - float y = ((float)iy + 0.5f) * dy; + float dx = 1.0f/(float)width, dy = 1.0f/(float)height; + float x = ((float)ix + 0.5f) * dx; + float y = ((float)iy + 0.5f) * dy; float t0, t1; - // x derivative - t0 = tex2D(texSource, x - 2.0f * dx, y); - t0 -= tex2D(texSource, x - 1.0f * dx, y) * 8.0f; - t0 += tex2D(texSource, x + 1.0f * dx, y) * 8.0f; - t0 -= tex2D(texSource, x + 2.0f * dx, y); - t0 /= 12.0f; - - t1 = tex2D(texTarget, x - 2.0f * dx, y); - t1 -= tex2D(texTarget, x - 1.0f * dx, y) * 8.0f; - t1 += tex2D(texTarget, x + 1.0f * dx, y) * 8.0f; - t1 -= tex2D(texTarget, x + 2.0f * dx, y); - t1 /= 12.0f; - + t0 = TEX_FETCH(texSource, x-2.f*dx, y); + t0 -= TEX_FETCH(texSource, x-1.f*dx, y) * 8.f; + t0 += TEX_FETCH(texSource, x+1.f*dx, y) * 8.f; + t0 -= TEX_FETCH(texSource, x+2.f*dx, y); + t0 /= 12.f; + + t1 = TEX_FETCH(texTarget, x-2.f*dx, y); + t1 -= TEX_FETCH(texTarget, x-1.f*dx, y) * 8.f; + t1 += TEX_FETCH(texTarget, x+1.f*dx, y) * 8.f; + t1 -= TEX_FETCH(texTarget, x+2.f*dx, y); + t1 /= 12.f; Ix[pos] = (t0 + t1) * 0.5f; - // t derivative - Iz[pos] = tex2D(texTarget, x, y) - tex2D(texSource, x, y); - - // y derivative - t0 = tex2D(texSource, x, y - 2.0f * dy); - t0 -= tex2D(texSource, x, y - 1.0f * dy) * 8.0f; - t0 += tex2D(texSource, x, y + 1.0f * dy) * 8.0f; - t0 -= tex2D(texSource, x, y + 2.0f * dy); - t0 /= 12.0f; + Iz[pos] = TEX_FETCH(texTarget, x, y) - TEX_FETCH(texSource, x, y); - t1 = tex2D(texTarget, x, y - 2.0f * dy); - t1 -= tex2D(texTarget, x, y - 1.0f * dy) * 8.0f; - t1 += tex2D(texTarget, x, y + 1.0f * dy) * 8.0f; - t1 -= tex2D(texTarget, x, y + 2.0f * dy); - t1 /= 12.0f; + t0 = TEX_FETCH(texSource, x, y-2.f*dy); + t0 -= TEX_FETCH(texSource, x, y-1.f*dy) * 8.f; + t0 += TEX_FETCH(texSource, x, y+1.f*dy) * 8.f; + t0 -= TEX_FETCH(texSource, x, y+2.f*dy); + t0 /= 12.f; + t1 = TEX_FETCH(texTarget, x, y-2.f*dy); + t1 -= TEX_FETCH(texTarget, x, y-1.f*dy) * 8.f; + t1 += TEX_FETCH(texTarget, x, y+1.f*dy) * 8.f; + t1 -= TEX_FETCH(texTarget, x, y+2.f*dy); + t1 /= 12.f; Iy[pos] = (t0 + t1) * 0.5f; } -/////////////////////////////////////////////////////////////////////////////// -/// \brief compute image derivatives -/// -/// \param[in] I0 source image -/// \param[in] I1 tracked image -/// \param[in] w image width -/// \param[in] h image height -/// \param[in] s image stride -/// \param[out] Ix x derivative -/// \param[out] Iy y derivative -/// \param[out] Iz temporal derivative -/////////////////////////////////////////////////////////////////////////////// -static void ComputeDerivatives(const float *I0, const float *I1, int w, int h, - int s, float *Ix, float *Iy, float *Iz) { +static void ComputeDerivatives(const float *I0, const float *I1, + int w, int h, int s, + float *Ix, float *Iy, float *Iz) { dim3 threads(32, 6); dim3 blocks(iDivUp(w, threads.x), iDivUp(h, threads.y)); - - cudaTextureObject_t texSource, texTarget; - cudaResourceDesc texRes; - memset(&texRes, 0, sizeof(cudaResourceDesc)); - - texRes.resType = cudaResourceTypePitch2D; - texRes.res.pitch2D.devPtr = (void *)I0; - texRes.res.pitch2D.desc = cudaCreateChannelDesc(); - texRes.res.pitch2D.width = w; - texRes.res.pitch2D.height = h; - texRes.res.pitch2D.pitchInBytes = s * sizeof(float); - - cudaTextureDesc texDescr; - memset(&texDescr, 0, sizeof(cudaTextureDesc)); - - texDescr.normalizedCoords = true; - texDescr.filterMode = cudaFilterModeLinear; - texDescr.addressMode[0] = cudaAddressModeMirror; - texDescr.addressMode[1] = cudaAddressModeMirror; - texDescr.readMode = cudaReadModeElementType; - - checkCudaErrors( - cudaCreateTextureObject(&texSource, &texRes, &texDescr, NULL)); - memset(&texRes, 0, sizeof(cudaResourceDesc)); - texRes.resType = cudaResourceTypePitch2D; - texRes.res.pitch2D.devPtr = (void *)I1; - texRes.res.pitch2D.desc = cudaCreateChannelDesc(); - texRes.res.pitch2D.width = w; - texRes.res.pitch2D.height = h; - texRes.res.pitch2D.pitchInBytes = s * sizeof(float); - checkCudaErrors( - cudaCreateTextureObject(&texTarget, &texRes, &texDescr, NULL)); - - ComputeDerivativesKernel<<>>(w, h, s, Ix, Iy, Iz, texSource, - texTarget); - + TEX_CREATE(texSource, I0, w, h, s); + TEX_CREATE(texTarget, I1, w, h, s); + ComputeDerivativesKernel<<>>(w, h, s, Ix, Iy, Iz, + texSource, texTarget); checkCudaErrors(cudaDeviceSynchronize()); - checkCudaErrors(cudaDestroyTextureObject(texSource)); - checkCudaErrors(cudaDestroyTextureObject(texTarget)); + TEX_DESTROY(texSource); + TEX_DESTROY(texTarget); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-cuda/downscaleKernel.cuh b/src/opticalFlow-cuda/downscaleKernel.cuh index 2c0d52844f..57f011e15f 100644 --- a/src/opticalFlow-cuda/downscaleKernel.cuh +++ b/src/opticalFlow-cuda/downscaleKernel.cuh @@ -1,101 +1,42 @@ -/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * * Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * * Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * * Neither the name of NVIDIA CORPORATION nor the names of its - * contributors may be used to endorse or promote products derived - * from this software without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY - * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR - * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY - * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) */ - #include "common.h" -/////////////////////////////////////////////////////////////////////////////// -/// \brief downscale image -/// -/// CUDA kernel, relies heavily on texture unit -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// -__global__ void DownscaleKernel(int width, int height, int stride, float *out, - cudaTextureObject_t texFine) { +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(tex, x, y) tex2D_sw(tex, x, y) +#define TEX_CREATE(texVar, src, width, height, stride) \ + Tex2DDesc texVar { src, width, height, stride } +#define TEX_DESTROY(texVar) /* no-op: no object to destroy */ + +__global__ void DownscaleKernel(int width, int height, int stride, + float *out, TEX_TYPE texFine) { const int ix = threadIdx.x + blockIdx.x * blockDim.x; const int iy = threadIdx.y + blockIdx.y * blockDim.y; - - if (ix >= width || iy >= height) { - return; - } + if (ix >= width || iy >= height) return; float dx = 1.0f / (float)width; float dy = 1.0f / (float)height; + float x = ((float)ix + 0.5f) * dx; + float y = ((float)iy + 0.5f) * dy; - float x = ((float)ix + 0.5f) * dx; - float y = ((float)iy + 0.5f) * dy; - - out[ix + iy * stride] = 0.25f * (tex2D(texFine, x - dx * 0.25f, y) + - tex2D(texFine, x + dx * 0.25f, y) + - tex2D(texFine, x, y - dy * 0.25f) + - tex2D(texFine, x, y + dy * 0.25f)); + out[ix + iy * stride] = 0.25f * (TEX_FETCH(texFine, x - dx*0.25f, y) + + TEX_FETCH(texFine, x + dx*0.25f, y) + + TEX_FETCH(texFine, x, y - dy*0.25f) + + TEX_FETCH(texFine, x, y + dy*0.25f)); } -/////////////////////////////////////////////////////////////////////////////// -/// \brief downscale image -/// -/// \param[in] src image to downscale -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// static void Downscale(const float *src, int width, int height, int stride, int newWidth, int newHeight, int newStride, float *out) { dim3 threads(32, 8); dim3 blocks(iDivUp(newWidth, threads.x), iDivUp(newHeight, threads.y)); - - cudaTextureObject_t texFine; - cudaResourceDesc texRes; - memset(&texRes, 0, sizeof(cudaResourceDesc)); - - texRes.resType = cudaResourceTypePitch2D; - texRes.res.pitch2D.devPtr = (void *)src; - texRes.res.pitch2D.desc = cudaCreateChannelDesc(); - texRes.res.pitch2D.width = width; - texRes.res.pitch2D.height = height; - texRes.res.pitch2D.pitchInBytes = stride * sizeof(float); - - cudaTextureDesc texDescr; - memset(&texDescr, 0, sizeof(cudaTextureDesc)); - - texDescr.normalizedCoords = true; - texDescr.filterMode = cudaFilterModeLinear; - texDescr.addressMode[0] = cudaAddressModeMirror; - texDescr.addressMode[1] = cudaAddressModeMirror; - texDescr.readMode = cudaReadModeElementType; - - checkCudaErrors(cudaCreateTextureObject(&texFine, &texRes, &texDescr, NULL)); - - DownscaleKernel<<>>(newWidth, newHeight, newStride, out, - texFine); - + TEX_CREATE(texFine, src, width, height, stride); + DownscaleKernel<<>>(newWidth, newHeight, newStride, out, texFine); checkCudaErrors(cudaDeviceSynchronize()); - checkCudaErrors(cudaDestroyTextureObject(texFine)); + TEX_DESTROY(texFine); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-cuda/tex2d_compat.h b/src/opticalFlow-cuda/tex2d_compat.h new file mode 100644 index 0000000000..efea3cb226 --- /dev/null +++ b/src/opticalFlow-cuda/tex2d_compat.h @@ -0,0 +1,41 @@ +// tex2d_compat.h — Portable tex2D + +#pragma once + +// Descriptor passed to kernels in place of hipTextureObject_t +struct Tex2DDesc { + const float* ptr; // pitched device pointer + int width; // image width in elements + int height; // image height in elements + int pitch; // row stride in elements (pitchInBytes / sizeof(float)) +}; + +// Mirror repeat address clamping +__device__ inline int tex2d_mirror_clamp(int i, int size) { + int mi = ((i % (2 * size)) + 2 * size) % (2 * size); + return (mi >= size) ? (2 * size - 1 - mi) : mi; +} + +// Software bilinear tex2D — normalized coords [0,1), mirror addressing, linear filter +// Matches: normalizedCoords=true, filterMode=Linear, addressMode=Mirror +__device__ inline float tex2D_sw(const Tex2DDesc& tex, float u, float v) { + // Convert normalized coords to pixel space + float fx = u * (float)tex.width - 0.5f; + float fy = v * (float)tex.height - 0.5f; + + int x0 = (int)floorf(fx), y0 = (int)floorf(fy); + float tx = fx - (float)x0, ty = fy - (float)y0; + + int cx0 = tex2d_mirror_clamp(x0, tex.width); + int cx1 = tex2d_mirror_clamp(x0+1, tex.width); + int cy0 = tex2d_mirror_clamp(y0, tex.height); + int cy1 = tex2d_mirror_clamp(y0+1, tex.height); + + float v00 = tex.ptr[cy0 * tex.pitch + cx0]; + float v10 = tex.ptr[cy0 * tex.pitch + cx1]; + float v01 = tex.ptr[cy1 * tex.pitch + cx0]; + float v11 = tex.ptr[cy1 * tex.pitch + cx1]; + + return (1.f-tx)*(1.f-ty)*v00 + tx*(1.f-ty)*v10 + + (1.f-tx)*ty *v01 + tx*ty *v11; +} diff --git a/src/opticalFlow-cuda/upscaleKernel.cuh b/src/opticalFlow-cuda/upscaleKernel.cuh index 838efbbe31..d9ec9d1cd8 100644 --- a/src/opticalFlow-cuda/upscaleKernel.cuh +++ b/src/opticalFlow-cuda/upscaleKernel.cuh @@ -1,100 +1,34 @@ -/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * * Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * * Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * * Neither the name of NVIDIA CORPORATION nor the names of its - * contributors may be used to endorse or promote products derived - * from this software without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY - * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR - * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY - * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) */ - #include "common.h" - -/////////////////////////////////////////////////////////////////////////////// -/// \brief upscale one component of a displacement field, CUDA kernel -/// \param[in] width field width -/// \param[in] height field height -/// \param[in] stride field stride -/// \param[in] scale scale factor (multiplier) -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// -__global__ void UpscaleKernel(int width, int height, int stride, float scale, - float *out, cudaTextureObject_t texCoarse) { +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(t,x,y) tex2D_sw(t,x,y) +#define TEX_CREATE(v,s,w,h,st) Tex2DDesc v{s,w,h,st} +#define TEX_DESTROY(v) + +__global__ void UpscaleKernel(int width, int height, int stride, + float scale, float *out, TEX_TYPE texCoarse) { const int ix = threadIdx.x + blockIdx.x * blockDim.x; const int iy = threadIdx.y + blockIdx.y * blockDim.y; - if (ix >= width || iy >= height) return; - float x = ((float)ix + 0.5f) / (float)width; float y = ((float)iy + 0.5f) / (float)height; - - // exploit hardware interpolation - // and scale interpolated vector to match next pyramid level resolution - out[ix + iy * stride] = tex2D(texCoarse, x, y) * scale; + out[ix + iy * stride] = TEX_FETCH(texCoarse, x, y) * scale; } -/////////////////////////////////////////////////////////////////////////////// -/// \brief upscale one component of a displacement field, kernel wrapper -/// \param[in] src field component to upscale -/// \param[in] width field current width -/// \param[in] height field current height -/// \param[in] stride field current stride -/// \param[in] newWidth field new width -/// \param[in] newHeight field new height -/// \param[in] newStride field new stride -/// \param[in] scale value scale factor (multiplier) -/// \param[out] out upscaled field component -/////////////////////////////////////////////////////////////////////////////// static void Upscale(const float *src, int width, int height, int stride, - int newWidth, int newHeight, int newStride, float scale, - float *out) { + int newWidth, int newHeight, int newStride, + float scale, float *out) { dim3 threads(32, 8); dim3 blocks(iDivUp(newWidth, threads.x), iDivUp(newHeight, threads.y)); - - cudaTextureObject_t texCoarse; - cudaResourceDesc texRes; - memset(&texRes, 0, sizeof(cudaResourceDesc)); - - texRes.resType = cudaResourceTypePitch2D; - texRes.res.pitch2D.devPtr = (void *)src; - texRes.res.pitch2D.desc = cudaCreateChannelDesc(); - texRes.res.pitch2D.width = width; - texRes.res.pitch2D.height = height; - texRes.res.pitch2D.pitchInBytes = stride * sizeof(float); - - cudaTextureDesc texDescr; - memset(&texDescr, 0, sizeof(cudaTextureDesc)); - - texDescr.normalizedCoords = true; - texDescr.filterMode = cudaFilterModeLinear; - texDescr.addressMode[0] = cudaAddressModeMirror; - texDescr.addressMode[1] = cudaAddressModeMirror; - texDescr.readMode = cudaReadModeElementType; - - checkCudaErrors( - cudaCreateTextureObject(&texCoarse, &texRes, &texDescr, NULL)); - - UpscaleKernel<<>>(newWidth, newHeight, newStride, scale, out, - texCoarse); - + TEX_CREATE(texCoarse, src, width, height, stride); + UpscaleKernel<<>>(newWidth, newHeight, newStride, scale, out, texCoarse); checkCudaErrors(cudaDeviceSynchronize()); - checkCudaErrors(cudaDestroyTextureObject(texCoarse)); - + TEX_DESTROY(texCoarse); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-cuda/warpingKernel.cuh b/src/opticalFlow-cuda/warpingKernel.cuh index 2ef86337a2..7c9037a422 100644 --- a/src/opticalFlow-cuda/warpingKernel.cuh +++ b/src/opticalFlow-cuda/warpingKernel.cuh @@ -1,103 +1,35 @@ -/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * * Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * * Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * * Neither the name of NVIDIA CORPORATION nor the names of its - * contributors may be used to endorse or promote products derived - * from this software without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY - * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR - * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY - * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) */ - #include "common.h" - -/////////////////////////////////////////////////////////////////////////////// -/// \brief warp image with a given displacement field, CUDA kernel. -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[in] u horizontal displacement -/// \param[in] v vertical displacement -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// -__global__ void WarpingKernel(int width, int height, int stride, const float *u, - const float *v, float *out, - cudaTextureObject_t texToWarp) { - const int ix = threadIdx.x + blockIdx.x * blockDim.x; - const int iy = threadIdx.y + blockIdx.y * blockDim.y; - +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(t,x,y) tex2D_sw(t,x,y) +#define TEX_CREATE(v,s,w,h,st) Tex2DDesc v{s,w,h,st} +#define TEX_DESTROY(v) + +__global__ void WarpingKernel(int width, int height, int stride, + const float *u, const float *v, + float *out, TEX_TYPE texToWarp) { + const int ix = threadIdx.x + blockIdx.x * blockDim.x; + const int iy = threadIdx.y + blockIdx.y * blockDim.y; const int pos = ix + iy * stride; - if (ix >= width || iy >= height) return; - float x = ((float)ix + u[pos] + 0.5f) / (float)width; float y = ((float)iy + v[pos] + 0.5f) / (float)height; - - out[pos] = tex2D(texToWarp, x, y); + out[pos] = TEX_FETCH(texToWarp, x, y); } -/////////////////////////////////////////////////////////////////////////////// -/// \brief warp image with provided vector field, CUDA kernel wrapper. -/// -/// For each output pixel there is a vector which tells which pixel -/// from a source image should be mapped to this particular output -/// pixel. -/// It is assumed that images and the vector field have the same stride and -/// resolution. -/// \param[in] src source image -/// \param[in] w width -/// \param[in] h height -/// \param[in] s stride -/// \param[in] u horizontal displacement -/// \param[in] v vertical displacement -/// \param[out] out warped image -/////////////////////////////////////////////////////////////////////////////// -static void WarpImage(const float *src, int w, int h, int s, const float *u, - const float *v, float *out) { +static void WarpImage(const float *src, int w, int h, int s, + const float *u, const float *v, float *out) { dim3 threads(32, 6); dim3 blocks(iDivUp(w, threads.x), iDivUp(h, threads.y)); - - cudaTextureObject_t texToWarp; - cudaResourceDesc texRes; - memset(&texRes, 0, sizeof(cudaResourceDesc)); - - texRes.resType = cudaResourceTypePitch2D; - texRes.res.pitch2D.devPtr = (void *)src; - texRes.res.pitch2D.desc = cudaCreateChannelDesc(); - texRes.res.pitch2D.width = w; - texRes.res.pitch2D.height = h; - texRes.res.pitch2D.pitchInBytes = s * sizeof(float); - - cudaTextureDesc texDescr; - memset(&texDescr, 0, sizeof(cudaTextureDesc)); - - texDescr.normalizedCoords = true; - texDescr.filterMode = cudaFilterModeLinear; - texDescr.addressMode[0] = cudaAddressModeMirror; - texDescr.addressMode[1] = cudaAddressModeMirror; - texDescr.readMode = cudaReadModeElementType; - - checkCudaErrors( - cudaCreateTextureObject(&texToWarp, &texRes, &texDescr, NULL)); - + TEX_CREATE(texToWarp, src, w, h, s); WarpingKernel<<>>(w, h, s, u, v, out, texToWarp); - checkCudaErrors(cudaDeviceSynchronize()); - checkCudaErrors(cudaDestroyTextureObject(texToWarp)); + TEX_DESTROY(texToWarp); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-hip/derivativesKernel.cuh b/src/opticalFlow-hip/derivativesKernel.cuh index f96a033f79..44723b5319 100644 --- a/src/opticalFlow-hip/derivativesKernel.cuh +++ b/src/opticalFlow-hip/derivativesKernel.cuh @@ -1,148 +1,70 @@ -/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * * Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * * Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * * Neither the name of NVIDIA CORPORATION nor the names of its - * contributors may be used to endorse or promote products derived - * from this software without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY - * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR - * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY - * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) */ - #include "common.h" -/////////////////////////////////////////////////////////////////////////////// -/// \brief compute image derivatives -/// -/// CUDA kernel, relies heavily on texture unit -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[out] Ix x derivative -/// \param[out] Iy y derivative -/// \param[out] Iz temporal derivative -/////////////////////////////////////////////////////////////////////////////// +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(t,x,y) tex2D_sw(t,x,y) +#define TEX_CREATE(v,s,w,h,st) Tex2DDesc v{s,w,h,st} +#define TEX_DESTROY(v) + __global__ void ComputeDerivativesKernel(int width, int height, int stride, float *Ix, float *Iy, float *Iz, - hipTextureObject_t texSource, - hipTextureObject_t texTarget) { - const int ix = threadIdx.x + blockIdx.x * blockDim.x; - const int iy = threadIdx.y + blockIdx.y * blockDim.y; - + TEX_TYPE texSource, TEX_TYPE texTarget) { + const int ix = threadIdx.x + blockIdx.x * blockDim.x; + const int iy = threadIdx.y + blockIdx.y * blockDim.y; const int pos = ix + iy * stride; - if (ix >= width || iy >= height) return; - float dx = 1.0f / (float)width; - float dy = 1.0f / (float)height; - - float x = ((float)ix + 0.5f) * dx; - float y = ((float)iy + 0.5f) * dy; + float dx = 1.0f/(float)width, dy = 1.0f/(float)height; + float x = ((float)ix + 0.5f) * dx; + float y = ((float)iy + 0.5f) * dy; float t0, t1; - // x derivative - t0 = tex2D(texSource, x - 2.0f * dx, y); - t0 -= tex2D(texSource, x - 1.0f * dx, y) * 8.0f; - t0 += tex2D(texSource, x + 1.0f * dx, y) * 8.0f; - t0 -= tex2D(texSource, x + 2.0f * dx, y); - t0 /= 12.0f; - - t1 = tex2D(texTarget, x - 2.0f * dx, y); - t1 -= tex2D(texTarget, x - 1.0f * dx, y) * 8.0f; - t1 += tex2D(texTarget, x + 1.0f * dx, y) * 8.0f; - t1 -= tex2D(texTarget, x + 2.0f * dx, y); - t1 /= 12.0f; - + t0 = TEX_FETCH(texSource, x-2.f*dx, y); + t0 -= TEX_FETCH(texSource, x-1.f*dx, y) * 8.f; + t0 += TEX_FETCH(texSource, x+1.f*dx, y) * 8.f; + t0 -= TEX_FETCH(texSource, x+2.f*dx, y); + t0 /= 12.f; + + t1 = TEX_FETCH(texTarget, x-2.f*dx, y); + t1 -= TEX_FETCH(texTarget, x-1.f*dx, y) * 8.f; + t1 += TEX_FETCH(texTarget, x+1.f*dx, y) * 8.f; + t1 -= TEX_FETCH(texTarget, x+2.f*dx, y); + t1 /= 12.f; Ix[pos] = (t0 + t1) * 0.5f; - // t derivative - Iz[pos] = tex2D(texTarget, x, y) - tex2D(texSource, x, y); - - // y derivative - t0 = tex2D(texSource, x, y - 2.0f * dy); - t0 -= tex2D(texSource, x, y - 1.0f * dy) * 8.0f; - t0 += tex2D(texSource, x, y + 1.0f * dy) * 8.0f; - t0 -= tex2D(texSource, x, y + 2.0f * dy); - t0 /= 12.0f; + Iz[pos] = TEX_FETCH(texTarget, x, y) - TEX_FETCH(texSource, x, y); - t1 = tex2D(texTarget, x, y - 2.0f * dy); - t1 -= tex2D(texTarget, x, y - 1.0f * dy) * 8.0f; - t1 += tex2D(texTarget, x, y + 1.0f * dy) * 8.0f; - t1 -= tex2D(texTarget, x, y + 2.0f * dy); - t1 /= 12.0f; + t0 = TEX_FETCH(texSource, x, y-2.f*dy); + t0 -= TEX_FETCH(texSource, x, y-1.f*dy) * 8.f; + t0 += TEX_FETCH(texSource, x, y+1.f*dy) * 8.f; + t0 -= TEX_FETCH(texSource, x, y+2.f*dy); + t0 /= 12.f; + t1 = TEX_FETCH(texTarget, x, y-2.f*dy); + t1 -= TEX_FETCH(texTarget, x, y-1.f*dy) * 8.f; + t1 += TEX_FETCH(texTarget, x, y+1.f*dy) * 8.f; + t1 -= TEX_FETCH(texTarget, x, y+2.f*dy); + t1 /= 12.f; Iy[pos] = (t0 + t1) * 0.5f; } -/////////////////////////////////////////////////////////////////////////////// -/// \brief compute image derivatives -/// -/// \param[in] I0 source image -/// \param[in] I1 tracked image -/// \param[in] w image width -/// \param[in] h image height -/// \param[in] s image stride -/// \param[out] Ix x derivative -/// \param[out] Iy y derivative -/// \param[out] Iz temporal derivative -/////////////////////////////////////////////////////////////////////////////// -static void ComputeDerivatives(const float *I0, const float *I1, int w, int h, - int s, float *Ix, float *Iy, float *Iz) { +static void ComputeDerivatives(const float *I0, const float *I1, + int w, int h, int s, + float *Ix, float *Iy, float *Iz) { dim3 threads(32, 6); dim3 blocks(iDivUp(w, threads.x), iDivUp(h, threads.y)); - - hipTextureObject_t texSource, texTarget; - hipResourceDesc texRes; - memset(&texRes, 0, sizeof(hipResourceDesc)); - - texRes.resType = hipResourceTypePitch2D; - texRes.res.pitch2D.devPtr = (void *)I0; - texRes.res.pitch2D.desc = hipCreateChannelDesc(); - texRes.res.pitch2D.width = w; - texRes.res.pitch2D.height = h; - texRes.res.pitch2D.pitchInBytes = s * sizeof(float); - - hipTextureDesc texDescr; - memset(&texDescr, 0, sizeof(hipTextureDesc)); - - texDescr.normalizedCoords = true; - texDescr.filterMode = hipFilterModeLinear; - texDescr.addressMode[0] = hipAddressModeMirror; - texDescr.addressMode[1] = hipAddressModeMirror; - texDescr.readMode = hipReadModeElementType; - - checkCudaErrors( - hipCreateTextureObject(&texSource, &texRes, &texDescr, NULL)); - memset(&texRes, 0, sizeof(hipResourceDesc)); - texRes.resType = hipResourceTypePitch2D; - texRes.res.pitch2D.devPtr = (void *)I1; - texRes.res.pitch2D.desc = hipCreateChannelDesc(); - texRes.res.pitch2D.width = w; - texRes.res.pitch2D.height = h; - texRes.res.pitch2D.pitchInBytes = s * sizeof(float); - checkCudaErrors( - hipCreateTextureObject(&texTarget, &texRes, &texDescr, NULL)); - - ComputeDerivativesKernel<<>>(w, h, s, Ix, Iy, Iz, texSource, - texTarget); - + TEX_CREATE(texSource, I0, w, h, s); + TEX_CREATE(texTarget, I1, w, h, s); + ComputeDerivativesKernel<<>>(w, h, s, Ix, Iy, Iz, + texSource, texTarget); checkCudaErrors(hipDeviceSynchronize()); - checkCudaErrors(hipDestroyTextureObject(texSource)); - checkCudaErrors(hipDestroyTextureObject(texTarget)); + TEX_DESTROY(texSource); + TEX_DESTROY(texTarget); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-hip/downscaleKernel.cuh b/src/opticalFlow-hip/downscaleKernel.cuh index ce728975d9..1dfc508b67 100644 --- a/src/opticalFlow-hip/downscaleKernel.cuh +++ b/src/opticalFlow-hip/downscaleKernel.cuh @@ -1,101 +1,42 @@ -/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * * Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * * Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * * Neither the name of NVIDIA CORPORATION nor the names of its - * contributors may be used to endorse or promote products derived - * from this software without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY - * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR - * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY - * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) */ - #include "common.h" -/////////////////////////////////////////////////////////////////////////////// -/// \brief downscale image -/// -/// CUDA kernel, relies heavily on texture unit -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// -__global__ void DownscaleKernel(int width, int height, int stride, float *out, - hipTextureObject_t texFine) { +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(tex, x, y) tex2D_sw(tex, x, y) +#define TEX_CREATE(texVar, src, width, height, stride) \ + Tex2DDesc texVar { src, width, height, stride } +#define TEX_DESTROY(texVar) /* no-op: no object to destroy */ + +__global__ void DownscaleKernel(int width, int height, int stride, + float *out, TEX_TYPE texFine) { const int ix = threadIdx.x + blockIdx.x * blockDim.x; const int iy = threadIdx.y + blockIdx.y * blockDim.y; - - if (ix >= width || iy >= height) { - return; - } + if (ix >= width || iy >= height) return; float dx = 1.0f / (float)width; float dy = 1.0f / (float)height; + float x = ((float)ix + 0.5f) * dx; + float y = ((float)iy + 0.5f) * dy; - float x = ((float)ix + 0.5f) * dx; - float y = ((float)iy + 0.5f) * dy; - - out[ix + iy * stride] = 0.25f * (tex2D(texFine, x - dx * 0.25f, y) + - tex2D(texFine, x + dx * 0.25f, y) + - tex2D(texFine, x, y - dy * 0.25f) + - tex2D(texFine, x, y + dy * 0.25f)); + out[ix + iy * stride] = 0.25f * (TEX_FETCH(texFine, x - dx*0.25f, y) + + TEX_FETCH(texFine, x + dx*0.25f, y) + + TEX_FETCH(texFine, x, y - dy*0.25f) + + TEX_FETCH(texFine, x, y + dy*0.25f)); } -/////////////////////////////////////////////////////////////////////////////// -/// \brief downscale image -/// -/// \param[in] src image to downscale -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// static void Downscale(const float *src, int width, int height, int stride, int newWidth, int newHeight, int newStride, float *out) { dim3 threads(32, 8); dim3 blocks(iDivUp(newWidth, threads.x), iDivUp(newHeight, threads.y)); - - hipTextureObject_t texFine; - hipResourceDesc texRes; - memset(&texRes, 0, sizeof(hipResourceDesc)); - - texRes.resType = hipResourceTypePitch2D; - texRes.res.pitch2D.devPtr = (void *)src; - texRes.res.pitch2D.desc = hipCreateChannelDesc(); - texRes.res.pitch2D.width = width; - texRes.res.pitch2D.height = height; - texRes.res.pitch2D.pitchInBytes = stride * sizeof(float); - - hipTextureDesc texDescr; - memset(&texDescr, 0, sizeof(hipTextureDesc)); - - texDescr.normalizedCoords = true; - texDescr.filterMode = hipFilterModeLinear; - texDescr.addressMode[0] = hipAddressModeMirror; - texDescr.addressMode[1] = hipAddressModeMirror; - texDescr.readMode = hipReadModeElementType; - - checkCudaErrors(hipCreateTextureObject(&texFine, &texRes, &texDescr, NULL)); - - DownscaleKernel<<>>(newWidth, newHeight, newStride, out, - texFine); - + TEX_CREATE(texFine, src, width, height, stride); + DownscaleKernel<<>>(newWidth, newHeight, newStride, out, texFine); checkCudaErrors(hipDeviceSynchronize()); - checkCudaErrors(hipDestroyTextureObject(texFine)); + TEX_DESTROY(texFine); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-hip/tex2d_compat.h b/src/opticalFlow-hip/tex2d_compat.h new file mode 100644 index 0000000000..ad3135ec11 --- /dev/null +++ b/src/opticalFlow-hip/tex2d_compat.h @@ -0,0 +1,42 @@ +// tex2d_compat.h — Portable tex2D for AMD GPUs + +#pragma once +#include + +// Descriptor passed to kernels in place of hipTextureObject_t +struct Tex2DDesc { + const float* ptr; // pitched device pointer + int width; // image width in elements + int height; // image height in elements + int pitch; // row stride in elements (pitchInBytes / sizeof(float)) +}; + +// Mirror repeat address clamping +__device__ inline int tex2d_mirror_clamp(int i, int size) { + int mi = ((i % (2 * size)) + 2 * size) % (2 * size); + return (mi >= size) ? (2 * size - 1 - mi) : mi; +} + +// Software bilinear tex2D — normalized coords [0,1), mirror addressing, linear filter +// Matches: normalizedCoords=true, filterMode=Linear, addressMode=Mirror +__device__ inline float tex2D_sw(const Tex2DDesc& tex, float u, float v) { + // Convert normalized coords to pixel space + float fx = u * (float)tex.width - 0.5f; + float fy = v * (float)tex.height - 0.5f; + + int x0 = (int)floorf(fx), y0 = (int)floorf(fy); + float tx = fx - (float)x0, ty = fy - (float)y0; + + int cx0 = tex2d_mirror_clamp(x0, tex.width); + int cx1 = tex2d_mirror_clamp(x0+1, tex.width); + int cy0 = tex2d_mirror_clamp(y0, tex.height); + int cy1 = tex2d_mirror_clamp(y0+1, tex.height); + + float v00 = tex.ptr[cy0 * tex.pitch + cx0]; + float v10 = tex.ptr[cy0 * tex.pitch + cx1]; + float v01 = tex.ptr[cy1 * tex.pitch + cx0]; + float v11 = tex.ptr[cy1 * tex.pitch + cx1]; + + return (1.f-tx)*(1.f-ty)*v00 + tx*(1.f-ty)*v10 + + (1.f-tx)*ty *v01 + tx*ty *v11; +} diff --git a/src/opticalFlow-hip/upscaleKernel.cuh b/src/opticalFlow-hip/upscaleKernel.cuh index 4ff2c5b915..c7d71cc51b 100644 --- a/src/opticalFlow-hip/upscaleKernel.cuh +++ b/src/opticalFlow-hip/upscaleKernel.cuh @@ -1,99 +1,34 @@ -/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * * Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * * Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * * Neither the name of NVIDIA CORPORATION nor the names of its - * contributors may be used to endorse or promote products derived - * from this software without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY - * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR - * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY - * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) */ - #include "common.h" - -/////////////////////////////////////////////////////////////////////////////// -/// \brief upscale one component of a displacement field, CUDA kernel -/// \param[in] width field width -/// \param[in] height field height -/// \param[in] stride field stride -/// \param[in] scale scale factor (multiplier) -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// -__global__ void UpscaleKernel(int width, int height, int stride, float scale, - float *out, hipTextureObject_t texCoarse) { +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(t,x,y) tex2D_sw(t,x,y) +#define TEX_CREATE(v,s,w,h,st) Tex2DDesc v{s,w,h,st} +#define TEX_DESTROY(v) + +__global__ void UpscaleKernel(int width, int height, int stride, + float scale, float *out, TEX_TYPE texCoarse) { const int ix = threadIdx.x + blockIdx.x * blockDim.x; const int iy = threadIdx.y + blockIdx.y * blockDim.y; - if (ix >= width || iy >= height) return; - float x = ((float)ix + 0.5f) / (float)width; float y = ((float)iy + 0.5f) / (float)height; - - // exploit hardware interpolation - // and scale interpolated vector to match next pyramid level resolution - out[ix + iy * stride] = tex2D(texCoarse, x, y) * scale; + out[ix + iy * stride] = TEX_FETCH(texCoarse, x, y) * scale; } -/////////////////////////////////////////////////////////////////////////////// -/// \brief upscale one component of a displacement field, kernel wrapper -/// \param[in] src field component to upscale -/// \param[in] width field current width -/// \param[in] height field current height -/// \param[in] stride field current stride -/// \param[in] newWidth field new width -/// \param[in] newHeight field new height -/// \param[in] newStride field new stride -/// \param[in] scale value scale factor (multiplier) -/// \param[out] out upscaled field component -/////////////////////////////////////////////////////////////////////////////// static void Upscale(const float *src, int width, int height, int stride, - int newWidth, int newHeight, int newStride, float scale, - float *out) { + int newWidth, int newHeight, int newStride, + float scale, float *out) { dim3 threads(32, 8); dim3 blocks(iDivUp(newWidth, threads.x), iDivUp(newHeight, threads.y)); - - hipTextureObject_t texCoarse; - hipResourceDesc texRes; - memset(&texRes, 0, sizeof(hipResourceDesc)); - - texRes.resType = hipResourceTypePitch2D; - texRes.res.pitch2D.devPtr = (void *)src; - texRes.res.pitch2D.desc = hipCreateChannelDesc(); - texRes.res.pitch2D.width = width; - texRes.res.pitch2D.height = height; - texRes.res.pitch2D.pitchInBytes = stride * sizeof(float); - - hipTextureDesc texDescr; - memset(&texDescr, 0, sizeof(hipTextureDesc)); - - texDescr.normalizedCoords = true; - texDescr.filterMode = hipFilterModeLinear; - texDescr.addressMode[0] = hipAddressModeMirror; - texDescr.addressMode[1] = hipAddressModeMirror; - texDescr.readMode = hipReadModeElementType; - - checkCudaErrors( - hipCreateTextureObject(&texCoarse, &texRes, &texDescr, NULL)); - - UpscaleKernel<<>>(newWidth, newHeight, newStride, scale, out, - texCoarse); - + TEX_CREATE(texCoarse, src, width, height, stride); + UpscaleKernel<<>>(newWidth, newHeight, newStride, scale, out, texCoarse); checkCudaErrors(hipDeviceSynchronize()); - checkCudaErrors(hipDestroyTextureObject(texCoarse)); + TEX_DESTROY(texCoarse); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-hip/warpingKernel.cuh b/src/opticalFlow-hip/warpingKernel.cuh index af97cf72db..32daef9c22 100644 --- a/src/opticalFlow-hip/warpingKernel.cuh +++ b/src/opticalFlow-hip/warpingKernel.cuh @@ -1,103 +1,35 @@ -/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * * Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * * Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * * Neither the name of NVIDIA CORPORATION nor the names of its - * contributors may be used to endorse or promote products derived - * from this software without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY - * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR - * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY - * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) */ - #include "common.h" - -/////////////////////////////////////////////////////////////////////////////// -/// \brief warp image with a given displacement field, CUDA kernel. -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[in] u horizontal displacement -/// \param[in] v vertical displacement -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// -__global__ void WarpingKernel(int width, int height, int stride, const float *u, - const float *v, float *out, - hipTextureObject_t texToWarp) { - const int ix = threadIdx.x + blockIdx.x * blockDim.x; - const int iy = threadIdx.y + blockIdx.y * blockDim.y; - +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(t,x,y) tex2D_sw(t,x,y) +#define TEX_CREATE(v,s,w,h,st) Tex2DDesc v{s,w,h,st} +#define TEX_DESTROY(v) + +__global__ void WarpingKernel(int width, int height, int stride, + const float *u, const float *v, + float *out, TEX_TYPE texToWarp) { + const int ix = threadIdx.x + blockIdx.x * blockDim.x; + const int iy = threadIdx.y + blockIdx.y * blockDim.y; const int pos = ix + iy * stride; - if (ix >= width || iy >= height) return; - float x = ((float)ix + u[pos] + 0.5f) / (float)width; float y = ((float)iy + v[pos] + 0.5f) / (float)height; - - out[pos] = tex2D(texToWarp, x, y); + out[pos] = TEX_FETCH(texToWarp, x, y); } -/////////////////////////////////////////////////////////////////////////////// -/// \brief warp image with provided vector field, CUDA kernel wrapper. -/// -/// For each output pixel there is a vector which tells which pixel -/// from a source image should be mapped to this particular output -/// pixel. -/// It is assumed that images and the vector field have the same stride and -/// resolution. -/// \param[in] src source image -/// \param[in] w width -/// \param[in] h height -/// \param[in] s stride -/// \param[in] u horizontal displacement -/// \param[in] v vertical displacement -/// \param[out] out warped image -/////////////////////////////////////////////////////////////////////////////// -static void WarpImage(const float *src, int w, int h, int s, const float *u, - const float *v, float *out) { +static void WarpImage(const float *src, int w, int h, int s, + const float *u, const float *v, float *out) { dim3 threads(32, 6); dim3 blocks(iDivUp(w, threads.x), iDivUp(h, threads.y)); - - hipTextureObject_t texToWarp; - hipResourceDesc texRes; - memset(&texRes, 0, sizeof(hipResourceDesc)); - - texRes.resType = hipResourceTypePitch2D; - texRes.res.pitch2D.devPtr = (void *)src; - texRes.res.pitch2D.desc = hipCreateChannelDesc(); - texRes.res.pitch2D.width = w; - texRes.res.pitch2D.height = h; - texRes.res.pitch2D.pitchInBytes = s * sizeof(float); - - hipTextureDesc texDescr; - memset(&texDescr, 0, sizeof(hipTextureDesc)); - - texDescr.normalizedCoords = true; - texDescr.filterMode = hipFilterModeLinear; - texDescr.addressMode[0] = hipAddressModeMirror; - texDescr.addressMode[1] = hipAddressModeMirror; - texDescr.readMode = hipReadModeElementType; - - checkCudaErrors( - hipCreateTextureObject(&texToWarp, &texRes, &texDescr, NULL)); - + TEX_CREATE(texToWarp, src, w, h, s); WarpingKernel<<>>(w, h, s, u, v, out, texToWarp); - checkCudaErrors(hipDeviceSynchronize()); - checkCudaErrors(hipDestroyTextureObject(texToWarp)); + TEX_DESTROY(texToWarp); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-omp/CMakeLists.txt b/src/opticalFlow-omp/CMakeLists.txt new file mode 100644 index 0000000000..d8c462f6e5 --- /dev/null +++ b/src/opticalFlow-omp/CMakeLists.txt @@ -0,0 +1,11 @@ +# opticalFlow-omp/CMakeLists.txt + +set(INC_DIR "${CMAKE_CURRENT_LIST_DIR}/../opticalFlow-cuda/Common") + +add_hecbench_benchmark( + NAME opticalFlow + MODEL omp + SOURCES flowGold.cpp flowOMP.cpp main.cpp + INCLUDE_DIRS ${INC_DIR} + CATEGORIES algorithms +) diff --git a/src/opticalFlow-omp/Makefile b/src/opticalFlow-omp/Makefile new file mode 100644 index 0000000000..9828b25693 --- /dev/null +++ b/src/opticalFlow-omp/Makefile @@ -0,0 +1,59 @@ +#=============================================================================== +# User Options +#=============================================================================== + +CC = icpx +OPTIMIZE = yes +DEBUG = no +DEVICE = gpu +LAUNCHER = + +#=============================================================================== +# Program name & source code list +#=============================================================================== + +program = main + +source = main.cpp flowOMP.cpp flowGold.cpp + +obj = $(source:.cpp=.o) + +#=============================================================================== +# Sets Flags +#=============================================================================== + +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../opticalFlow-cuda/Common + +LDFLAGS = + +ifeq ($(DEBUG),yes) + CFLAGS += -g + LDFLAGS += -g +endif + +ifeq ($(OPTIMIZE),yes) + CFLAGS += -O3 +endif + +ifeq ($(DEVICE),gpu) + # OpenMP 4.5 target offload (map/is_device_ptr; no omp_target_alloc) + CFLAGS += -fiopenmp -fopenmp-targets=spir64 -fopenmp-version=45 -D__STRICT_ANSI__ +else + CFLAGS += -qopenmp +endif + +#=============================================================================== +# Targets to Build +#=============================================================================== + +$(program): $(obj) + $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) + +%.o: %.cpp *.h *.hpp Makefile + $(CC) $(CFLAGS) -c $< -o $@ + +clean: + rm -rf $(program) $(obj) *.flo + +run: $(program) + $(LAUNCHER) ./$(program) ../opticalFlow-cuda/data/frame10.ppm ../opticalFlow-cuda/data/frame11.ppm diff --git a/src/opticalFlow-omp/Makefile.aomp b/src/opticalFlow-omp/Makefile.aomp new file mode 100644 index 0000000000..130564d49a --- /dev/null +++ b/src/opticalFlow-omp/Makefile.aomp @@ -0,0 +1,64 @@ +#=============================================================================== +# User Options +#=============================================================================== + +# Compiler can be set below, or via environment variable +CC = clang++ +OPTIMIZE = yes +DEBUG = no +DEVICE = gpu +ARCH = gfx906 +LAUNCHER = + +#=============================================================================== +# Program name & source code list +#=============================================================================== + +program = main + +source = main.cpp flowOMP.cpp flowGold.cpp + +obj = $(source:.cpp=.o) + +#=============================================================================== +# Sets Flags +#=============================================================================== + +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../opticalFlow-cuda/Common + +LDFLAGS = + +ifeq ($(DEBUG),yes) + CFLAGS += -g + LDFLAGS += -g +endif + +ifeq ($(OPTIMIZE),yes) + CFLAGS += -O3 +endif + +ifeq ($(DEVICE),gpu) + CFLAGS += -target x86_64-pc-linux-gnu \ + -fopenmp \ + -fopenmp-targets=amdgcn-amd-amdhsa \ + -Xopenmp-target=amdgcn-amd-amdhsa \ + -march=$(ARCH) +else + CFLAGS += -fopenmp +endif + +#=============================================================================== +# Targets to Build +#=============================================================================== + +$(program): $(obj) + $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) + +%.o: %.cpp *.h *.hpp Makefile.aomp + $(CC) $(CFLAGS) -c $< -o $@ + +clean: + rm -rf $(program) $(obj) *.flo + +run: $(program) + $(LAUNCHER) ./$(program) ../opticalFlow-cuda/data/frame10.ppm ../opticalFlow-cuda/data/frame11.ppm diff --git a/src/opticalFlow-omp/Makefile.nvc b/src/opticalFlow-omp/Makefile.nvc new file mode 100644 index 0000000000..55fede7766 --- /dev/null +++ b/src/opticalFlow-omp/Makefile.nvc @@ -0,0 +1,60 @@ +#=============================================================================== +# User Options +#=============================================================================== + +# Compiler can be set below, or via environment variable +CC = nvc++ +OPTIMIZE = yes +DEBUG = no +DEVICE = gpu +SM = cc70 +LAUNCHER = + +#=============================================================================== +# Program name & source code list +#=============================================================================== + +program = main + +source = main.cpp flowOMP.cpp flowGold.cpp + +obj = $(source:.cpp=.o) + +#=============================================================================== +# Sets Flags +#=============================================================================== + +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../opticalFlow-cuda/Common + +LDFLAGS = + +ifeq ($(DEBUG),yes) + CFLAGS += -g + LDFLAGS += -g +endif + +ifeq ($(OPTIMIZE),yes) + CFLAGS += -O3 +endif + +ifeq ($(DEVICE),gpu) + CFLAGS += -Minfo -mp=gpu -gpu=$(SM) +else + CFLAGS += -qopenmp +endif + +#=============================================================================== +# Targets to Build +#=============================================================================== + +$(program): $(obj) + $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) + +%.o: %.cpp *.h *.hpp Makefile.nvc + $(CC) $(CFLAGS) -c $< -o $@ + +clean: + rm -rf $(program) $(obj) *.flo + +run: $(program) + $(LAUNCHER) ./$(program) ../opticalFlow-cuda/data/frame10.ppm ../opticalFlow-cuda/data/frame11.ppm diff --git a/src/opticalFlow-omp/addKernel.hpp b/src/opticalFlow-omp/addKernel.hpp new file mode 100644 index 0000000000..7d51e1e25b --- /dev/null +++ b/src/opticalFlow-omp/addKernel.hpp @@ -0,0 +1,10 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. + */ +#include "common.h" + +static void Add(const float *op1, const float *op2, int count, float *sum) { + #pragma omp target teams distribute parallel for + for (int pos = 0; pos < count; ++pos) { + sum[pos] = op1[pos] + op2[pos]; + } +} diff --git a/src/opticalFlow-omp/common.h b/src/opticalFlow-omp/common.h new file mode 100644 index 0000000000..34ac9ac938 --- /dev/null +++ b/src/opticalFlow-omp/common.h @@ -0,0 +1,30 @@ +// Header for common includes and utility functions + +#ifndef COMMON_H +#define COMMON_H + +#include +#include +#include +#include +#include + +const int StrideAlignment = 256; + +inline int iAlignUp(int n, int m = StrideAlignment) { + int mod = n % m; + if (mod) + return n + m - mod; + return n; +} + +inline int iDivUp(int n, int m) { return (n + m - 1) / m; } + +template +inline void Swap(T &a, T &b) { + T t = a; + a = b; + b = t; +} + +#endif diff --git a/src/opticalFlow-omp/derivativesKernel.hpp b/src/opticalFlow-omp/derivativesKernel.hpp new file mode 100644 index 0000000000..387b6fc870 --- /dev/null +++ b/src/opticalFlow-omp/derivativesKernel.hpp @@ -0,0 +1,53 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) + */ +#include "common.h" +#include "tex2d_compat.h" + +static void ComputeDerivatives(const float *I0, const float *I1, + int w, int h, int s, + float *Ix, float *Iy, float *Iz) { + TEX_CREATE(texSource, w, h, s); + TEX_CREATE(texTarget, w, h, s); + + #pragma omp target teams distribute parallel for collapse(2) map(to: texSource, texTarget) + for (int iy = 0; iy < h; ++iy) { + for (int ix = 0; ix < w; ++ix) { + const int pos = ix + iy * s; + float dx = 1.0f / (float)w, dy = 1.0f / (float)h; + float x = ((float)ix + 0.5f) * dx; + float y = ((float)iy + 0.5f) * dy; + + float t0, t1; + t0 = TEX_FETCH(texSource, I0, x-2.f*dx, y); + t0 -= TEX_FETCH(texSource, I0, x-1.f*dx, y) * 8.f; + t0 += TEX_FETCH(texSource, I0, x+1.f*dx, y) * 8.f; + t0 -= TEX_FETCH(texSource, I0, x+2.f*dx, y); + t0 /= 12.f; + + t1 = TEX_FETCH(texTarget, I1, x-2.f*dx, y); + t1 -= TEX_FETCH(texTarget, I1, x-1.f*dx, y) * 8.f; + t1 += TEX_FETCH(texTarget, I1, x+1.f*dx, y) * 8.f; + t1 -= TEX_FETCH(texTarget, I1, x+2.f*dx, y); + t1 /= 12.f; + Ix[pos] = (t0 + t1) * 0.5f; + + Iz[pos] = TEX_FETCH(texTarget, I1, x, y) - TEX_FETCH(texSource, I0, x, y); + + t0 = TEX_FETCH(texSource, I0, x, y-2.f*dy); + t0 -= TEX_FETCH(texSource, I0, x, y-1.f*dy) * 8.f; + t0 += TEX_FETCH(texSource, I0, x, y+1.f*dy) * 8.f; + t0 -= TEX_FETCH(texSource, I0, x, y+2.f*dy); + t0 /= 12.f; + + t1 = TEX_FETCH(texTarget, I1, x, y-2.f*dy); + t1 -= TEX_FETCH(texTarget, I1, x, y-1.f*dy) * 8.f; + t1 += TEX_FETCH(texTarget, I1, x, y+1.f*dy) * 8.f; + t1 -= TEX_FETCH(texTarget, I1, x, y+2.f*dy); + t1 /= 12.f; + Iy[pos] = (t0 + t1) * 0.5f; + } + + } + TEX_DESTROY(texSource); + TEX_DESTROY(texTarget); +} diff --git a/src/opticalFlow-omp/downscaleKernel.hpp b/src/opticalFlow-omp/downscaleKernel.hpp new file mode 100644 index 0000000000..0524ae2091 --- /dev/null +++ b/src/opticalFlow-omp/downscaleKernel.hpp @@ -0,0 +1,26 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) + */ +#include "common.h" +#include "tex2d_compat.h" + +static void Downscale(const float *src, int width, int height, int stride, + int newWidth, int newHeight, int newStride, float *out) { + + TEX_CREATE(texFine, width, height, stride); + + #pragma omp target teams distribute parallel for collapse(2) map(to:texFine) + for (int iy = 0; iy < newHeight; ++iy) { + for (int ix = 0; ix < newWidth; ++ix) { + float dx = 1.0f / (float)newWidth; + float dy = 1.0f / (float)newHeight; + float x = ((float)ix + 0.5f) * dx; + float y = ((float)iy + 0.5f) * dy; + + out[ix + iy * newStride] = 0.25f * (TEX_FETCH(texFine, src, x - dx*0.25f, y) + + TEX_FETCH(texFine, src, x + dx*0.25f, y) + + TEX_FETCH(texFine, src, x, y - dy*0.25f) + + TEX_FETCH(texFine, src, x, y + dy*0.25f)); + } + } + TEX_DESTROY(texFine); +} diff --git a/src/opticalFlow-omp/flowGold.cpp b/src/opticalFlow-omp/flowGold.cpp new file mode 100644 index 0000000000..98a72129e5 --- /dev/null +++ b/src/opticalFlow-omp/flowGold.cpp @@ -0,0 +1,287 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. + */ +#include "common.h" +#include "flowGold.h" + +inline float Tex2D(const float *t, int w, int h, int s, float x, float y) { + float intPartX, intPartY; + float dx = fabsf(modff(x, &intPartX)); + float dy = fabsf(modff(y, &intPartY)); + + int ix0 = (int)intPartX; + int iy0 = (int)intPartY; + + if (ix0 < 0) ix0 = abs(ix0 + 1); + if (iy0 < 0) iy0 = abs(iy0 + 1); + if (ix0 >= w) ix0 = w * 2 - ix0 - 1; + if (iy0 >= h) iy0 = h * 2 - iy0 - 1; + + int ix1 = ix0 + 1; + int iy1 = iy0 + 1; + + if (ix1 >= w) ix1 = w * 2 - ix1 - 1; + if (iy1 >= h) iy1 = h * 2 - iy1 - 1; + + float res = t[ix0 + iy0 * s] * (1.0f - dx) * (1.0f - dy); + res += t[ix1 + iy0 * s] * dx * (1.0f - dy); + res += t[ix0 + iy1 * s] * (1.0f - dx) * dy; + res += t[ix1 + iy1 * s] * dx * dy; + + return res; +} + +inline float Tex2Di(const float *src, int w, int h, int s, int x, int y) { + if (x < 0) x = abs(x + 1); + if (y < 0) y = abs(y + 1); + if (x >= w) x = w * 2 - x - 1; + if (y >= h) y = h * 2 - y - 1; + + return src[x + y * s]; +} + +static void Downscale(const float *src, int width, int height, int stride, + int newWidth, int newHeight, int newStride, float *out) { + for (int i = 0; i < newHeight; ++i) { + for (int j = 0; j < newWidth; ++j) { + const int srcX = j * 2; + const int srcY = i * 2; + float sum; + sum = Tex2Di(src, width, height, stride, srcX + 0, srcY + 0); + sum += Tex2Di(src, width, height, stride, srcX + 0, srcY + 1); + sum += Tex2Di(src, width, height, stride, srcX + 1, srcY + 0); + sum += Tex2Di(src, width, height, stride, srcX + 1, srcY + 1); + sum *= 0.25f; + out[j + i * newStride] = sum; + } + } +} + +static void Upscale(const float *src, int width, int height, int stride, + int newWidth, int newHeight, int newStride, float scale, + float *out) { + for (int i = 0; i < newHeight; ++i) { + for (int j = 0; j < newWidth; ++j) { + float x = ((float)j - 0.5f) * 0.5f; + float y = ((float)i - 0.5f) * 0.5f; + out[j + i * newStride] = Tex2D(src, width, height, stride, x, y) * scale; + } + } +} + +static void WarpImage(const float *src, int w, int h, int s, const float *u, + const float *v, float *out) { + for (int i = 0; i < h; ++i) { + for (int j = 0; j < w; ++j) { + const int pos = j + i * s; + float x = (float)j + u[pos]; + float y = (float)i + v[pos]; + out[pos] = Tex2D(src, w, h, s, x, y); + } + } +} + +static void ComputeDerivatives(const float *I0, const float *I1, int w, int h, + int s, float *Ix, float *Iy, float *Iz) { + for (int i = 0; i < h; ++i) { + for (int j = 0; j < w; ++j) { + const int pos = j + i * s; + float t0, t1; + + t0 = Tex2Di(I0, w, h, s, j - 2, i); + t0 -= Tex2Di(I0, w, h, s, j - 1, i) * 8.0f; + t0 += Tex2Di(I0, w, h, s, j + 1, i) * 8.0f; + t0 -= Tex2Di(I0, w, h, s, j + 2, i); + t0 /= 12.0f; + + t1 = Tex2Di(I1, w, h, s, j - 2, i); + t1 -= Tex2Di(I1, w, h, s, j - 1, i) * 8.0f; + t1 += Tex2Di(I1, w, h, s, j + 1, i) * 8.0f; + t1 -= Tex2Di(I1, w, h, s, j + 2, i); + t1 /= 12.0f; + + Ix[pos] = (t0 + t1) * 0.5f; + Iz[pos] = I1[pos] - I0[pos]; + + t0 = Tex2Di(I0, w, h, s, j, i - 2); + t0 -= Tex2Di(I0, w, h, s, j, i - 1) * 8.0f; + t0 += Tex2Di(I0, w, h, s, j, i + 1) * 8.0f; + t0 -= Tex2Di(I0, w, h, s, j, i + 2); + t0 /= 12.0f; + + t1 = Tex2Di(I1, w, h, s, j, i - 2); + t1 -= Tex2Di(I1, w, h, s, j, i - 1) * 8.0f; + t1 += Tex2Di(I1, w, h, s, j, i + 1) * 8.0f; + t1 -= Tex2Di(I1, w, h, s, j, i + 2); + t1 /= 12.0f; + + Iy[pos] = (t0 + t1) * 0.5f; + } + } +} + +static void SolveForUpdate(const float *du0, const float *dv0, const float *Ix, + const float *Iy, const float *Iz, int w, int h, + int s, float alpha, float *du1, float *dv1) { + for (int i = 0; i < h; ++i) { + for (int j = 0; j < w; ++j) { + const int pos = j + i * s; + int left, right, up, down; + + if (j != 0) + left = pos - 1; + else + left = pos; + + if (j != w - 1) + right = pos + 1; + else + right = pos; + + if (i != 0) + down = pos - s; + else + down = pos; + + if (i != h - 1) + up = pos + s; + else + up = pos; + + float sumU = (du0[left] + du0[right] + du0[up] + du0[down]) * 0.25f; + float sumV = (dv0[left] + dv0[right] + dv0[up] + dv0[down]) * 0.25f; + + float frac = (Ix[pos] * sumU + Iy[pos] * sumV + Iz[pos]) / + (Ix[pos] * Ix[pos] + Iy[pos] * Iy[pos] + alpha); + + du1[pos] = sumU - Ix[pos] * frac; + dv1[pos] = sumV - Iy[pos] * frac; + } + } +} + +void ComputeFlowGold(const float *I0, const float *I1, int width, int height, + int stride, float alpha, int nLevels, int nWarpIters, + int nSolverIters, float *u, float *v) { + printf("Computing optical flow on CPU...\n"); + + float *u0 = u; + float *v0 = v; + + const float **pI0 = new const float *[nLevels]; + const float **pI1 = new const float *[nLevels]; + + int *pW = new int[nLevels]; + int *pH = new int[nLevels]; + int *pS = new int[nLevels]; + + const int pixelCountAligned = height * stride; + + float *tmp = new float[pixelCountAligned]; + float *du0 = new float[pixelCountAligned]; + float *dv0 = new float[pixelCountAligned]; + float *du1 = new float[pixelCountAligned]; + float *dv1 = new float[pixelCountAligned]; + float *Ix = new float[pixelCountAligned]; + float *Iy = new float[pixelCountAligned]; + float *Iz = new float[pixelCountAligned]; + float *nu = new float[pixelCountAligned]; + float *nv = new float[pixelCountAligned]; + + int currentLevel = nLevels - 1; + pI0[currentLevel] = I0; + pI1[currentLevel] = I1; + + pW[currentLevel] = width; + pH[currentLevel] = height; + pS[currentLevel] = stride; + + for (; currentLevel > 0; --currentLevel) { + int nw = pW[currentLevel] / 2; + int nh = pH[currentLevel] / 2; + int ns = iAlignUp(nw); + pI0[currentLevel - 1] = new float[ns * nh]; + pI1[currentLevel - 1] = new float[ns * nh]; + + Downscale(pI0[currentLevel], pW[currentLevel], pH[currentLevel], + pS[currentLevel], nw, nh, ns, (float *)pI0[currentLevel - 1]); + Downscale(pI1[currentLevel], pW[currentLevel], pH[currentLevel], + pS[currentLevel], nw, nh, ns, (float *)pI1[currentLevel - 1]); + + pW[currentLevel - 1] = nw; + pH[currentLevel - 1] = nh; + pS[currentLevel - 1] = ns; + } + + memset(u, 0, stride * height * sizeof(float)); + memset(v, 0, stride * height * sizeof(float)); + + for (; currentLevel < nLevels; ++currentLevel) { + for (int warpIter = 0; warpIter < nWarpIters; ++warpIter) { + memset(du0, 0, pixelCountAligned * sizeof(float)); + memset(dv0, 0, pixelCountAligned * sizeof(float)); + memset(du1, 0, pixelCountAligned * sizeof(float)); + memset(dv1, 0, pixelCountAligned * sizeof(float)); + + WarpImage(pI1[currentLevel], pW[currentLevel], pH[currentLevel], + pS[currentLevel], u, v, tmp); + + ComputeDerivatives(pI0[currentLevel], tmp, pW[currentLevel], + pH[currentLevel], pS[currentLevel], Ix, Iy, Iz); + + for (int iter = 0; iter < nSolverIters; ++iter) { + SolveForUpdate(du0, dv0, Ix, Iy, Iz, pW[currentLevel], pH[currentLevel], + pS[currentLevel], alpha, du1, dv1); + Swap(du0, du1); + Swap(dv0, dv1); + } + + for (int i = 0; i < pH[currentLevel] * pS[currentLevel]; ++i) { + u[i] += du0[i]; + v[i] += dv0[i]; + } + } + + if (currentLevel != nLevels - 1) { + float scaleX = (float)pW[currentLevel + 1] / (float)pW[currentLevel]; + Upscale(u, pW[currentLevel], pH[currentLevel], pS[currentLevel], + pW[currentLevel + 1], pH[currentLevel + 1], pS[currentLevel + 1], + scaleX, nu); + + float scaleY = (float)pH[currentLevel + 1] / (float)pH[currentLevel]; + Upscale(v, pW[currentLevel], pH[currentLevel], pS[currentLevel], + pW[currentLevel + 1], pH[currentLevel + 1], pS[currentLevel + 1], + scaleY, nv); + + Swap(u, nu); + Swap(v, nv); + } + } + + if (u != u0) { + memcpy(u0, u, pixelCountAligned * sizeof(float)); + memcpy(v0, v, pixelCountAligned * sizeof(float)); + Swap(u, nu); + Swap(v, nv); + } + + for (int i = 0; i < nLevels - 1; ++i) { + delete[] pI0[i]; + delete[] pI1[i]; + } + + delete[] pI0; + delete[] pI1; + delete[] pW; + delete[] pH; + delete[] pS; + delete[] tmp; + delete[] du0; + delete[] dv0; + delete[] du1; + delete[] dv1; + delete[] Ix; + delete[] Iy; + delete[] Iz; + delete[] nu; + delete[] nv; +} diff --git a/src/opticalFlow-omp/flowGold.h b/src/opticalFlow-omp/flowGold.h new file mode 100644 index 0000000000..371c784bbb --- /dev/null +++ b/src/opticalFlow-omp/flowGold.h @@ -0,0 +1,8 @@ +#ifndef FLOW_GOLD_H +#define FLOW_GOLD_H + +void ComputeFlowGold(const float *I0, const float *I1, int width, int height, + int stride, float alpha, int nLevels, int nWarpIters, + int nIters, float *u, float *v); + +#endif diff --git a/src/opticalFlow-omp/flowOMP.cpp b/src/opticalFlow-omp/flowOMP.cpp new file mode 100644 index 0000000000..28ba1762e3 --- /dev/null +++ b/src/opticalFlow-omp/flowOMP.cpp @@ -0,0 +1,177 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. + */ +#include "common.h" + +#include "addKernel.hpp" +#include "derivativesKernel.hpp" +#include "downscaleKernel.hpp" +#include "solverKernel.hpp" +#include "upscaleKernel.hpp" +#include "warpingKernel.hpp" + +void ComputeFlowOMP(const float *I0, const float *I1, int width, int height, + int stride, float alpha, int nLevels, int nWarpIters, + int nSolverIters, float *d_u, float *d_v) { + printf("Computing optical flow on GPU (OpenMP target)...\n"); + + float **pI0 = new float *[nLevels]; + float **pI1 = new float *[nLevels]; + + int *pW = new int[nLevels]; + int *pH = new int[nLevels]; + int *pS = new int[nLevels]; + + const int pixelCount = stride * height; + + float *d_tmp = new float[pixelCount]; + float *d_du0 = new float[pixelCount]; + float *d_dv0 = new float[pixelCount]; + float *d_du1 = new float[pixelCount]; + float *d_dv1 = new float[pixelCount]; + float *d_Ix = new float[pixelCount]; + float *d_Iy = new float[pixelCount]; + float *d_Iz = new float[pixelCount]; + //float *d_u = new float[pixelCount]; + //float *d_v = new float[pixelCount]; + float *d_nu = new float[pixelCount]; + float *d_nv = new float[pixelCount]; + + #pragma omp target enter data map(alloc : d_tmp[0:pixelCount]) + #pragma omp target enter data map(alloc : d_du0[0:pixelCount]) + #pragma omp target enter data map(alloc : d_dv0[0:pixelCount]) + #pragma omp target enter data map(alloc : d_du1[0:pixelCount]) + #pragma omp target enter data map(alloc : d_dv1[0:pixelCount]) + #pragma omp target enter data map(alloc : d_Ix[0:pixelCount]) + #pragma omp target enter data map(alloc : d_Iy[0:pixelCount]) + #pragma omp target enter data map(alloc : d_Iz[0:pixelCount]) + #pragma omp target enter data map(alloc : d_u[0:pixelCount]) + #pragma omp target enter data map(alloc : d_v[0:pixelCount]) + #pragma omp target enter data map(alloc : d_nu[0:pixelCount]) + #pragma omp target enter data map(alloc : d_nv[0:pixelCount]) + + int currentLevel = nLevels - 1; + + pI0[currentLevel] = new float [pixelCount]; + pI1[currentLevel] = new float [pixelCount]; + + memcpy(pI0[currentLevel], I0, pixelCount * sizeof(float)); + memcpy(pI1[currentLevel], I1, pixelCount * sizeof(float)); + + #pragma omp target enter data map(to : pI0[currentLevel][0:pixelCount]) + #pragma omp target enter data map(to : pI1[currentLevel][0:pixelCount]) + + pW[currentLevel] = width; + pH[currentLevel] = height; + pS[currentLevel] = stride; + + for (; currentLevel > 0; --currentLevel) { + int nw = pW[currentLevel] / 2; + int nh = pH[currentLevel] / 2; + int ns = iAlignUp(nw); + + pI0[currentLevel-1] = new float [ns * nh]; + pI1[currentLevel-1] = new float [ns * nh]; + #pragma omp target enter data map(alloc : pI0[currentLevel-1][0:ns*nh]) + #pragma omp target enter data map(alloc : pI1[currentLevel-1][0:ns*nh]) + + Downscale(pI0[currentLevel], pW[currentLevel], pH[currentLevel], + pS[currentLevel], nw, nh, ns, pI0[currentLevel-1]); + Downscale(pI1[currentLevel], pW[currentLevel], pH[currentLevel], + pS[currentLevel], nw, nh, ns, pI1[currentLevel-1]); + pW[currentLevel - 1] = nw; + pH[currentLevel - 1] = nh; + pS[currentLevel - 1] = ns; + } + + #pragma omp target teams distribute parallel for + for (int i = 0; i < pixelCount; i++) { + d_u[i] = 0; + d_v[i] = 0; + } + + for (; currentLevel < nLevels; ++currentLevel) { + for (int warpIter = 0; warpIter < nWarpIters; ++warpIter) { + #pragma omp target teams distribute parallel for + for (int i = 0; i < pixelCount; i++) { + d_du0[i] = 0; + d_dv0[i] = 0; + d_du1[i] = 0; + d_dv1[i] = 0; + } + + WarpImage(pI1[currentLevel], pW[currentLevel], pH[currentLevel], pS[currentLevel], d_u, d_v, d_tmp); + + ComputeDerivatives(pI0[currentLevel], d_tmp, pW[currentLevel], pH[currentLevel], + pS[currentLevel], d_Ix, d_Iy, d_Iz); + + for (int iter = 0; iter < nSolverIters; ++iter) { + SolveForUpdate(d_du0, d_dv0, d_Ix, d_Iy, d_Iz, pW[currentLevel], + pH[currentLevel], pS[currentLevel], alpha, d_du1, d_dv1); + Swap(d_du0, d_du1); + Swap(d_dv0, d_dv1); + } + + Add(d_u, d_du0, pH[currentLevel] * pS[currentLevel], d_u); + Add(d_v, d_dv0, pH[currentLevel] * pS[currentLevel], d_v); + } + + if (currentLevel != nLevels - 1) { + float scaleX = (float)pW[currentLevel + 1] / (float)pW[currentLevel]; + + Upscale(d_u, pW[currentLevel], pH[currentLevel], pS[currentLevel], + pW[currentLevel + 1], pH[currentLevel + 1], pS[currentLevel + 1], + scaleX, d_nu); + + float scaleY = (float)pH[currentLevel + 1] / (float)pH[currentLevel]; + + Upscale(d_v, pW[currentLevel], pH[currentLevel], pS[currentLevel], + pW[currentLevel + 1], pH[currentLevel + 1], pS[currentLevel + 1], + scaleY, d_nv); + + Swap(d_u, d_nu); + Swap(d_v, d_nv); + } + } + + #pragma omp target update from (d_u[0:pixelCount]) + #pragma omp target update from (d_v[0:pixelCount]) + + for (int l = 0; l < nLevels; ++l) { + #pragma omp target exit data map(delete : pI0[l][0:pS[l]*pH[l]]) + #pragma omp target exit data map(delete : pI1[l][0:pS[l]*pH[l]]) + delete(pI0[l]); + delete(pI1[l]); + } + + #pragma omp target exit data map(delete : d_tmp[0:pixelCount]) + #pragma omp target exit data map(delete : d_du0[0:pixelCount]) + #pragma omp target exit data map(delete : d_dv0[0:pixelCount]) + #pragma omp target exit data map(delete : d_du1[0:pixelCount]) + #pragma omp target exit data map(delete : d_dv1[0:pixelCount]) + #pragma omp target exit data map(delete : d_Ix[0:pixelCount]) + #pragma omp target exit data map(delete : d_Iy[0:pixelCount]) + #pragma omp target exit data map(delete : d_Iz[0:pixelCount]) + #pragma omp target exit data map(delete : d_u[0:pixelCount]) + #pragma omp target exit data map(delete : d_v[0:pixelCount]) + #pragma omp target exit data map(delete : d_nu[0:pixelCount]) + #pragma omp target exit data map(delete : d_nv[0:pixelCount]) + + delete[] pI0; + delete[] pI1; + delete[] pW; + delete[] pH; + delete[] pS; + + delete[] d_tmp; + delete[] d_du0; + delete[] d_dv0; + delete[] d_du1; + delete[] d_dv1; + delete[] d_Ix; + delete[] d_Iy; + delete[] d_Iz; + //delete[] d_u; + //delete[] d_v; + delete[] d_nu; + delete[] d_nv; +} diff --git a/src/opticalFlow-omp/flowOMP.h b/src/opticalFlow-omp/flowOMP.h new file mode 100644 index 0000000000..562d6d9afb --- /dev/null +++ b/src/opticalFlow-omp/flowOMP.h @@ -0,0 +1,8 @@ +#ifndef FLOW_OMP_H +#define FLOW_OMP_H + +void ComputeFlowOMP(const float *I0, const float *I1, int width, int height, + int stride, float alpha, int nLevels, int nWarpIters, + int nSolverIters, float *u, float *v); + +#endif diff --git a/src/opticalFlow-omp/main.cpp b/src/opticalFlow-omp/main.cpp new file mode 100644 index 0000000000..d94e6885e6 --- /dev/null +++ b/src/opticalFlow-omp/main.cpp @@ -0,0 +1,158 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. + */ +const static char *const sSDKsample = "HSOpticalFlow"; + +const float THRESHOLD = 0.05f; + +#include "common.h" +#include "flowGold.h" +#include "flowOMP.h" +#include "helper_functions.h" +#include +#include + +using Time = std::chrono::steady_clock; +using float_ms = std::chrono::duration; + +void WriteFloFile(const char *name, int w, int h, int s, const float *u, + const float *v) { + FILE *stream = fopen(name, "wb"); + + if (stream == 0) { + printf("Could not save flow to \"%s\"\n", name); + return; + } + + float data = 202021.25f; + fwrite(&data, sizeof(float), 1, stream); + fwrite(&w, sizeof(w), 1, stream); + fwrite(&h, sizeof(h), 1, stream); + + for (int i = 0; i < h; ++i) { + for (int j = 0; j < w; ++j) { + const int pos = j + i * s; + fwrite(u + pos, sizeof(float), 1, stream); + fwrite(v + pos, sizeof(float), 1, stream); + } + } + + fclose(stream); +} + +bool LoadImageAsFP32(float *&img_data, int &img_w, int &img_h, int &img_s, + const char *name, const char *exePath) { + printf("Loading \"%s\" ...\n", name); + char *name_ = sdkFindFilePath(name, exePath); + + if (!name_) { + printf("File not found\n"); + return false; + } + + unsigned char *data = 0; + unsigned int w = 0, h = 0; + bool result = sdkLoadPPM4ub(name_, &data, &w, &h); + + if (result == false) { + printf("Invalid file format\n"); + return false; + } + + img_w = w; + img_h = h; + img_s = iAlignUp(img_w); + + img_data = new float[img_s * h]; + + const int widthStep = 4 * img_w; + + for (int i = 0; i < img_h; ++i) { + for (int j = 0; j < img_w; ++j) { + img_data[j + i * img_s] = ((float)data[j * 4 + i * widthStep]) / 255.0f; + } + } + + return true; +} + +bool CompareWithGold(int width, int height, int stride, const float *h_uGold, + const float *h_vGold, const float *h_u, const float *h_v) { + float error = 0.0f; + + for (int i = 0; i < height; ++i) { + for (int j = 0; j < width; ++j) { + const int pos = j + i * stride; + error += fabsf(h_u[pos] - h_uGold[pos]) + fabsf(h_v[pos] - h_vGold[pos]); + } + } + + error /= (float)(width * height); + + printf("L1 error : %.6f\n", error); + + return (error < THRESHOLD); +} + +int main(int argc, char **argv) { + printf("%s Starting...\n\n", sSDKsample); + + const char *const sourceFrameName = argv[1]; + const char *const targetFrameName = argv[2]; + + int width; + int height; + int stride; + + float *h_source; + float *h_target; + + if (!LoadImageAsFP32(h_source, width, height, stride, sourceFrameName, + argv[0])) { + exit(EXIT_FAILURE); + } + + if (!LoadImageAsFP32(h_target, width, height, stride, targetFrameName, + argv[0])) { + exit(EXIT_FAILURE); + } + + float *h_uGold = new float[stride * height]; + float *h_vGold = new float[stride * height]; + float *h_u = new float[stride * height]; + float *h_v = new float[stride * height]; + + const float alpha = 0.2f; + const int nLevels = 5; + const int nSolverIters = 500; + const int nWarpIters = 3; + + auto start = Time::now(); + ComputeFlowGold(h_source, h_target, width, height, stride, alpha, nLevels, + nWarpIters, nSolverIters, h_uGold, h_vGold); + auto stop = Time::now(); + + auto duration = std::chrono::duration_cast(stop - start).count(); + printf("Processing time on CPU: %f (ms)\n", duration); + + start = Time::now(); + ComputeFlowOMP(h_source, h_target, width, height, stride, alpha, nLevels, + nWarpIters, nSolverIters, h_u, h_v); + stop = Time::now(); + duration = std::chrono::duration_cast(stop - start).count(); + printf("Processing time on Device: %f (ms)\n", duration); + + bool status = + CompareWithGold(width, height, stride, h_uGold, h_vGold, h_u, h_v); + + WriteFloFile("FlowGPU.flo", width, height, stride, h_u, h_v); + WriteFloFile("FlowCPU.flo", width, height, stride, h_uGold, h_vGold); + + delete[] h_uGold; + delete[] h_vGold; + delete[] h_u; + delete[] h_v; + delete[] h_source; + delete[] h_target; + + exit(status ? EXIT_SUCCESS : EXIT_FAILURE); +} diff --git a/src/opticalFlow-omp/solverKernel.hpp b/src/opticalFlow-omp/solverKernel.hpp new file mode 100644 index 0000000000..45e79b1af3 --- /dev/null +++ b/src/opticalFlow-omp/solverKernel.hpp @@ -0,0 +1,38 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. + */ +#include "common.h" + +static void SolveForUpdate(const float *du0, const float *dv0, const float *Ix, + const float *Iy, const float *Iz, int w, + int h, int s, float alpha, float *du1, float *dv1) { + #pragma omp target teams distribute parallel for collapse(2) + for (int iy = 0; iy < h; ++iy) { + for (int ix = 0; ix < w; ++ix) { + const int pos = iy * s + ix; + + // Clamp-to-edge neighbour indices + const int x_left = std::max(ix - 1, 0); + const int x_right = std::min(ix + 1, w - 1); + const int y_up = std::min(iy + 1, h - 1); + const int y_down = std::max(iy - 1, 0); + + // Neighbour linear addresses + const int pos_left = iy * s + x_left; + const int pos_right = iy * s + x_right; + const int pos_up = y_up * s + ix; + const int pos_down = y_down * s + ix; + + // Jacobi stencil + const float sumU = (du0[pos_left] + du0[pos_right] + + du0[pos_up] + du0[pos_down]) * 0.25f; + const float sumV = (dv0[pos_left] + dv0[pos_right] + + dv0[pos_up] + dv0[pos_down]) * 0.25f; + + const float frac = (Ix[pos] * sumU + Iy[pos] * sumV + Iz[pos]) / + (Ix[pos] * Ix[pos] + Iy[pos] * Iy[pos] + alpha); + + du1[pos] = sumU - Ix[pos] * frac; + dv1[pos] = sumV - Iy[pos] * frac; + } + } +} diff --git a/src/opticalFlow-omp/tex2d_compat.h b/src/opticalFlow-omp/tex2d_compat.h new file mode 100644 index 0000000000..3b2d19e971 --- /dev/null +++ b/src/opticalFlow-omp/tex2d_compat.h @@ -0,0 +1,43 @@ +// tex2d_compat.h — Portable tex2D for OpenMP target + +#pragma once + +struct Tex2DDesc { + int width; + int height; + int pitch; +}; + +#pragma omp declare target +inline int tex2d_mirror_clamp(int i, int size) { + int mi = ((i % (2 * size)) + 2 * size) % (2 * size); + return (mi >= size) ? (2 * size - 1 - mi) : mi; +} + +inline float tex2D_sw(const Tex2DDesc &tex, const float *ptr, float u, float v) { + float fx = u * (float)tex.width - 0.5f; + float fy = v * (float)tex.height - 0.5f; + + int x0 = (int)floorf(fx), y0 = (int)floorf(fy); + float tx = fx - (float)x0, ty = fy - (float)y0; + + int cx0 = tex2d_mirror_clamp(x0, tex.width); + int cx1 = tex2d_mirror_clamp(x0 + 1, tex.width); + int cy0 = tex2d_mirror_clamp(y0, tex.height); + int cy1 = tex2d_mirror_clamp(y0 + 1, tex.height); + + float v00 = ptr[cy0 * tex.pitch + cx0]; + float v10 = ptr[cy0 * tex.pitch + cx1]; + float v01 = ptr[cy1 * tex.pitch + cx0]; + float v11 = ptr[cy1 * tex.pitch + cx1]; + + return (1.f - tx) * (1.f - ty) * v00 + tx * (1.f - ty) * v10 + + (1.f - tx) * ty * v01 + tx * ty * v11; +} + +#pragma omp end declare target + +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(t,p,x,y) tex2D_sw(t,p,x,y) +#define TEX_CREATE(v,w,h,s) Tex2DDesc v{w,h,s} +#define TEX_DESTROY(v) diff --git a/src/opticalFlow-omp/upscaleKernel.hpp b/src/opticalFlow-omp/upscaleKernel.hpp new file mode 100644 index 0000000000..6449b97029 --- /dev/null +++ b/src/opticalFlow-omp/upscaleKernel.hpp @@ -0,0 +1,20 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) + */ +#include "common.h" +#include "tex2d_compat.h" + +static void Upscale(const float *src, int width, int height, int stride, + int newWidth, int newHeight, int newStride, + float scale, float *out) { + TEX_CREATE(texCoarse, width, height, stride); + + #pragma omp target teams distribute parallel for collapse(2) map(to: texCoarse) + for (int iy = 0; iy < newHeight; ++iy) { + for (int ix = 0; ix < newWidth; ++ix) { + float x = ((float)ix + 0.5f) / (float)newWidth; + float y = ((float)iy + 0.5f) / (float)newHeight; + out[ix + iy * newStride] = TEX_FETCH(texCoarse, src, x, y) * scale; + } + } + TEX_DESTROY(texCoarse); +} diff --git a/src/opticalFlow-omp/warpingKernel.hpp b/src/opticalFlow-omp/warpingKernel.hpp new file mode 100644 index 0000000000..13a8863b29 --- /dev/null +++ b/src/opticalFlow-omp/warpingKernel.hpp @@ -0,0 +1,20 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) + */ +#include "common.h" +#include "tex2d_compat.h" + +static void WarpImage(const float *src, int w, int h, int s, + const float *u, const float *v, float *out) { + + TEX_CREATE(texToWarp, w, h, s); + #pragma omp target teams distribute parallel for collapse(2) map(to: texToWarp) + for (int iy = 0; iy < h; ++iy) { + for (int ix = 0; ix < w; ++ix) { + const int pos = ix + iy * s; + float x = ((float)ix + u[pos] + 0.5f) / (float)w; + float y = ((float)iy + v[pos] + 0.5f) / (float)h; + out[pos] = TEX_FETCH(texToWarp, src, x, y); + } + } + TEX_DESTROY(texToWarp); +} diff --git a/src/opticalFlow-sycl/COPYRIGHT b/src/opticalFlow-sycl/COPYRIGHT deleted file mode 100644 index 39418d7917..0000000000 --- a/src/opticalFlow-sycl/COPYRIGHT +++ /dev/null @@ -1,35 +0,0 @@ -//========================================================= -// Modifications Copyright © 2022 Intel Corporation -// -// SPDX-License-Identifier: BSD-3-Clause -//========================================================= - - -/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * * Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * * Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * * Neither the name of NVIDIA CORPORATION nor the names of its - * contributors may be used to endorse or promote products derived - * from this software without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY - * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR - * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR - * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY - * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ - - diff --git a/src/opticalFlow-sycl/Makefile b/src/opticalFlow-sycl/Makefile index 40a865cb8f..0d30ef74c1 100644 --- a/src/opticalFlow-sycl/Makefile +++ b/src/opticalFlow-sycl/Makefile @@ -72,7 +72,7 @@ endif $(program): $(obj) $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) -%.o: %.cpp *.h +%.o: %.cpp *.h *.hpp $(CC) $(CFLAGS) -c $< -o $@ clean: diff --git a/src/opticalFlow-sycl/addKernel.hpp b/src/opticalFlow-sycl/addKernel.hpp index 0603026b85..600e52db4c 100644 --- a/src/opticalFlow-sycl/addKernel.hpp +++ b/src/opticalFlow-sycl/addKernel.hpp @@ -1,3 +1,30 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * * Neither the name of NVIDIA CORPORATION nor the names of its + * contributors may be used to endorse or promote products derived + * from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY + * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY + * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + #include "common.h" /////////////////////////////////////////////////////////////////////////////// @@ -11,8 +38,7 @@ /////////////////////////////////////////////////////////////////////////////// void AddKernel(const float *op1, const float *op2, int count, float *sum, const sycl::nd_item<3> &item) { - const int pos = item.get_local_id(2) + - item.get_group(2) * item.get_local_range(2); + const int pos = item.get_global_id(2); if (pos >= count) return; diff --git a/src/opticalFlow-sycl/common.h b/src/opticalFlow-sycl/common.h index 935610e3fe..03d3115c67 100644 --- a/src/opticalFlow-sycl/common.h +++ b/src/opticalFlow-sycl/common.h @@ -18,7 +18,7 @@ /////////////////////////////////////////////////////////////////////////////// // Common constants /////////////////////////////////////////////////////////////////////////////// -const int StrideAlignment = 8; +const int StrideAlignment = 256; /////////////////////////////////////////////////////////////////////////////// // Common functions diff --git a/src/opticalFlow-sycl/derivativesKernel.hpp b/src/opticalFlow-sycl/derivativesKernel.hpp index d203af116b..a16ad7dfca 100644 --- a/src/opticalFlow-sycl/derivativesKernel.hpp +++ b/src/opticalFlow-sycl/derivativesKernel.hpp @@ -1,152 +1,77 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) + */ #include "common.h" -/////////////////////////////////////////////////////////////////////////////// -/// \brief compute image derivatives -/// -/// CUDA kernel, relies heavily on texture unit -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[out] Ix x derivative -/// \param[out] Iy y derivative -/// \param[out] Iz temporal derivative -/////////////////////////////////////////////////////////////////////////////// -void ComputeDerivativesKernel(int width, int height, int stride, float *Ix, - float *Iy, float *Iz, - sycl::accessor texSource, - sycl::accessor texTarget, - sycl::sampler texDesc, - const sycl::nd_item<3> &item) { +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(t,x,y) tex2D_sw(t,x,y) +#define TEX_CREATE(v,s,w,h,st) Tex2DDesc v{s,w,h,st} +#define TEX_DESTROY(v) + +void ComputeDerivativesKernel(int width, int height, int stride, + float *Ix, float *Iy, float *Iz, + TEX_TYPE texSource, TEX_TYPE texTarget, + const sycl::nd_item<3> &item) { const int ix = item.get_global_id(2); const int iy = item.get_global_id(1); - const int pos = ix + iy * stride; - if (ix >= width || iy >= height) return; - float t0, t1; - - auto x_inputCoords1 = sycl::float2(ix - 2.0f, iy); - auto x_inputCoords2 = sycl::float2(ix - 1.0f, iy); - auto x_inputCoords3 = sycl::float2(ix + 1.0f, iy); - auto x_inputCoords4 = sycl::float2(ix + 2.0f, iy); - - t0 = texSource.read(x_inputCoords1, texDesc)[0]; - t0 -= texSource.read(x_inputCoords2, texDesc)[0] * 8.0f; - t0 += texSource.read(x_inputCoords3, texDesc)[0] * 8.0f; - t0 -= texSource.read(x_inputCoords4, texDesc)[0]; - t0 /= 12.0f; - - t1 = texTarget.read(x_inputCoords1, texDesc)[0]; - t1 -= texTarget.read(x_inputCoords2, texDesc)[0] * 8.0f; - t1 += texTarget.read(x_inputCoords3, texDesc)[0] * 8.0f; - t1 -= texTarget.read(x_inputCoords4, texDesc)[0]; - t1 /= 12.0f; + float dx = 1.0f/(float)width, dy = 1.0f/(float)height; + float x = ((float)ix + 0.5f) * dx; + float y = ((float)iy + 0.5f) * dy; + float t0, t1; + t0 = TEX_FETCH(texSource, x-2.f*dx, y); + t0 -= TEX_FETCH(texSource, x-1.f*dx, y) * 8.f; + t0 += TEX_FETCH(texSource, x+1.f*dx, y) * 8.f; + t0 -= TEX_FETCH(texSource, x+2.f*dx, y); + t0 /= 12.f; + + t1 = TEX_FETCH(texTarget, x-2.f*dx, y); + t1 -= TEX_FETCH(texTarget, x-1.f*dx, y) * 8.f; + t1 += TEX_FETCH(texTarget, x+1.f*dx, y) * 8.f; + t1 -= TEX_FETCH(texTarget, x+2.f*dx, y); + t1 /= 12.f; Ix[pos] = (t0 + t1) * 0.5f; - // t derivative - auto inputCoord = sycl::float2(ix, iy); - Iz[pos] = texTarget.read(inputCoord, texDesc)[0] - - texSource.read(inputCoord, texDesc)[0]; - - // y derivative - auto y_inputCoords1 = sycl::float2(ix, iy - 2.0f); - auto y_inputCoords2 = sycl::float2(ix, iy - 1.0f); - auto y_inputCoords3 = sycl::float2(ix, iy + 1.0f); - auto y_inputCoords4 = sycl::float2(ix, iy + 2.0f); - - t0 = texSource.read(y_inputCoords1, texDesc)[0]; - t0 -= texSource.read(y_inputCoords2, texDesc)[0] * 8.0f; - t0 += texSource.read(y_inputCoords3, texDesc)[0] * 8.0f; - t0 -= texSource.read(y_inputCoords4, texDesc)[0]; - t0 /= 12.0f; + Iz[pos] = TEX_FETCH(texTarget, x, y) - TEX_FETCH(texSource, x, y); - t1 = texTarget.read(y_inputCoords1, texDesc)[0]; - t1 -= texTarget.read(y_inputCoords2, texDesc)[0] * 8.0f; - t1 += texTarget.read(y_inputCoords3, texDesc)[0] * 8.0f; - t1 -= texTarget.read(y_inputCoords4, texDesc)[0]; - t1 /= 12.0f; + t0 = TEX_FETCH(texSource, x, y-2.f*dy); + t0 -= TEX_FETCH(texSource, x, y-1.f*dy) * 8.f; + t0 += TEX_FETCH(texSource, x, y+1.f*dy) * 8.f; + t0 -= TEX_FETCH(texSource, x, y+2.f*dy); + t0 /= 12.f; + t1 = TEX_FETCH(texTarget, x, y-2.f*dy); + t1 -= TEX_FETCH(texTarget, x, y-1.f*dy) * 8.f; + t1 += TEX_FETCH(texTarget, x, y+1.f*dy) * 8.f; + t1 -= TEX_FETCH(texTarget, x, y+2.f*dy); + t1 /= 12.f; Iy[pos] = (t0 + t1) * 0.5f; } -/////////////////////////////////////////////////////////////////////////////// -/// \brief compute image derivatives -/// -/// \param[in] I0 source image -/// \param[in] I1 tracked image -/// \param[in] w image width -/// \param[in] h image height -/// \param[in] s image stride -/// \param[out] Ix x derivative -/// \param[out] Iy y derivative -/// \param[out] Iz temporal derivative -/////////////////////////////////////////////////////////////////////////////// -static void ComputeDerivatives(const float *I0, const float *I1, float *pI0_h, - float *pI1_h, float *I0_h, float *I1_h, - float *src_d0, float *src_d1, int w, int h, - int s, float *Ix, float *Iy, float *Iz, sycl::queue q) { +static void ComputeDerivatives(const float *I0, const float *I1, + int w, int h, int s, + float *Ix, float *Iy, float *Iz, + sycl::queue &q) +{ sycl::range<3> threads(1, 6, 32); sycl::range<3> blocks(1, iDivUp(h, threads[1]), iDivUp(w, threads[2])); - int dataSize = s * h * sizeof(float); - - q.memcpy(I0_h, I0, dataSize); - q.memcpy(I1_h, I1, dataSize); - q.wait(); - - for (int i = 0; i < h; i++) { - for (int j = 0; j < w; j++) { - int index = i * s + j; - pI0_h[index * 4 + 0] = I0_h[index]; - pI0_h[index * 4 + 1] = pI0_h[index * 4 + 2] = pI0_h[index * 4 + 3] = 0.f; - } - } - - for (int i = 0; i < h; i++) { - for (int j = 0; j < w; j++) { - int index = i * s + j; - pI1_h[index * 4 + 0] = I1_h[index]; - pI1_h[index * 4 + 1] = pI1_h[index * 4 + 2] = pI1_h[index * 4 + 3] = 0.f; - } - } - - q.memcpy(src_d0, pI0_h, s * h * sizeof(sycl::float4)); - q.memcpy(src_d1, pI1_h, s * h * sizeof(sycl::float4)); - - q.wait(); - - auto texDescr = sycl::sampler( - sycl::coordinate_normalization_mode::unnormalized, - sycl::addressing_mode::clamp_to_edge, sycl::filtering_mode::nearest); - - auto texSource = - sycl::image<2>(src_d0, sycl::image_channel_order::rgba, - sycl::image_channel_type::fp32, sycl::range<2>(w, h), - sycl::range<1>(s * sizeof(sycl::float4))); - - auto texTarget = - sycl::image<2>(src_d1, sycl::image_channel_order::rgba, - sycl::image_channel_type::fp32, sycl::range<2>(w, h), - sycl::range<1>(s * sizeof(sycl::float4))); - - - q.submit([&](sycl::handler &cgh) { - auto texSource_acc = - texSource.template get_access(cgh); - auto texTarget_acc = - texTarget.template get_access(cgh); - - cgh.parallel_for( - sycl::nd_range<3>(blocks * threads, threads), - [=](sycl::nd_item<3> item) { - ComputeDerivativesKernel( - w, h, s, Ix, Iy, Iz, - texSource_acc, texTarget_acc, - texDescr, item); - }); - }); + TEX_CREATE(texSource, I0, w, h, s); + TEX_CREATE(texTarget, I1, w, h, s); + q.parallel_for( + sycl::nd_range<3>(blocks * threads, threads), + [=](sycl::nd_item<3> item) { + ComputeDerivativesKernel(w, h, s, Ix, Iy, Iz, texSource, texTarget, + item); + }).wait(); + TEX_DESTROY(texSource); + TEX_DESTROY(texTarget); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-sycl/downscaleKernel.hpp b/src/opticalFlow-sycl/downscaleKernel.hpp index 38f8707a4e..f075a818ff 100644 --- a/src/opticalFlow-sycl/downscaleKernel.hpp +++ b/src/opticalFlow-sycl/downscaleKernel.hpp @@ -1,87 +1,49 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) + */ #include "common.h" -/////////////////////////////////////////////////////////////////////////////// -/// \brief downscale image -/// -/// CUDA kernel, relies heavily on texture unit -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// -void DownscaleKernel(int width, int height, int stride, float *out, - sycl::accessor tex_acc, - sycl::sampler texDesc, - const sycl::nd_item<3> &item) { +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(tex, x, y) tex2D_sw(tex, x, y) +#define TEX_CREATE(texVar, src, width, height, stride) \ + Tex2DDesc texVar { src, width, height, stride } +#define TEX_DESTROY(texVar) /* no-op: no object to destroy */ + +void DownscaleKernel(int width, int height, int stride, + float *out, TEX_TYPE texFine, + const sycl::nd_item<3> &item) { const int ix = item.get_global_id(2); const int iy = item.get_global_id(1); + if (ix >= width || iy >= height) return; - if (ix >= width || iy >= height) { - return; - } - - int srcx = ix * 2; - int srcy = iy * 2; - - auto inputCoords1 = sycl::float2(srcx + 0, srcy + 0); - auto inputCoords2 = sycl::float2(srcx + 0, srcy + 1); - auto inputCoords3 = sycl::float2(srcx + 1, srcy + 0); - auto inputCoords4 = sycl::float2(srcx + 1, srcy + 1); + float dx = 1.0f / (float)width; + float dy = 1.0f / (float)height; + float x = ((float)ix + 0.5f) * dx; + float y = ((float)iy + 0.5f) * dy; - out[ix + iy * stride] = 0.25f * (tex_acc.read(inputCoords1, texDesc)[0] + - tex_acc.read(inputCoords2, texDesc)[0] + - tex_acc.read(inputCoords3, texDesc)[0] + - tex_acc.read(inputCoords4, texDesc)[0]); + out[ix + iy * stride] = 0.25f * (TEX_FETCH(texFine, x - dx*0.25f, y) + + TEX_FETCH(texFine, x + dx*0.25f, y) + + TEX_FETCH(texFine, x, y - dy*0.25f) + + TEX_FETCH(texFine, x, y + dy*0.25f)); } -/////////////////////////////////////////////////////////////////////////////// -/// \brief downscale image -/// -/// \param[in] src image to downscale -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// -static void Downscale(const float *src, float *pI0_h, float *I0_h, float *src_p, int width, int height, int stride, - int newWidth, int newHeight, int newStride, float *out, sycl::queue q) { +static void Downscale(const float *src, int width, int height, int stride, + int newWidth, int newHeight, int newStride, float *out, + sycl::queue &q) +{ sycl::range<3> threads(1, 8, 32); sycl::range<3> blocks(1, iDivUp(newHeight, threads[1]), iDivUp(newWidth, threads[2])); + TEX_CREATE(texFine, src, width, height, stride); - int dataSize = height * stride * sizeof(float); - - q.memcpy(I0_h, src, dataSize).wait(); - - for (int i = 0; i < height; i++) { - for (int j = 0; j < width; j++) { - int index = i * stride + j; - pI0_h[index * 4 + 0] = I0_h[index]; - pI0_h[index * 4 + 1] = pI0_h[index * 4 + 2] = pI0_h[index * 4 + 3] = 0.f; - } - } - - q.memcpy(src_p, pI0_h, height * width * sizeof(sycl::float4)).wait(); - - auto texDescr = sycl::sampler( - sycl::coordinate_normalization_mode::unnormalized, - sycl::addressing_mode::clamp_to_edge, sycl::filtering_mode::nearest); - - auto texFine = sycl::image<2>(src_p, sycl::image_channel_order::rgba, - sycl::image_channel_type::fp32, - sycl::range<2>(width, height), - sycl::range<1>(stride * sizeof(sycl::float4))); - - q.submit([&](sycl::handler &cgh) { - auto tex_acc = - texFine.template get_access(cgh); - - cgh.parallel_for(sycl::nd_range<3>(blocks * threads, threads), - [=](sycl::nd_item<3> item) { - DownscaleKernel(newWidth, newHeight, newStride, out, - tex_acc, texDescr, item); - }); - }); + q.parallel_for(sycl::nd_range<3>(blocks * threads, threads), + [=](sycl::nd_item<3> item) { + DownscaleKernel(newWidth, newHeight, newStride, out, texFine, item); + }).wait(); + TEX_DESTROY(texFine); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-sycl/flowGold.cpp b/src/opticalFlow-sycl/flowGold.cpp index aa292441e5..9ebc4eba52 100644 --- a/src/opticalFlow-sycl/flowGold.cpp +++ b/src/opticalFlow-sycl/flowGold.cpp @@ -1,3 +1,30 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * * Neither the name of NVIDIA CORPORATION nor the names of its + * contributors may be used to endorse or promote products derived + * from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY + * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY + * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + #include "common.h" #include "flowGold.h" #include @@ -19,9 +46,6 @@ inline float Tex2D(const float *t, int w, int h, int s, float x, float y) { // integer parts in floating point format float intPartX, intPartY; - x -= 0.5f; - y -= 0.5f; - // get fractional parts of coordinates float dx = fabsf(modff(x, &intPartX)); float dy = fabsf(modff(y, &intPartY)); @@ -32,21 +56,21 @@ inline float Tex2D(const float *t, int w, int h, int s, float x, float y) { int iy0 = (int)intPartY; // mirror out-of-range position - if (ix0 < 0) ix0 = 0; + if (ix0 < 0) ix0 = abs(ix0 + 1); - if (iy0 < 0) iy0 = 0; + if (iy0 < 0) iy0 = abs(iy0 + 1); - if (ix0 >= w) ix0 = w - 1; + if (ix0 >= w) ix0 = w * 2 - ix0 - 1; - if (iy0 >= h) iy0 = h - 1; + if (iy0 >= h) iy0 = h * 2 - iy0 - 1; // corner which is opposite to (ix0, iy0) int ix1 = ix0 + 1; int iy1 = iy0 + 1; - if (ix1 >= w) ix1 = w - 1; + if (ix1 >= w) ix1 = w * 2 - ix1 - 1; - if (iy1 >= h) iy1 = h - 1; + if (iy1 >= h) iy1 = h * 2 - iy1 - 1; float res = t[ix0 + iy0 * s] * (1.0f - dx) * (1.0f - dy); res += t[ix1 + iy0 * s] * dx * (1.0f - dy); @@ -70,13 +94,13 @@ inline float Tex2D(const float *t, int w, int h, int s, float x, float y) { /// \return fetched value /////////////////////////////////////////////////////////////////////////////// inline float Tex2Di(const float *src, int w, int h, int s, int x, int y) { - if (x < 0) x = 0; + if (x < 0) x = abs(x + 1); - if (y < 0) y = 0; + if (y < 0) y = abs(y + 1); - if (x >= w) x = w - 1; + if (x >= w) x = w * 2 - x - 1; - if (y >= h) y = h - 1; + if (y >= h) y = h * 2 - y - 1; return src[x + y * s]; } diff --git a/src/opticalFlow-sycl/flowGold.h b/src/opticalFlow-sycl/flowGold.h index 8337af3736..0dc0297019 100644 --- a/src/opticalFlow-sycl/flowGold.h +++ b/src/opticalFlow-sycl/flowGold.h @@ -1,9 +1,3 @@ -//========================================================= -// Modifications Copyright © 2022 Intel Corporation -// -// SPDX-License-Identifier: BSD-3-Clause -//========================================================= - /* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. * * Redistribution and use in source and binary forms, with or without diff --git a/src/opticalFlow-sycl/flowSYCL.cpp b/src/opticalFlow-sycl/flowSYCL.cpp index 4a9731728d..dcd4a5d0ba 100644 --- a/src/opticalFlow-sycl/flowSYCL.cpp +++ b/src/opticalFlow-sycl/flowSYCL.cpp @@ -1,3 +1,30 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * * Neither the name of NVIDIA CORPORATION nor the names of its + * contributors may be used to endorse or promote products derived + * from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY + * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY + * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + #include "common.h" // include kernels @@ -8,7 +35,6 @@ #include "solverKernel.hpp" #include "addKernel.hpp" - /////////////////////////////////////////////////////////////////////////////// /// \brief method logic /// @@ -25,15 +51,17 @@ /// \param[out] u horizontal displacement /// \param[out] v vertical displacement /////////////////////////////////////////////////////////////////////////////// -void ComputeFlowSYCL(const float *I0, const float *I1, int width, int height, +void ComputeFlowCUDA(const float *I0, const float *I1, int width, int height, int stride, float alpha, int nLevels, int nWarpIters, int nSolverIters, float *u, float *v) { - printf("Computing optical flow on Device...\n"); - sycl::queue q{aspect_selector(sycl::aspect::ext_intel_legacy_image), sycl::property::queue::in_order()}; +#ifdef USE_GPU + sycl::queue q(sycl::gpu_selector_v, sycl::property::queue::in_order()); +#else + sycl::queue q(sycl::cpu_selector_v, sycl::property::queue::in_order()); +#endif + printf("Computing optical flow on Device...\n"); - std::cout << "\nRunning on " - << q.get_device().get_info() << "\n"; // pI0 and pI1 will hold device pointers const float **pI0 = new const float *[nLevels]; const float **pI1 = new const float *[nLevels]; @@ -82,22 +110,8 @@ void ComputeFlowSYCL(const float *I0, const float *I1, int width, int height, CHECK_ERROR(*(pI0 + currentLevel) = (const float *)sycl::malloc_device(dataSize, q)); CHECK_ERROR(*(pI1 + currentLevel) = (const float *)sycl::malloc_device(dataSize, q)); - float *pI0_h = (float *)sycl::malloc_host(stride * height * sizeof(sycl::float4), q); - float *I0_h = (float *)sycl::malloc_host(dataSize, q); - - float *pI1_h = (float *)sycl::malloc_host(stride * height * sizeof(sycl::float4), q); - float *I1_h = (float *)sycl::malloc_host(dataSize, q); - - float *src_d0 = (float *)sycl::malloc_device(stride * height * sizeof(sycl::float4), q); - float *src_d1 = (float *)sycl::malloc_device(stride * height * sizeof(sycl::float4), q); - - q.memcpy((void *)I0_h, I0, dataSize); - q.memcpy((void *)I1_h, I1, dataSize); - - q.memcpy((void *)pI0[currentLevel], I0, dataSize); - q.memcpy((void *)pI1[currentLevel], I1, dataSize); - - q.wait(); + CHECK_ERROR(q.memcpy((void *)pI0[currentLevel], I0, dataSize)); + CHECK_ERROR(q.memcpy((void *)pI1[currentLevel], I1, dataSize)); pW[currentLevel] = width; pH[currentLevel] = height; @@ -108,20 +122,16 @@ void ComputeFlowSYCL(const float *I0, const float *I1, int width, int height, int nh = pH[currentLevel] / 2; int ns = iAlignUp(nw); - CHECK_ERROR( - *(pI0 + currentLevel - 1) = (const float *)sycl::malloc_device( - ns * nh * sizeof(float), q)); - CHECK_ERROR( - *(pI1 + currentLevel - 1) = (const float *)sycl::malloc_device( - ns * nh * sizeof(float), q)); + CHECK_ERROR(*(pI0 + currentLevel - 1) = (const float *)sycl::malloc_device( + ns * nh * sizeof(float), q)); + CHECK_ERROR(*(pI1 + currentLevel - 1) = (const float *)sycl::malloc_device( + ns * nh * sizeof(float), q)); - Downscale(pI0[currentLevel], pI0_h, I0_h, src_d0, pW[currentLevel], - pH[currentLevel], pS[currentLevel], nw, nh, ns, - (float *)pI0[currentLevel - 1], q); + Downscale(pI0[currentLevel], pW[currentLevel], pH[currentLevel], + pS[currentLevel], nw, nh, ns, (float *)pI0[currentLevel - 1], q); - Downscale(pI1[currentLevel], pI0_h, I0_h, src_d0, pW[currentLevel], - pH[currentLevel], pS[currentLevel], nw, nh, ns, - (float *)pI1[currentLevel - 1], q); + Downscale(pI1[currentLevel], pW[currentLevel], pH[currentLevel], + pS[currentLevel], nw, nh, ns, (float *)pI1[currentLevel - 1], q); pW[currentLevel - 1] = nw; pH[currentLevel - 1] = nh; @@ -143,11 +153,10 @@ void ComputeFlowSYCL(const float *I0, const float *I1, int width, int height, // on current level we compute optical flow // between frame 0 and warped frame 1 - WarpImage(pI1[currentLevel], pI0_h, I0_h, src_d0, pW[currentLevel], pH[currentLevel], + WarpImage(pI1[currentLevel], pW[currentLevel], pH[currentLevel], pS[currentLevel], d_u, d_v, d_tmp, q); - ComputeDerivatives(pI0[currentLevel], d_tmp, pI0_h, pI1_h, I0_h, I1_h, - src_d0, src_d1, pW[currentLevel], + ComputeDerivatives(pI0[currentLevel], d_tmp, pW[currentLevel], pH[currentLevel], pS[currentLevel], d_Ix, d_Iy, d_Iz, q); for (int iter = 0; iter < nSolverIters; ++iter) { @@ -167,13 +176,13 @@ void ComputeFlowSYCL(const float *I0, const float *I1, int width, int height, // prolongate solution float scaleX = (float)pW[currentLevel + 1] / (float)pW[currentLevel]; - Upscale(d_u, pI0_h, I0_h, src_d0, pW[currentLevel], pH[currentLevel], pS[currentLevel], + Upscale(d_u, pW[currentLevel], pH[currentLevel], pS[currentLevel], pW[currentLevel + 1], pH[currentLevel + 1], pS[currentLevel + 1], scaleX, d_nu, q); float scaleY = (float)pH[currentLevel + 1] / (float)pH[currentLevel]; - Upscale(d_v, pI0_h, I0_h, src_d0, pW[currentLevel], pH[currentLevel], pS[currentLevel], + Upscale(d_v, pW[currentLevel], pH[currentLevel], pS[currentLevel], pW[currentLevel + 1], pH[currentLevel + 1], pS[currentLevel + 1], scaleY, d_nv, q); diff --git a/src/opticalFlow-sycl/flowSYCL.h b/src/opticalFlow-sycl/flowSYCL.h index 7c902f38d6..ec727901b4 100644 --- a/src/opticalFlow-sycl/flowSYCL.h +++ b/src/opticalFlow-sycl/flowSYCL.h @@ -1,9 +1,3 @@ -//========================================================= -// Modifications Copyright © 2022 Intel Corporation -// -// SPDX-License-Identifier: BSD-3-Clause -//========================================================= - /* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -34,7 +28,7 @@ #ifndef FLOW_SYCL_H #define FLOW_SYCL_H -void ComputeFlowSYCL( +void ComputeFlowCUDA( const float *I0, // source frame const float *I1, // tracked frame int width, // frame width diff --git a/src/opticalFlow-sycl/main.cpp b/src/opticalFlow-sycl/main.cpp index 83e8f7d722..5dcd3dff2f 100644 --- a/src/opticalFlow-sycl/main.cpp +++ b/src/opticalFlow-sycl/main.cpp @@ -1,13 +1,38 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * * Neither the name of NVIDIA CORPORATION nor the names of its + * contributors may be used to endorse or promote products derived + * from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY + * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY + * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + const static char *const sSDKsample = "HSOpticalFlow"; // CPU-GPU discrepancy threshold for self-test const float THRESHOLD = 0.05f; - #include "common.h" #include "flowGold.h" #include "flowSYCL.h" - #include "helper_functions.h" #include #include @@ -125,7 +150,7 @@ bool CompareWithGold(int width, int height, int stride, const float *h_uGold, printf("L1 error : %.6f\n", error); - return (error < 1.0f); + return (error < THRESHOLD); } /////////////////////////////////////////////////////////////////////////////// @@ -189,18 +214,20 @@ int main(int argc, char **argv) { auto duration = std::chrono::duration_cast(stop - start).count(); printf("Processing time on CPU: %f (ms)\n", duration); - + start = Time::now(); - ComputeFlowSYCL(h_source, h_target, width, height, stride, alpha, nLevels, + ComputeFlowCUDA(h_source, h_target, width, height, stride, alpha, nLevels, nWarpIters, nSolverIters, h_u, h_v); stop = Time::now(); duration = std::chrono::duration_cast(stop - start).count(); printf("Processing time on Device: %f (ms)\n", duration); // compare results (L1 norm) - bool status = CompareWithGold(width, height, stride, h_uGold, h_vGold, h_u, h_v); + bool status = + CompareWithGold(width, height, stride, h_uGold, h_vGold, h_u, h_v); WriteFloFile("FlowGPU.flo", width, height, stride, h_u, h_v); + WriteFloFile("FlowCPU.flo", width, height, stride, h_uGold, h_vGold); // free resources diff --git a/src/opticalFlow-sycl/solverKernel.hpp b/src/opticalFlow-sycl/solverKernel.hpp index f95219b92d..87940d5bfb 100644 --- a/src/opticalFlow-sycl/solverKernel.hpp +++ b/src/opticalFlow-sycl/solverKernel.hpp @@ -1,7 +1,34 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * * Neither the name of NVIDIA CORPORATION nor the names of its + * contributors may be used to endorse or promote products derived + * from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY + * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY + * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + #include "common.h" /////////////////////////////////////////////////////////////////////////////// -/// \brief one iteration of classical Horn-Schunck method, CUDA kernel. +/// \brief one iteration of classical Horn-Schunck method /// /// It is one iteration of Jacobi method for a corresponding linear system. /// Template parameters are describe CTA size @@ -24,9 +51,6 @@ void JacobiIteration(const float *du0, const float *dv0, float alpha, float *du1, float *dv1, const sycl::nd_item<3> &item, volatile float *du, volatile float *dv) { - // Handle to thread block group - auto cta = item.get_group(); - const int ix = item.get_global_id(2); const int iy = item.get_global_id(1); @@ -140,24 +164,26 @@ void JacobiIteration(const float *du0, const float *dv0, /////////////////////////////////////////////////////////////////////////////// static void SolveForUpdate(const float *du0, const float *dv0, const float *Ix, const float *Iy, const float *Iz, int w, int h, - int s, float alpha, float *du1, float *dv1, sycl::queue &q) { + int s, float alpha, float *du1, float *dv1, + sycl::queue &q) { // CTA size sycl::range<3> threads(1, 6, 32); // grid size sycl::range<3> blocks(1, iDivUp(h, threads[1]), iDivUp(w, threads[2])); q.submit([&](sycl::handler &cgh) { - sycl::local_accessor du_acc_ct1( + sycl::local_accessor du_acc( sycl::range<1>((32 + 2) * (6 + 2)), cgh); - sycl::local_accessor dv_acc_ct1( + sycl::local_accessor dv_acc( sycl::range<1>((32 + 2) * (6 + 2)), cgh); - cgh.parallel_for(sycl::nd_range<3>(blocks * threads, threads), - [=](sycl::nd_item<3> item) { - JacobiIteration<32, 6>(du0, dv0, Ix, Iy, Iz, w, h, s, - alpha, du1, dv1, item, - du_acc_ct1.get_multi_ptr().get(), - dv_acc_ct1.get_multi_ptr().get()); - }); + cgh.parallel_for( + sycl::nd_range<3>(blocks * threads, threads), + [=](sycl::nd_item<3> item) { + JacobiIteration<32, 6>( + du0, dv0, Ix, Iy, Iz, w, h, s, alpha, du1, dv1, item, + du_acc.get_multi_ptr().get(), + dv_acc.get_multi_ptr().get()); + }); }); } diff --git a/src/opticalFlow-sycl/tex2d_compat.h b/src/opticalFlow-sycl/tex2d_compat.h new file mode 100644 index 0000000000..4d8ef19421 --- /dev/null +++ b/src/opticalFlow-sycl/tex2d_compat.h @@ -0,0 +1,41 @@ +// tex2d_compat.h — Portable tex2D + +#pragma once + +// Descriptor passed to kernels in place of hipTextureObject_t +struct Tex2DDesc { + const float* ptr; // pitched device pointer + int width; // image width in elements + int height; // image height in elements + int pitch; // row stride in elements (pitchInBytes / sizeof(float)) +}; + +// Mirror repeat address clamping +inline int tex2d_mirror_clamp(int i, int size) { + int mi = ((i % (2 * size)) + 2 * size) % (2 * size); + return (mi >= size) ? (2 * size - 1 - mi) : mi; +} + +// Software bilinear tex2D — normalized coords [0,1), mirror addressing, linear filter +// Matches: normalizedCoords=true, filterMode=Linear, addressMode=Mirror +inline float tex2D_sw(const Tex2DDesc& tex, float u, float v) { + // Convert normalized coords to pixel space + float fx = u * (float)tex.width - 0.5f; + float fy = v * (float)tex.height - 0.5f; + + int x0 = (int)sycl::floor(fx), y0 = (int)sycl::floor(fy); + float tx = fx - (float)x0, ty = fy - (float)y0; + + int cx0 = tex2d_mirror_clamp(x0, tex.width); + int cx1 = tex2d_mirror_clamp(x0+1, tex.width); + int cy0 = tex2d_mirror_clamp(y0, tex.height); + int cy1 = tex2d_mirror_clamp(y0+1, tex.height); + + float v00 = tex.ptr[cy0 * tex.pitch + cx0]; + float v10 = tex.ptr[cy0 * tex.pitch + cx1]; + float v01 = tex.ptr[cy1 * tex.pitch + cx0]; + float v11 = tex.ptr[cy1 * tex.pitch + cx1]; + + return (1.f-tx)*(1.f-ty)*v00 + tx*(1.f-ty)*v10 + + (1.f-tx)*ty *v01 + tx*ty *v11; +} diff --git a/src/opticalFlow-sycl/upscaleKernel.hpp b/src/opticalFlow-sycl/upscaleKernel.hpp index 24eb7a5359..7a6c330f19 100644 --- a/src/opticalFlow-sycl/upscaleKernel.hpp +++ b/src/opticalFlow-sycl/upscaleKernel.hpp @@ -1,82 +1,39 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) + */ #include "common.h" - -/////////////////////////////////////////////////////////////////////////////// -/// \brief upscale one component of a displacement field, CUDA kernel -/// \param[in] width field width -/// \param[in] height field height -/// \param[in] stride field stride -/// \param[in] scale scale factor (multiplier) -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// -void UpscaleKernel(int width, int height, int stride, float scale, float *out, - sycl::accessor texCoarse_acc, - sycl::sampler texDesc, - const sycl::nd_item<3> &item) { +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(t,x,y) tex2D_sw(t,x,y) +#define TEX_CREATE(v,s,w,h,st) Tex2DDesc v{s,w,h,st} +#define TEX_DESTROY(v) + +void UpscaleKernel(int width, int height, int stride, + float scale, float *out, TEX_TYPE texCoarse, + const sycl::nd_item<3> &item) { const int ix = item.get_global_id(2); const int iy = item.get_global_id(1); - if (ix >= width || iy >= height) return; - - float x = ((float)ix - 0.5f) * 0.5f; - float y = ((float)iy - 0.5f) * 0.5f; - - auto inputCoord = sycl::float2(x, y); - - // exploit hardware interpolation - // and scale interpolated vector to match next pyramid level resolution - out[ix + iy * stride] = texCoarse_acc.read(inputCoord, texDesc)[0] * scale; + float x = ((float)ix + 0.5f) / (float)width; + float y = ((float)iy + 0.5f) / (float)height; + out[ix + iy * stride] = TEX_FETCH(texCoarse, x, y) * scale; } -/////////////////////////////////////////////////////////////////////////////// -/// \brief upscale one component of a displacement field, kernel wrapper -/// \param[in] src field component to upscale -/// \param[in] width field current width -/// \param[in] height field current height -/// \param[in] stride field current stride -/// \param[in] newWidth field new width -/// \param[in] newHeight field new height -/// \param[in] newStride field new stride -/// \param[in] scale value scale factor (multiplier) -/// \param[out] out upscaled field component -/////////////////////////////////////////////////////////////////////////////// -static void Upscale(const float *src, float *pI0_h, float *I0_h, float *src_p, int width, int height, int stride, - int newWidth, int newHeight, int newStride, float scale, - float *out, sycl::queue &q) { +static void Upscale(const float *src, int width, int height, int stride, + int newWidth, int newHeight, int newStride, + float scale, float *out, sycl::queue &q) { sycl::range<3> threads(1, 8, 32); sycl::range<3> blocks(1, iDivUp(newHeight, threads[1]), iDivUp(newWidth, threads[2])); - - int dataSize = stride * height * sizeof(float); - q.memcpy(I0_h, src, dataSize).wait(); - - for (int i = 0; i < height; i++) { - for (int j = 0; j < width; j++) { - int index = i * stride + j; - pI0_h[index * 4 + 0] = I0_h[index]; - pI0_h[index * 4 + 1] = pI0_h[index * 4 + 2] = pI0_h[index * 4 + 3] = 0.f; - } - } - q.memcpy(src_p, pI0_h, height * stride * sizeof(sycl::float4)).wait(); - - auto texDescr = sycl::sampler( - sycl::coordinate_normalization_mode::unnormalized, - sycl::addressing_mode::clamp_to_edge, sycl::filtering_mode::linear); - - auto texCoarse = sycl::image<2>( - src_p, sycl::image_channel_order::rgba, - sycl::image_channel_type::fp32, sycl::range<2>(width, height), - sycl::range<1>(stride * sizeof(sycl::float4))); - - q.submit([&](sycl::handler &cgh) { - auto texCoarse_acc = - texCoarse.template get_access(cgh); - - cgh.parallel_for(sycl::nd_range<3>(blocks * threads, threads), - [=](sycl::nd_item<3> item) { - UpscaleKernel(newWidth, newHeight, newStride, scale, out, - texCoarse_acc, texDescr, item); - }); - }); + TEX_CREATE(texCoarse, src, width, height, stride); + q.parallel_for( + sycl::nd_range<3>(blocks * threads, threads), + [=](sycl::nd_item<3> item) { + UpscaleKernel(newWidth, newHeight, newStride, scale, out, texCoarse, item); + }).wait(); + TEX_DESTROY(texCoarse); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/opticalFlow-sycl/warpingKernel.hpp b/src/opticalFlow-sycl/warpingKernel.hpp index d553f8f7d2..680447825a 100644 --- a/src/opticalFlow-sycl/warpingKernel.hpp +++ b/src/opticalFlow-sycl/warpingKernel.hpp @@ -1,87 +1,39 @@ +/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. (BSD license) + */ #include "common.h" - -/////////////////////////////////////////////////////////////////////////////// -/// \brief warp image with a given displacement field, CUDA kernel. -/// \param[in] width image width -/// \param[in] height image height -/// \param[in] stride image stride -/// \param[in] u horizontal displacement -/// \param[in] v vertical displacement -/// \param[out] out result -/////////////////////////////////////////////////////////////////////////////// -void WarpingKernel(int width, int height, int stride, const float *u, - const float *v, float *out, - sycl::accessor - texToWarp, - sycl::sampler texDesc, - const sycl::nd_item<3> &item) { +#include "tex2d_compat.h" +#define TEX_TYPE Tex2DDesc +#define TEX_FETCH(t,x,y) tex2D_sw(t,x,y) +#define TEX_CREATE(v,s,w,h,st) Tex2DDesc v{s,w,h,st} +#define TEX_DESTROY(v) + +void WarpingKernel(int width, int height, int stride, + const float *u, const float *v, + float *out, TEX_TYPE texToWarp, + const sycl::nd_item<3> &item) { const int ix = item.get_global_id(2); const int iy = item.get_global_id(1); - const int pos = ix + iy * stride; - if (ix >= width || iy >= height) return; - - float x = ((float)ix + u[pos]); - float y = ((float)iy + v[pos]); - - auto inputCoord = sycl::float2(x, y); - - out[pos] = texToWarp.read(inputCoord, texDesc)[0]; + float x = ((float)ix + u[pos] + 0.5f) / (float)width; + float y = ((float)iy + v[pos] + 0.5f) / (float)height; + out[pos] = TEX_FETCH(texToWarp, x, y); } -/////////////////////////////////////////////////////////////////////////////// -/// \brief warp image with provided vector field, CUDA kernel wrapper. -/// -/// For each output pixel there is a vector which tells which pixel -/// from a source image should be mapped to this particular output -/// pixel. -/// It is assumed that images and the vector field have the same stride and -/// resolution. -/// \param[in] src source image -/// \param[in] w width -/// \param[in] h height -/// \param[in] s stride -/// \param[in] u horizontal displacement -/// \param[in] v vertical displacement -/// \param[out] out warped image -/////////////////////////////////////////////////////////////////////////////// -static void WarpImage(const float *src, float *pI0_h, float *I0_h, float *src_p, int w, int h, int s, const float *u, - const float *v, float *out, sycl::queue &q) { +static void WarpImage(const float *src, int w, int h, int s, + const float *u, const float *v, float *out, sycl::queue &q) { sycl::range<3> threads(1, 6, 32); sycl::range<3> blocks(1, iDivUp(h, threads[1]), iDivUp(w, threads[2])); - - int dataSize = s * h * sizeof(float); - q.memcpy(I0_h, src, dataSize).wait(); - - for (int i = 0; i < h; i++) { - for (int j = 0; j < w; j++) { - int index = i * s + j; - pI0_h[index * 4 + 0] = I0_h[index]; - pI0_h[index * 4 + 1] = pI0_h[index * 4 + 2] = pI0_h[index * 4 + 3] = 0.f; - } - } - q.memcpy(src_p, pI0_h, s * h * sizeof(sycl::float4)).wait(); - - auto texDescr = sycl::sampler( - sycl::coordinate_normalization_mode::unnormalized, - sycl::addressing_mode::clamp_to_edge, sycl::filtering_mode::linear); - - auto texToWarp = - sycl::image<2>(src_p, sycl::image_channel_order::rgba, - sycl::image_channel_type::fp32, sycl::range<2>(w, h), - sycl::range<1>(s * sizeof(sycl::float4))); - - q.submit([&](sycl::handler &cgh) { - auto texToWarp_acc = - texToWarp.template get_access(cgh); - - cgh.parallel_for(sycl::nd_range<3>(blocks * threads, threads), - [=](sycl::nd_item<3> item) { - WarpingKernel(w, h, s, u, v, out, - texToWarp_acc, texDescr, item); - }); - }); + TEX_CREATE(texToWarp, src, w, h, s); + q.parallel_for( + sycl::nd_range<3>(blocks * threads, threads), + [=](sycl::nd_item<3> item) { + WarpingKernel(w, h, s, u, v, out, texToWarp, item); + }).wait(); + TEX_DESTROY(texToWarp); } + +#undef TEX_TYPE +#undef TEX_FETCH +#undef TEX_CREATE +#undef TEX_DESTROY diff --git a/src/softmax-cuda/main.cu b/src/softmax-cuda/main.cu index 5bb28c61d0..8c058210b8 100644 --- a/src/softmax-cuda/main.cu +++ b/src/softmax-cuda/main.cu @@ -26,6 +26,7 @@ void softMax_cpu(const int numSlice, const int sliceSize, const float* src, floa } } +__launch_bounds__(BLOCK_SIZE) __global__ void softMax (const int numSlice, const int sliceSize, const float* src, float* dest) @@ -45,6 +46,7 @@ void softMax (const int numSlice, const int sliceSize, } } +__launch_bounds__(BLOCK_SIZE) __global__ void softMax2 (const int numSlice, const int sliceSize, const float* src, float* dest) diff --git a/src/softmax-hip/main.cu b/src/softmax-hip/main.cu index 67b0969c3d..e2b141ff28 100644 --- a/src/softmax-hip/main.cu +++ b/src/softmax-hip/main.cu @@ -26,6 +26,7 @@ void softMax_cpu(const int numSlice, const int sliceSize, const float* src, floa } } +__launch_bounds__(BLOCK_SIZE) __global__ void softMax (const int numSlice, const int sliceSize, const float* src, float* dest) @@ -45,6 +46,7 @@ void softMax (const int numSlice, const int sliceSize, } } +__launch_bounds__(BLOCK_SIZE) __global__ void softMax2 (const int numSlice, const int sliceSize, const float* src, float* dest) diff --git a/src/softmax-sycl/main.cpp b/src/softmax-sycl/main.cpp index 7d866d440b..66bbe05f40 100644 --- a/src/softmax-sycl/main.cpp +++ b/src/softmax-sycl/main.cpp @@ -73,7 +73,8 @@ int main(int argc, char* argv[]) { for (int n = 0; n < repeat; n++) { q.submit([&](sycl::handler &h) { h.parallel_for( - sycl::nd_range<1>(gws, lws), [=](sycl::nd_item<1> item) + sycl::nd_range<1>(gws, lws), + [=] [[sycl::reqd_work_group_size(BLOCK_SIZE)]] (sycl::nd_item<1> item) //[[sycl::reqd_sub_group_size(warpSize)]] { sycl::sub_group warp = item.get_sub_group(); @@ -111,7 +112,8 @@ int main(int argc, char* argv[]) { for (int n = 0; n < repeat; n++) { q.submit([&](sycl::handler &h) { h.parallel_for( - sycl::nd_range<1>(gws, lws), [=](sycl::nd_item<1> item) { + sycl::nd_range<1>(gws, lws), + [=] [[sycl::reqd_work_group_size(BLOCK_SIZE)]] (sycl::nd_item<1> item) { int i = item.get_global_id(0); if (i >= numSlice) return; float max_ = d_input[i * sliceSize];