SHOGUN  4.0.0
Convolve.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) The Shogun Machine Learning Toolbox
3  * Written (w) 2014 Khaled Nasr
4  * All rights reserved.
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions are met:
8  *
9  * 1. Redistributions of source code must retain the above copyright notice, this
10  * list of conditions and the following disclaimer.
11  * 2. Redistributions in binary form must reproduce the above copyright notice,
12  * this list of conditions and the following disclaimer in the documentation
13  * and/or other materials provided with the distribution.
14  *
15  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
19  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25  *
26  * The views and conclusions contained in the software and documentation are those
27  * of the authors and should not be interpreted as representing official policies,
28  * either expressed or implied, of the Shogun Development Team.
29  */
30 
31 #ifndef CONVOLVE_H_
32 #define CONVOLVE_H_
33 
34 #include <shogun/lib/config.h>
35 #include <shogun/lib/SGMatrix.h>
37 
38 #include <shogun/io/SGIO.h>
39 
40 #ifdef HAVE_EIGEN3
42 #endif // HAVE_EIGEN3
43 
44 #ifdef HAVE_VIENNACL
45 #include <shogun/lib/GPUMatrix.h>
47 #endif // HAVE_VIENNACL
48 
49 namespace shogun
50 {
51 
52 namespace linalg
53 {
54 
55 namespace implementation
56 {
57 
61 template <enum Backend, class Matrix>
62 struct convolve
63 {
64  typedef typename Matrix::Scalar T;
65 
82  static void compute(Matrix X, Matrix W, Matrix Y, bool flip ,
83  bool overwrite, int32_t stride_x, int32_t stride_y);
84 };
85 
86 #ifdef HAVE_EIGEN3
87 
89 template <> template <class Matrix>
90 struct convolve<Backend::EIGEN3, Matrix>
91 {
92  typedef typename Matrix::Scalar T;
95 
109  static void compute(SGMatrix<T> X, SGMatrix<T> W, SGMatrix<T> Y, bool flip ,
110  bool overwrite, int32_t stride_x, int32_t stride_y)
111  {
112  int32_t width = X.num_cols;
113  int32_t height = X.num_rows;
114 
115  int32_t kx = W.num_cols;
116  int32_t ky = W.num_rows;
117 
118  int32_t rx = (kx-1)/2;
119  int32_t ry = (ky-1)/2;
120 
121  for (int32_t x=0; x<width; x+=stride_x)
122  {
123  int32_t xout = x/stride_x;
124 
125  for (int32_t y=0; y<height; y+=stride_y)
126  {
127  int32_t yout = y/stride_y;
128 
129  T sum = overwrite ? 0 : Y(yout,xout);
130  for (int32_t x1=x-rx; x1<=x+rx; x1++)
131  {
132  int32_t wx = flip ? x1-x+rx : rx-x1+x;
133  for (int32_t y1=y-ry; y1<=y+ry; y1++)
134  {
135  if (x1>=0 && y1>=0 && x1<width && y1<height)
136  {
137  if (flip)
138  sum += W(y1-y+ry,wx)*X(y1,x1);
139  else
140  sum += W(ry-y1+y,wx)*X(y1,x1);
141  }
142  }
143  }
144  Y(yout,xout) = sum;
145  }
146  }
147  }
148 };
149 #endif // HAVE_EIGEN3
150 
151 #ifdef HAVE_VIENNACL
152 
154 template <> template <class Matrix>
155 struct convolve<Backend::VIENNACL, Matrix>
156 {
157  typedef typename Matrix::Scalar T;
158 
160  template <class T>
161  static viennacl::ocl::kernel& generate_kernel_unity_stride(
162  int32_t radius_x, int32_t radius_y, bool flip, bool overwrite)
163  {
164  std::string kernel_name =
165  "convolve_unity_stride_" + ocl::get_type_string<T>() + "_" +
166  std::to_string(radius_x) + "_" + std::to_string(radius_y);
167 
168  if (flip) kernel_name.append("_flip");
169  if (overwrite) kernel_name.append("_overwrite");
170 
171  if (ocl::kernel_exists(kernel_name))
172  return ocl::get_kernel(kernel_name);
173 
174  std::string source = ocl::generate_kernel_preamble<T>(kernel_name);
175 
176  if (flip) source.append("#define FLIP\n");
177  if (overwrite) source.append("#define OVERWRITE\n");
178 
179  source.append("#define RADIUS_X " + std::to_string(radius_x) + "\n");
180  source.append("#define RADIUS_Y " + std::to_string(radius_y) + "\n");
181 
182  source.append(
183  R"(
184  #define W_WIDTH (2*RADIUS_X+1)
185  #define W_HEIGHT (2*RADIUS_Y+1)
186 
187  #define X_LOCAL_WIDTH (WORK_GROUP_SIZE_2D+2*RADIUS_X)
188  #define X_LOCAL_HEIGHT (WORK_GROUP_SIZE_2D+2*RADIUS_Y)
189 
190  inline DATATYPE readX(read_only __global DATATYPE* X, int x, int y,
191  int X_width, int X_height, int X_offset)
192  {
193  if (x>=0 && y>=0 && x<X_width && y<X_height)
194  return X[y + x*X_height + X_offset];
195  else
196  return 0;
197  }
198 
199  __kernel void KERNEL_NAME(
200  read_only __global DATATYPE* X, int X_width, int X_height, int X_offset,
201  __constant DATATYPE* W, int W_offset,
202  __global DATATYPE* Y, int Y_offset)
203  {
204  __local DATATYPE X_local[X_LOCAL_WIDTH][X_LOCAL_HEIGHT];
205 
206  int x = get_global_id(0);
207  int y = get_global_id(1);
208 
209  int xl = get_local_id(0);
210  int yl = get_local_id(1);
211 
212  if (xl==WORK_GROUP_SIZE_2D-1 && yl == WORK_GROUP_SIZE_2D-1)
213  {
214  for (int rx=0; rx<=2*RADIUS_X; rx++)
215  for (int ry=0; ry<=2*RADIUS_Y; ry++)
216  X_local[xl+rx][yl+ry] = readX(X, x-RADIUS_X+rx, y-RADIUS_Y+ry, X_width, X_height, X_offset);
217  }
218  else if (xl==WORK_GROUP_SIZE_2D-1)
219  {
220  for (int rx=0; rx<=2*RADIUS_X; rx++)
221  X_local[xl+rx][yl] = readX(X, x-RADIUS_X+rx, y-RADIUS_Y, X_width, X_height, X_offset);
222  }
223  else if (yl == WORK_GROUP_SIZE_2D-1)
224  {
225  for (int ry=0; ry<=2*RADIUS_Y; ry++)
226  X_local[xl][yl+ry] = readX(X, x-RADIUS_X, y-RADIUS_Y+ry, X_width, X_height, X_offset);
227  }
228  else
229  X_local[xl][yl] = readX(X, x-RADIUS_X, y-RADIUS_Y, X_width, X_height, X_offset);
230 
231  barrier(CLK_LOCAL_MEM_FENCE);
232 
233  if (x>=X_width || y>=X_height)
234  return;
235 
236  DATATYPE sum = 0;
237  for (int x1=0; x1<W_WIDTH; x1++)
238  {
239  #ifdef FLIP
240  int wx = x1*W_HEIGHT+W_offset;
241  #else
242  int wx = (2*RADIUS_X-x1)*W_HEIGHT+W_offset;
243  #endif
244  int inx = x1+xl;
245  for (int y1=0; y1<W_HEIGHT; y1++)
246  {
247  int iny = y1+yl;
248  #ifdef FLIP
249  sum += W[y1+wx]*X_local[inx][iny];
250  #else
251  sum += W[2*RADIUS_Y-y1+wx]*X_local[inx][iny];
252  #endif
253  }
254  }
255  #ifdef OVERWRITE
256  Y[y+X_height*x + Y_offset] = sum;
257  #else
258  Y[y+X_height*x + Y_offset] += sum;
259  #endif
260  }
261  )"
262  );
263 
264  viennacl::ocl::kernel& kernel = ocl::compile_kernel(kernel_name, source);
265 
266  kernel.local_work_size(0, OCL_WORK_GROUP_SIZE_2D);
267  kernel.local_work_size(1, OCL_WORK_GROUP_SIZE_2D);
268 
269  return kernel;
270  }
271 
273  template <class T>
274  static viennacl::ocl::kernel& generate_kernel_arbitrary_stride(
275  int32_t radius_x, int32_t radius_y, bool flip, bool overwrite)
276  {
277  std::string kernel_name =
278  "convolve_arbitrary_stride_" + ocl::get_type_string<T>() + "_" +
279  std::to_string(radius_x) + "_" + std::to_string(radius_y);
280 
281  if (flip) kernel_name.append("_flip");
282  if (overwrite) kernel_name.append("_overwrite");
283 
284  if (ocl::kernel_exists(kernel_name))
285  return ocl::get_kernel(kernel_name);
286 
287  std::string source = ocl::generate_kernel_preamble<T>(kernel_name);
288 
289  if (flip) source.append("#define FLIP\n");
290  if (overwrite) source.append("#define OVERWRITE\n");
291 
292  source.append("#define RADIUS_X " + std::to_string(radius_x) + "\n");
293  source.append("#define RADIUS_Y " + std::to_string(radius_y) + "\n");
294 
295  source.append(
296  R"(
297  #define W_WIDTH (2*RADIUS_X+1)
298  #define W_HEIGHT (2*RADIUS_Y+1)
299 
300  #define X_LOCAL_WIDTH (WORK_GROUP_SIZE_2D+2*RADIUS_X)
301  #define X_LOCAL_HEIGHT (WORK_GROUP_SIZE_2D+2*RADIUS_Y)
302 
303  __kernel void KERNEL_NAME(
304  read_only __global DATATYPE* X, int X_width, int X_height, int X_offset,
305  __constant DATATYPE* W, int W_offset,
306  __global DATATYPE* Y, int Y_offset,
307  int stride_x, int stride_y)
308  {
309  __local DATATYPE X_local[WORK_GROUP_SIZE_2D][WORK_GROUP_SIZE_2D];
310 
311  int x = get_global_id(0)*stride_x;
312  int y = get_global_id(1)*stride_y;
313 
314  int Y_width = X_width/stride_x;
315  int Y_height = X_height/stride_y;
316 
317  if (get_global_id(0)>=Y_width || get_global_id(1)>=Y_height)
318  return;
319 
320  DATATYPE sum = 0;
321  for (int x1=0; x1<W_WIDTH; x1++)
322  {
323  #ifdef FLIP
324  int wx = x1*W_HEIGHT+W_offset;
325  #else
326  int wx = (2*RADIUS_X-x1)*W_HEIGHT+W_offset;
327  #endif
328  int inx = x1+x-RADIUS_X;
329  for (int y1=0; y1<W_HEIGHT; y1++)
330  {
331  int iny = y1+y-RADIUS_Y;
332  if (inx>=0 && iny>=0 && inx<X_width && iny<X_height)
333  {
334  #ifdef FLIP
335  sum += W[y1+wx]*X[iny+inx*X_height+X_offset];
336  #else
337  sum += W[2*RADIUS_Y-y1+wx]*X[iny+inx*X_height+X_offset];
338  #endif
339  }
340  }
341  }
342  #ifdef OVERWRITE
343  Y[get_global_id(1)+Y_height*get_global_id(0) + Y_offset] = sum;
344  #else
345  Y[get_global_id(1)+Y_height*get_global_id(0) + Y_offset] += sum;
346  #endif
347  }
348  )"
349  );
350 
351  viennacl::ocl::kernel& kernel = ocl::compile_kernel(kernel_name, source);
352 
353  kernel.local_work_size(0, OCL_WORK_GROUP_SIZE_2D);
354  kernel.local_work_size(1, OCL_WORK_GROUP_SIZE_2D);
355 
356  return kernel;
357  }
358 
375  static void compute(CGPUMatrix<T> X, CGPUMatrix<T> W, CGPUMatrix<T> Y, bool flip ,
376  bool overwrite, int32_t stride_x, int32_t stride_y)
377  {
378  if (stride_x==1 && stride_y==1)
379  {
380  viennacl::ocl::kernel& kernel = generate_kernel_unity_stride<T>(
381  (W.num_cols-1)/2, (W.num_rows-1)/2, flip, overwrite);
382 
383  kernel.global_work_size(0, ocl::align_to_multiple_2d(Y.num_cols));
384  kernel.global_work_size(1, ocl::align_to_multiple_2d(Y.num_rows));
385 
386  viennacl::ocl::enqueue(kernel(
387  X.vcl_matrix(), cl_int(X.num_cols), cl_int(X.num_rows), cl_int(X.offset),
388  W.vcl_matrix(), cl_int(W.offset),
389  Y.vcl_matrix(), cl_int(Y.offset)));
390  }
391  else
392  {
393  viennacl::ocl::kernel& kernel = generate_kernel_arbitrary_stride<T>(
394  (W.num_cols-1)/2, (W.num_rows-1)/2, flip, overwrite);
395 
396  kernel.global_work_size(0, ocl::align_to_multiple_2d(Y.num_cols));
397  kernel.global_work_size(1, ocl::align_to_multiple_2d(Y.num_rows));
398 
399  viennacl::ocl::enqueue(kernel(
400  X.vcl_matrix(), cl_int(X.num_cols), cl_int(X.num_rows), cl_int(X.offset),
401  W.vcl_matrix(), cl_int(W.offset),
402  Y.vcl_matrix(), cl_int(Y.offset),
403  cl_int(stride_x), cl_int(stride_y)));
404  }
405  }
406 };
407 
408 #endif // HAVE_VIENNACL
409 
410 }
411 
412 }
413 
414 }
415 #endif // CONVOLVE_H_
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > MatrixXt
Definition: Convolve.h:93
static void compute(SGMatrix< T > X, SGMatrix< T > W, SGMatrix< T > Y, bool flip, bool overwrite, int32_t stride_x, int32_t stride_y)
Definition: Convolve.h:109
Generic class sum which provides a static compute method. This class is specialized for different typ...
Definition: Sum.h:71
Matrix::Scalar sum(Matrix m, bool no_diag=false)
Definition: Redux.h:70
index_t num_rows
Definition: SGMatrix.h:329
index_t num_cols
Definition: SGMatrix.h:331
static void compute(Matrix X, Matrix W, Matrix Y, bool flip, bool overwrite, int32_t stride_x, int32_t stride_y)
shogun matrix
Definition: Parameter.h:26
void convolve(Matrix X, Matrix W, Matrix Y, bool flip=false, bool overwrite=true, int32_t stride_x=1, int32_t stride_y=1)
Definition: Core.h:139
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18

SHOGUN Machine Learning Toolbox - Documentation