OmniSciDB  6686921089
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ExampleFunctions.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2021 OmniSci, Inc.
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 #ifdef HAVE_SYSTEM_TFS
18 #ifndef __CUDACC__
19 
20 #ifdef HAVE_TBB
21 #include <tbb/parallel_for.h>
22 #endif
23 #endif
24 
25 #include "ExampleFunctions.h"
26 
27 namespace Mandelbrot {
28 // Modified from original code by Akanksha Jolly,
29 // originally posted at https://www.geeksforgeeks.org/fractals-in-cc/,
30 // Author profile: https://auth.geeksforgeeks.org/user/akankshajolly/articles
31 // Used under Creative Commons License as allowed by GeeksForGeeks,
32 // see https://www.geeksforgeeks.org/copyright-information/.
33 template <typename T>
34 TEMPLATE_INLINE int32_t mandelbrot_pixel(const T cx,
35  const T cy,
36  const int32_t max_iterations) {
37  // z_real
38  T zx = 0;
39  // z_imaginary
40  T zy = 0;
41  int32_t num_iterations = 0;
42  // Calculate whether c(c_real + c_imaginary) belongs
43  // to the Mandelbrot set or not and draw a pixel
44  // at coordinates (x, y) accordingly
45  // If you reach the Maximum number of iterations
46  // and If the distance from the origin is
47  // greater terthan 2 exit the loop
48  while ((zx * zx + zy * zy < 4) && (num_iterations < max_iterations)) {
49  // Calculate Mandelbrot function
50  // z = z*z + c where z is a complex number
51  // tempx = z_real*_real - z_imaginary*z_imaginary + c_real
52  const T temp_x = zx * zx - zy * zy + cx;
53  // 2*z_real*z_imaginary + c_imaginary
54  zy = 2 * zx * zy + cy;
55  // Updating z_real = tempx
56  zx = temp_x;
57 
58  // Increment counter
59  ++num_iterations;
60  }
61  return num_iterations;
62 }
63 
64 ALWAYS_INLINE DEVICE double get_scale(const double domain_min,
65  const double domain_max,
66  const int32_t num_bins) {
67  return (domain_max - domain_min) / num_bins;
68 }
69 
70 template <typename T>
71 TEMPLATE_NOINLINE __attribute__((__used__)) void mandelbrot_impl(
72  const int32_t x_pixels,
73  const int32_t y_begin,
74  const int32_t y_end,
75  const T x_min,
76  const T y_min,
77  const T x_scale,
78  const T y_scale,
79  const int32_t max_iterations,
80  Column<T>& output_x,
81  Column<T>& output_y,
82  Column<int32_t>& output_num_iterations) {
83  // scanning every point in that rectangular area.
84  // each point represents a complex number (x + yi).
85  // iterate that complex number
86 
87  for (int32_t y = y_begin; y < y_end; ++y) {
88  // c_imaginary
89  const T cy = y * y_scale + y_min;
90  for (int32_t x = 0; x < x_pixels; ++x) {
91  // c_real
92  const T cx = x * x_scale + x_min;
93  const int32_t output_pixel = y * x_pixels + x;
94  output_x[output_pixel] = cx;
95  output_y[output_pixel] = cy;
96  output_num_iterations[output_pixel] = mandelbrot_pixel(cx, cy, max_iterations);
97  }
98  }
99 }
100 
101 #ifndef __CUDACC__
102 
103 template <typename T>
104 TEMPLATE_NOINLINE __attribute__((__used__)) int32_t mandelbrot_cpu_template(
106  const int32_t x_pixels,
107  const int32_t y_pixels,
108  const T x_min,
109  const T x_max,
110  const T y_min,
111  const T y_max,
112  const int32_t max_iterations,
113  Column<T>& output_x,
114  Column<T>& output_y,
115  Column<int32_t>& output_num_iterations) {
116  if (max_iterations < 1) {
117  return mgr.ERROR_MESSAGE("max_iterations must be a positive integer");
118  }
119  if (x_pixels < 1 || y_pixels < 1) {
120  return mgr.ERROR_MESSAGE("x_pixels and y_pixels must be positive integers");
121  }
122  if (x_max <= x_min || y_max <= y_min) {
123  return mgr.ERROR_MESSAGE(
124  "Max x and y bounds must be greater than min x and y bounds");
125  }
126 
127  const T x_scale = get_scale(x_min, x_max, x_pixels);
128  const T y_scale = get_scale(y_min, y_max, y_pixels);
129 
130  const int32_t num_pixels = x_pixels * y_pixels;
131  mgr.set_output_row_size(num_pixels);
132 #ifdef HAVE_TBB
133  tbb::parallel_for(tbb::blocked_range<int32_t>(0, y_pixels),
134  [&](const tbb::blocked_range<int32_t>& y_itr) {
135  const int32_t y_begin = y_itr.begin();
136  const int32_t y_end = y_itr.end();
137 #else
138  const int32_t y_begin = 0;
139  const int32_t y_end = y_pixels;
140 #endif
141  mandelbrot_impl(x_pixels,
142  y_begin,
143  y_end,
144  x_min,
145  y_min,
146  x_scale,
147  y_scale,
148  max_iterations,
149  output_x,
150  output_y,
151  output_num_iterations);
152 #ifdef HAVE_TBB
153  });
154 #endif
155  return num_pixels;
156 }
157 
158 #else // #ifndef __CUDACC__
159 
160 template <typename T>
161 TEMPLATE_NOINLINE __attribute__((__used__)) int32_t mandelbrot_cuda_template__gpu_(
162  const int32_t x_pixels,
163  const int32_t y_pixels,
164  const T x_min,
165  const T x_max,
166  const T y_min,
167  const T y_max,
168  const int32_t max_iterations,
169  const int32_t output_size,
170  Column<T>& output_x,
171  Column<T>& output_y,
172  Column<int32_t>& output_num_iterations) {
173  const T x_scale = get_scale(x_min, x_max, x_pixels);
174  const T y_scale = get_scale(y_min, y_max, y_pixels);
175  const int32_t num_pixels = x_pixels * y_pixels;
176 
177  int32_t start = threadIdx.x + blockDim.x * blockIdx.x;
178  int32_t step = blockDim.x * gridDim.x;
179 
180  for (int32_t output_pixel = start; output_pixel < num_pixels; output_pixel += step) {
181  const int32_t y = output_pixel / x_pixels;
182  const int32_t x = output_pixel % x_pixels;
183  const T cy = y * y_scale + y_min;
184  const T cx = x * x_scale + x_min;
185  T zx = 0;
186  T zy = 0;
187  int32_t num_iterations = 0;
188  for (; num_iterations < max_iterations; ++num_iterations) {
189  const T temp_x = zx * zx - zy * zy + cx;
190  zy = 2 * zx * zy + cy;
191  zx = temp_x;
192  if (zx * zx + zy * zy > 4.0) {
193  break;
194  }
195  }
196  output_x[output_pixel] = cx;
197  output_y[output_pixel] = cy;
198  output_num_iterations[output_pixel] = num_iterations;
199  }
200  return output_size;
201 }
202 #endif // #ifndef __CUDACC__
203 
204 } // namespace Mandelbrot
205 
206 #ifndef __CUDACC__
207 
209 __attribute__((used)) int32_t tf_mandelbrot__cpu_(
211  const int32_t x_pixels,
212  const int32_t y_pixels,
213  const double x_min,
214  const double x_max,
215  const double y_min,
216  const double y_max,
217  const int32_t max_iterations,
218  Column<double>& output_x,
219  Column<double>& output_y,
220  Column<int32_t>& output_num_iterations) {
221  return Mandelbrot::mandelbrot_cpu_template<double>(mgr,
222  x_pixels,
223  y_pixels,
224  x_min,
225  x_max,
226  y_min,
227  y_max,
228  max_iterations,
229  output_x,
230  output_y,
231  output_num_iterations);
232 }
233 
235 __attribute__((used)) int32_t tf_mandelbrot_float__cpu_(
237  const int32_t x_pixels,
238  const int32_t y_pixels,
239  const float x_min,
240  const float x_max,
241  const float y_min,
242  const float y_max,
243  const int32_t max_iterations,
244  Column<float>& output_x,
245  Column<float>& output_y,
246  Column<int32_t>& output_num_iterations) {
247  return Mandelbrot::mandelbrot_cpu_template<float>(mgr,
248  x_pixels,
249  y_pixels,
250  x_min,
251  x_max,
252  y_min,
253  y_max,
254  max_iterations,
255  output_x,
256  output_y,
257  output_num_iterations);
258 }
259 
260 #else // #ifndef __CUDACC__
261 
263 __attribute__((__used__)) int32_t tf_mandelbrot_cuda__gpu_(
264  const int32_t x_pixels,
265  const int32_t y_pixels,
266  const double x_min,
267  const double x_max,
268  const double y_min,
269  const double y_max,
270  const int32_t max_iterations,
271  const int32_t output_size,
272  Column<double>& output_x,
273  Column<double>& output_y,
274  Column<int32_t>& output_num_iterations) {
275  return Mandelbrot::mandelbrot_cuda_template__gpu_(x_pixels,
276  y_pixels,
277  x_min,
278  x_max,
279  y_min,
280  y_max,
281  max_iterations,
282  output_size,
283  output_x,
284  output_y,
285  output_num_iterations);
286 }
287 
289 __attribute__((__used__)) int32_t tf_mandelbrot_cuda_float__gpu_(
290  const int32_t x_pixels,
291  const int32_t y_pixels,
292  const float x_min,
293  const float x_max,
294  const float y_min,
295  const float y_max,
296  const int32_t max_iterations,
297  const int32_t output_size,
298  Column<float>& output_x,
299  Column<float>& output_y,
300  Column<int32_t>& output_num_iterations) {
301  return Mandelbrot::mandelbrot_cuda_template__gpu_(x_pixels,
302  y_pixels,
303  x_min,
304  x_max,
305  y_min,
306  y_max,
307  max_iterations,
308  output_size,
309  output_x,
310  output_y,
311  output_num_iterations);
312 }
313 
314 #endif // #ifndef __CUDACC__
315 #endif // HAVE_SYSTEM_TFS
#define const
#define TEMPLATE_NOINLINE
Definition: OmniSciTypes.h:30
#define DEVICE
#define EXTENSION_NOINLINE
Definition: OmniSciTypes.h:28
FORCE_INLINE T __attribute__((__may_alias__))*may_alias_ptr(T *ptr)
Definition: TypePunning.h:32
void parallel_for(const blocked_range< Int > &range, const Body &body, const Partitioner &p=Partitioner())
#define ALWAYS_INLINE
#define TEMPLATE_INLINE
Definition: OmniSciTypes.h:29