TrioCFD 1.9.9_beta
TrioCFD documentation
Loading...
Searching...
No Matches
TRUSTVect_tools.cpp
1/****************************************************************************
2* Copyright (c) 2026, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. 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.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "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.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER 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;
12* 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.
13*
14*****************************************************************************/
15
16#include <TRUSTVect.h>
17#include <TRUSTVect_tools.tpp>
18#include <TRUSTTabs.h>
19#ifndef LATATOOLS
20#include <View_Types.h>
21#include <TRUSTTrav.h>
22#endif
23#include <MD_Vector_seq.h>
24
25#ifndef LATATOOLS
26#include <Perf_counters.h>
27#endif
28
29/**************************************************************************************/
30/* Warning ! This kernels are critical for performance into several TRUST applications !
31 * Do not change implementation without using performance regression testing !
32 * You are warned.
33 */
34
35/*! @brief Determine which blocks of indices should be used to perform an operation.
36 * @param opt option specifying which vector items to process (all, sequential, real)
37 * @param md the parallel metadata vector describing the distributed structure
38 * @param vect_size_tot total size of the vector
39 * @param line_size number of components per item (line size)
40 * @param nblocs_left output: number of blocks to iterate over
41 * @return a block iterator over the selected index ranges
42 */
43template <typename _SIZE_>
44Block_Iter<_SIZE_> determine_blocks(Mp_vect_options opt, const MD_Vector& md, const _SIZE_ vect_size_tot, const int line_size, int& nblocs_left)
45{
46 nblocs_left = 1;
47#ifndef LATATOOLS
48 // Should we use the members bloc_items_to_* of the MD_Vector_* classes ?
49 // (this must be avoided when
50 // - md.valeur() is not defined (md is nul)
51 // - we want all items (VECT_SEQUENTIAL_ITEMS)
52 // - md.valeur() is a MD_Vector_seq
53 // - or md.valeur() is a MD_Vector_composite when it is an aggregation of several MD_Vector_seq)
54 const bool use_blocks = (opt != VECT_ALL_ITEMS && md && md->use_blocks());
55
56 if (use_blocks)
57 {
58 assert(opt == VECT_SEQUENTIAL_ITEMS || opt == VECT_REAL_ITEMS);
59#if INT_is_64_ == 2
60 // Should never use parallel patterns in 64b:
61 assert( (!std::is_same<_SIZE_,std::int64_t>::value) );
62#endif
63 const ArrOfInt& items_blocs = (opt == VECT_SEQUENTIAL_ITEMS) ? md->get_blocs_items_to_sum() : md->get_blocs_items_to_compute();
64 const ArrOfInt& items = (opt == VECT_SEQUENTIAL_ITEMS) ? md->get_items_to_sum() : md->get_items_to_compute();
65 assert(items_blocs.size_array() % 2 == 0);
66 nblocs_left = items_blocs.size_array() >> 1;
67 return Block_Iter<_SIZE_>(items_blocs, items); // iterator on int, but operator*() will cast and return a _SIZE_
68 }
69 else
70#endif
71 if (vect_size_tot > 0)
72 {
73 // Warning, if vect_size_tot is 0, line_size might be 0 too
74 // Compute all data, in the vector (including virtual data), build a big bloc:
75 nblocs_left = 1;
76 return Block_Iter<_SIZE_>(0, vect_size_tot / line_size); // iterator on a single (big) block
77 }
78 return Block_Iter<_SIZE_>();
79}
80
81// Explicit instanciations
82template Block_Iter<int> determine_blocks(Mp_vect_options opt, const MD_Vector& md, const int vect_size_tot, const int line_size, int& nblocs_left);
83#if INT_is_64_ == 2
84template Block_Iter<trustIdType> determine_blocks(Mp_vect_options opt, const MD_Vector& md, const trustIdType vect_size_tot, const int line_size, int& nblocs_left);
85#endif
86
87
88template<typename _TYPE_, typename _SIZE_>
89void ajoute_produit_scalaire(TRUSTVect<_TYPE_,_SIZE_>& resu, _TYPE_ alpha, const TRUSTVect<_TYPE_,_SIZE_>& vx, const TRUSTVect<_TYPE_,_SIZE_>& vy, Mp_vect_options opt)
90{
91#ifndef LATATOOLS
92 ToDo_Kokkos("critical"); // Ne semble pas utilise...
93 resu.ensureDataOnHost();
94 vx.ensureDataOnHost();
95 vy.ensureDataOnHost();
96 // Master vect donne la structure de reference, les autres vecteurs doivent avoir la meme structure.
97 const TRUSTVect<_TYPE_,_SIZE_>& master_vect = resu;
98 const int line_size = master_vect.line_size(), vect_size_tot = master_vect.size_totale();
99 const MD_Vector& md = master_vect.get_md_vector();
100 assert(vx.line_size() == line_size && vy.line_size() == line_size);
101 assert(vx.size_totale() == vect_size_tot && vy.size_totale() == vect_size_tot); // this test is necessary if md is null
102#ifndef LATATOOLS
103 assert(vx.get_md_vector() == md && vy.get_md_vector() == md);
104#endif
105 // Determine blocs of data to process, depending on " opt"
106 int nblocs_left;
107 Block_Iter<_SIZE_> bloc_itr = ::determine_blocks(opt, md, vect_size_tot, line_size, nblocs_left);
108 // Shortcut for empty arrays (avoid case line_size == 0)
109 if (bloc_itr.empty()) return;
110
111 _TYPE_ *resu_base = resu.addr();
112 const _TYPE_ *x_base = vx.addr();
113 const _TYPE_ *y_base = vy.addr();
114 for (; nblocs_left; nblocs_left--)
115 {
116 // Get index of next bloc start:
117 const int begin_bloc = (*(bloc_itr++)) * line_size, end_bloc = (*(bloc_itr++)) * line_size;
118 assert(begin_bloc >= 0 && end_bloc <= vect_size_tot && end_bloc >= begin_bloc);
119 _TYPE_ *resu_ptr = resu_base + begin_bloc;
120 const _TYPE_ *x_ptr = x_base + begin_bloc;
121 const _TYPE_ *y_ptr = y_base + begin_bloc;
122 int count = end_bloc - begin_bloc;
123 for (; count; count--)
124 {
125 const _TYPE_ x = *x_ptr;
126 const _TYPE_ y = *(y_ptr++);
127 _TYPE_& p_resu = *(resu_ptr++);
128 p_resu += alpha * x * y;
129 x_ptr++;
130 }
131 }
132 // In debug mode, put invalid values where data has not been computed
133#ifndef NDEBUG
134 invalidate_data(resu, opt);
135#endif
136 return;
137#endif // LATATOOLS
138}
139
140// Explicit instanciation for templates:
141template void ajoute_produit_scalaire<double, int>(TRUSTVect<double, int>& resu, double alpha, const TRUSTVect<double, int>& vx, const TRUSTVect<double, int>& vy, Mp_vect_options opt);
142template void ajoute_produit_scalaire<float, int>(TRUSTVect<float, int>& resu, float alpha, const TRUSTVect<float, int>& vx, const TRUSTVect<float, int>& vy, Mp_vect_options opt);
143
144
145//Process bloc function used below in operation_speciale_tres_generic
146//It is templated as a function of the in/out view location and execution spaces (Device/Host)
147#ifndef LATATOOLS
148namespace
149{
150template<typename ExecSpace, typename _TYPE_, typename _SIZE_, bool IS_MUL>
151void operation_speciale_tres_generic_kernel(TRUSTVect<_TYPE_, _SIZE_>& resu, const TRUSTVect<_TYPE_, _SIZE_>& vx, int nblocs_left,
152 Block_Iter<_SIZE_>& bloc_itr, const int line_size_vx, const _SIZE_ vect_size_tot, const int delta_line_size)
153{
154 auto vx_view= vx.template view_ro<1, ExecSpace>().data();
155 auto resu_view= resu.template view_rw<1, ExecSpace>().data();
156#ifdef TRUST_USE_GPU
157 if (nblocs_left>3) ToDo_Kokkos("nblocs_left too high, optimize by rewriting as local_operations_vect_bis_generic_kernel");
158#endif
159 for (; nblocs_left; nblocs_left--)
160 {
161 // Get index of next bloc start:
162 const int begin_bloc = (*(bloc_itr++)) * line_size_vx;
163 const int end_bloc = (*(bloc_itr++)) * line_size_vx;
164
165 assert(begin_bloc >= 0 && end_bloc <= vect_size_tot && end_bloc >= begin_bloc);
166
167 // Adjust pointers to indices
168 const int resu_start_idx = begin_bloc * delta_line_size;
169
170 Kokkos::RangePolicy<ExecSpace> policy(begin_bloc, end_bloc);
171 if (statistics().get_use_gpu()) start_gpu_timer(__KERNEL_NAME__);
172 Kokkos::parallel_for(policy, KOKKOS_LAMBDA(const int i)
173 {
174 const _TYPE_ x = vx_view[i];
175
176 //The // for could be also placed there
177 for (int j = 0; j < delta_line_size; ++j)
178 {
179 const int resu_idx = resu_start_idx + i * delta_line_size + j;
180 if (IS_MUL)
181 resu_view[resu_idx] *= x;
182 else //If it's not MUL, it's DIV
183 resu_view[resu_idx] *= ((_TYPE_)1 / x);
184 }
185 });
186 if (statistics().get_use_gpu()) end_gpu_timer(__KERNEL_NAME__, is_default_exec_space<ExecSpace>);
187 }
188}
189}
190#endif
191
192template<TYPE_OPERATION_VECT_SPEC_GENERIC _TYPE_OP_, typename _TYPE_, typename _SIZE_>
193void operation_speciale_tres_generic(TRUSTVect<_TYPE_, _SIZE_>& resu, const TRUSTVect<_TYPE_,_SIZE_>& vx, Mp_vect_options opt)
194{
195#ifndef LATATOOLS
196
197 // Check the nature of the operation
198 static constexpr bool IS_MUL = (_TYPE_OP_ == TYPE_OPERATION_VECT_SPEC_GENERIC::MUL_); //it's either MUL or DIV
199
200 // get info for computation
201 const int line_size = resu.line_size(), line_size_vx = vx.line_size(), vect_size_tot = resu.size_totale();
202 const MD_Vector& md = resu.get_md_vector();
203 // The line_size of vector resu must be a multiple of the line_size of vector vx
204 assert(line_size > 0 && line_size_vx > 0 && line_size % line_size_vx == 0);
205 const int delta_line_size = line_size / line_size_vx;
206 assert(vx.size_totale() * delta_line_size == vect_size_tot); // this test is necessary if md is null
207 assert(vx.get_md_vector() == md);
208
209 // Determine blocs of data to process, depending on " opt"
210 int nblocs_left;
211 Block_Iter<_SIZE_> bloc_itr = ::determine_blocks(opt, md, vect_size_tot, line_size, nblocs_left);
212 // Shortcut for empty arrays (avoid case line_size == 0)
213 if (bloc_itr.empty())
214 return;
215
216 bool kernelOnDevice = resu.checkDataOnDevice(vx);
217
218 //Lauch computation with the execution space and view types as (template) parameters
219 if (kernelOnDevice)
220 operation_speciale_tres_generic_kernel<Kokkos::DefaultExecutionSpace, _TYPE_, _SIZE_, IS_MUL>(resu, vx, nblocs_left, bloc_itr, line_size_vx, vect_size_tot, delta_line_size);
221 else
222 operation_speciale_tres_generic_kernel<Kokkos::DefaultHostExecutionSpace, _TYPE_, _SIZE_, IS_MUL>(resu, vx, nblocs_left, bloc_itr, line_size_vx, vect_size_tot, delta_line_size);
223
224#ifndef NDEBUG
225 // In debug mode, put invalid values where data has not been computed
226 invalidate_data(resu, opt);
227#endif
228 return;
229#else
230 Cerr << "Error! operation_speciale_tres_generic can't be called in your project!" << finl;
231 Process::exit();
232#endif
233}
234
235
236// Explicit instanciation for templates:
237template void operation_speciale_tres_generic<TYPE_OPERATION_VECT_SPEC_GENERIC::MUL_, double, int>(TRUSTVect<double, int>& resu, const TRUSTVect<double, int>& vx, Mp_vect_options opt);
238template void operation_speciale_tres_generic<TYPE_OPERATION_VECT_SPEC_GENERIC::MUL_, float, int>(TRUSTVect<float, int>& resu, const TRUSTVect<float, int>& vx, Mp_vect_options opt);
239template void operation_speciale_tres_generic<TYPE_OPERATION_VECT_SPEC_GENERIC::DIV_, double, int>(TRUSTVect<double, int>& resu, const TRUSTVect<double, int>& vx, Mp_vect_options opt);
240template void operation_speciale_tres_generic<TYPE_OPERATION_VECT_SPEC_GENERIC::DIV_, float, int>(TRUSTVect<float, int>& resu, const TRUSTVect<float, int>& vx, Mp_vect_options opt);
241
242#ifndef LATATOOLS
243namespace
244{
245template<typename ExecSpace, typename _TYPE_, typename _SIZE_, bool IS_ADD>
246void operation_speciale_generic_kernel(TRUSTVect<_TYPE_, _SIZE_>& resu, const TRUSTVect<_TYPE_, _SIZE_>& vx, _TYPE_ alpha, int nblocs_left,
247 Block_Iter<_SIZE_>& bloc_itr, const _SIZE_ vect_size_tot, const int line_size)
248{
249 auto vx_view= vx.template view_ro<1, ExecSpace>().data();
250 auto resu_view= resu.template view_rw<1, ExecSpace>().data();
251#ifdef TRUST_USE_GPU
252 if (nblocs_left>3) ToDo_Kokkos("nblocs_left too high, optimize by rewriting as local_operations_vect_bis_generic_kernel");
253#endif
254 for (; nblocs_left; nblocs_left--)
255 {
256 // Get index of next bloc start:
257 const _SIZE_ begin_bloc = (*(bloc_itr++)) * line_size;
258 const _SIZE_ end_bloc = (*(bloc_itr++)) * line_size;
259
260 assert(begin_bloc >= 0 && end_bloc <= vect_size_tot && end_bloc >= begin_bloc);
261
262 Kokkos::RangePolicy<ExecSpace> policy(begin_bloc, end_bloc);
263 if (statistics().get_use_gpu()) start_gpu_timer(__KERNEL_NAME__);
264 Kokkos::parallel_for(policy, KOKKOS_LAMBDA(const int i)
265 {
266 const _TYPE_ x = vx_view[i];
267
268 if (IS_ADD) //done at compile time
269 resu_view[i] += alpha * x;
270 else //If it's not ADD, it's SQUARE
271 resu_view[i] += alpha * x * x;
272 });
273 if (statistics().get_use_gpu()) end_gpu_timer(__KERNEL_NAME__, is_default_exec_space<ExecSpace>);
274 }
275}
276}
277#endif
278
279template <TYPE_OPERATION_VECT_SPEC _TYPE_OP_ ,typename _TYPE_, typename _SIZE_>
280void ajoute_operation_speciale_generic(TRUSTVect<_TYPE_,_SIZE_>& resu, _TYPE_ alpha,
281 const TRUSTVect<_TYPE_,_SIZE_>& vx, Mp_vect_options opt)
282{
283#ifndef LATATOOLS
284 static constexpr bool IS_ADD = (_TYPE_OP_ == TYPE_OPERATION_VECT_SPEC::ADD_);
285
286 const int line_size = resu.line_size();
287 const _SIZE_ vect_size_tot = resu.size_totale();
288 const MD_Vector& md = resu.get_md_vector();
289 assert(vx.line_size() == line_size);
290 assert(vx.size_totale() == vect_size_tot); // this test is necessary if md is null
291 // Determine blocs of data to process, depending on " opt"
292 int nblocs_left;
293 Block_Iter<_SIZE_> bloc_itr = ::determine_blocks(opt, md, vect_size_tot, line_size, nblocs_left);
294 // Shortcut for empty arrays (avoid case line_size == 0)
295 if (bloc_itr.empty()) return;
296
297 bool kernelOnDevice = resu.checkDataOnDevice(vx);
298
299 if (kernelOnDevice)
300 operation_speciale_generic_kernel<Kokkos::DefaultExecutionSpace, _TYPE_, _SIZE_, IS_ADD>(resu, vx, alpha, nblocs_left, bloc_itr, vect_size_tot, line_size);
301 else
302 operation_speciale_generic_kernel<Kokkos::DefaultHostExecutionSpace, _TYPE_, _SIZE_, IS_ADD>(resu, vx, alpha, nblocs_left, bloc_itr, vect_size_tot, line_size);
303
304
305#ifndef NDEBUG
306 invalidate_data(resu, opt);
307#endif
308 return;
309#else
310 Cerr << "Error! ajoute_operation_speciale_generic can't be called in your project!" << finl;
311 Process::exit();
312#endif
313}
314
315// Explicit instanciation for templates:
316template void ajoute_operation_speciale_generic<TYPE_OPERATION_VECT_SPEC::ADD_, double, int>(TRUSTVect<double, int>& resu, double alpha, const TRUSTVect<double, int>& vx, Mp_vect_options opt);
317template void ajoute_operation_speciale_generic<TYPE_OPERATION_VECT_SPEC::ADD_, float, int>(TRUSTVect<float, int>& resu, float alpha, const TRUSTVect<float, int>& vx, Mp_vect_options opt);
318template void ajoute_operation_speciale_generic<TYPE_OPERATION_VECT_SPEC::SQUARE_, double, int>(TRUSTVect<double, int>& resu, double alpha, const TRUSTVect<double, int>& vx, Mp_vect_options opt);
319template void ajoute_operation_speciale_generic<TYPE_OPERATION_VECT_SPEC::SQUARE_, float, int>(TRUSTVect<float, int>& resu, float alpha, const TRUSTVect<float, int>& vx, Mp_vect_options opt);
320
321#ifndef LATATOOLS
322namespace
323{
324template<typename ExecSpace, typename _TYPE_, typename _SIZE_, TYPE_OPERATOR_VECT _TYPE_OP_>
325void operator_vect_vect_generic_kernel(TRUSTVect<_TYPE_,_SIZE_>& resu, const TRUSTVect<_TYPE_, _SIZE_>& vx, int nblocs_left,
326 Block_Iter<_SIZE_>& bloc_itr, const _SIZE_ vect_size_tot, const int line_size)
327{
328 static constexpr bool IS_ADD = (_TYPE_OP_ == TYPE_OPERATOR_VECT::ADD_), IS_SUB = (_TYPE_OP_ == TYPE_OPERATOR_VECT::SUB_),
329 IS_MULT = (_TYPE_OP_ == TYPE_OPERATOR_VECT::MULT_), IS_DIV = (_TYPE_OP_ == TYPE_OPERATOR_VECT::DIV_),
330 IS_EGAL = (_TYPE_OP_ == TYPE_OPERATOR_VECT::EGAL_);
331
332#ifdef TRUST_USE_GPU
333 auto vx_view= vx.template view_ro<1, ExecSpace>().data();
334 auto resu_view= resu.template view_rw<1, ExecSpace>().data();
335#ifdef TRUST_USE_GPU
336 if (nblocs_left>3) ToDo_Kokkos("nblocs_left too high, optimize by rewriting as local_operations_vect_bis_generic_kernel");
337#endif
338 for (; nblocs_left; nblocs_left--)
339 {
340 // Get index of next bloc start:
341 const _SIZE_ begin_bloc = (*(bloc_itr++)) * line_size;
342 const _SIZE_ end_bloc = (*(bloc_itr++)) * line_size;
343
344 assert(begin_bloc >= 0 && end_bloc <= vect_size_tot && end_bloc >= begin_bloc);
345 Kokkos::RangePolicy<ExecSpace> policy(begin_bloc, end_bloc);
346 if (statistics().get_use_gpu()) start_gpu_timer(__KERNEL_NAME__);
347 Kokkos::parallel_for(policy, KOKKOS_LAMBDA(const _SIZE_ i)
348 {
349 const _TYPE_ x = vx_view[i];
350 if (IS_ADD) resu_view[i] += x;
351 if (IS_SUB) resu_view[i] -= x;
352 if (IS_MULT) resu_view[i] *= x;
353 if (IS_DIV) resu_view[i] /= x;
354 if (IS_EGAL) resu_view[i] = x;
355 });
356 if (statistics().get_use_gpu()) end_gpu_timer(__KERNEL_NAME__, is_default_exec_space<ExecSpace>);
357 }
358#else
359 // Need to keep C++ optimized (pointer) implementation for PolyMAC_CDO in Flica5
360 _TYPE_ *resu_base = resu.data();
361 const _TYPE_ *x_base = vx.data();
362 for (; nblocs_left; nblocs_left--)
363 {
364 // Get index of next bloc start:
365 const _SIZE_ begin_bloc = (*(bloc_itr++)) * line_size, end_bloc = (*(bloc_itr++)) * line_size;
366 assert(begin_bloc >= 0 && end_bloc <= vect_size_tot && end_bloc >= begin_bloc);
367 _TYPE_ *resu_ptr = resu_base + begin_bloc;
368 const _TYPE_ *x_ptr = x_base + begin_bloc;
369 for (_SIZE_ count = 0; count < end_bloc - begin_bloc ; count++)
370 {
371 const _TYPE_& x = x_ptr[count];
372 _TYPE_ &p_resu = resu_ptr[count];
373 if (IS_ADD) p_resu += x;
374 if (IS_SUB) p_resu -= x;
375 if (IS_MULT) p_resu *= x;
376 if (IS_EGAL) p_resu = x;
377 if (IS_DIV) p_resu /= x;
378 }
379 }
380#endif
381}
382}
383#endif
384
385
386template <typename _TYPE_, typename _SIZE_, TYPE_OPERATOR_VECT _TYPE_OP_>
387void operator_vect_vect_generic(TRUSTVect<_TYPE_,_SIZE_>& resu, const TRUSTVect<_TYPE_,_SIZE_>& vx, Mp_vect_options opt)
388{
389#ifndef LATATOOLS
390 const int line_size = resu.line_size();
391 const _SIZE_ vect_size_tot = resu.size_totale();
392 const MD_Vector& md = resu.get_md_vector();
393 assert(vx.line_size() == line_size);
394 assert(vx.size_totale() == vect_size_tot); // this test is necessary if md is null
395 // Determine blocs of data to process, depending on " opt"
396 int nblocs_left;
397 Block_Iter<_SIZE_> bloc_itr = ::determine_blocks(opt, md, vect_size_tot, line_size, nblocs_left);
398 // Shortcut for empty arrays (avoid case line_size == 0)
399 if (bloc_itr.empty()) return;
400
401 bool kernelOnDevice = resu.checkDataOnDevice(vx);
402
403 if (kernelOnDevice)
404 operator_vect_vect_generic_kernel<Kokkos::DefaultExecutionSpace, _TYPE_, _SIZE_, _TYPE_OP_>(resu, vx, nblocs_left, bloc_itr, vect_size_tot, line_size);
405 else
406 operator_vect_vect_generic_kernel<Kokkos::DefaultHostExecutionSpace, _TYPE_, _SIZE_, _TYPE_OP_>(resu, vx, nblocs_left, bloc_itr, vect_size_tot, line_size);
407 // In debug mode, put invalid values where data has not been computed
408#ifndef NDEBUG
409 invalidate_data(resu, opt);
410#endif
411 return;
412#else
413 Cerr << "Error! operator_vect_vect_generic can't be called in your project!" << finl;
414 Process::exit();
415#endif
416}
417// Explicit instanciation for templates:
418template void operator_vect_vect_generic<double, int, TYPE_OPERATOR_VECT::ADD_>(TRUSTVect<double, int>& resu, const TRUSTVect<double, int>& vx, Mp_vect_options opt);
419template void operator_vect_vect_generic<int, int, TYPE_OPERATOR_VECT::ADD_>(TRUSTVect<int, int>& resu, const TRUSTVect<int, int>& vx, Mp_vect_options opt);
420template void operator_vect_vect_generic<float, int, TYPE_OPERATOR_VECT::ADD_>(TRUSTVect<float, int>& resu, const TRUSTVect<float, int>& vx, Mp_vect_options opt);
421template void operator_vect_vect_generic<double, int, TYPE_OPERATOR_VECT::SUB_>(TRUSTVect<double, int>& resu, const TRUSTVect<double, int>& vx, Mp_vect_options opt);
422template void operator_vect_vect_generic<int, int, TYPE_OPERATOR_VECT::SUB_>(TRUSTVect<int, int>& resu, const TRUSTVect<int, int>& vx, Mp_vect_options opt);
423template void operator_vect_vect_generic<float, int, TYPE_OPERATOR_VECT::SUB_>(TRUSTVect<float, int>& resu, const TRUSTVect<float, int>& vx, Mp_vect_options opt);
424template void operator_vect_vect_generic<double, int, TYPE_OPERATOR_VECT::MULT_>(TRUSTVect<double, int>& resu, const TRUSTVect<double, int>& vx, Mp_vect_options opt);
425template void operator_vect_vect_generic<int, int, TYPE_OPERATOR_VECT::MULT_>(TRUSTVect<int, int>& resu, const TRUSTVect<int, int>& vx, Mp_vect_options opt);
426template void operator_vect_vect_generic<float, int, TYPE_OPERATOR_VECT::MULT_>(TRUSTVect<float, int>& resu, const TRUSTVect<float, int>& vx, Mp_vect_options opt);
427template void operator_vect_vect_generic<double, int, TYPE_OPERATOR_VECT::DIV_>(TRUSTVect<double, int>& resu, const TRUSTVect<double, int>& vx, Mp_vect_options opt);
428template void operator_vect_vect_generic<int, int, TYPE_OPERATOR_VECT::DIV_>(TRUSTVect<int, int>& resu, const TRUSTVect<int, int>& vx, Mp_vect_options opt);
429template void operator_vect_vect_generic<float, int, TYPE_OPERATOR_VECT::DIV_>(TRUSTVect<float, int>& resu, const TRUSTVect<float, int>& vx, Mp_vect_options opt);
430template void operator_vect_vect_generic<double, int, TYPE_OPERATOR_VECT::EGAL_>(TRUSTVect<double, int>& resu, const TRUSTVect<double, int>& vx, Mp_vect_options opt);
431template void operator_vect_vect_generic<int, int, TYPE_OPERATOR_VECT::EGAL_>(TRUSTVect<int, int>& resu, const TRUSTVect<int, int>& vx, Mp_vect_options opt);
432template void operator_vect_vect_generic<float, int, TYPE_OPERATOR_VECT::EGAL_>(TRUSTVect<float, int>& resu, const TRUSTVect<float, int>& vx, Mp_vect_options opt);
433
434#ifndef LATATOOLS
435namespace
436{
437template<typename ExecSpace, typename _TYPE_, typename _SIZE_, TYPE_OPERATOR_SINGLE _TYPE_OP_>
438void operator_vect_single_generic_kernel(TRUSTVect<_TYPE_,_SIZE_>& resu, const _TYPE_ x, int nblocs_left,
439 Block_Iter<_SIZE_>& bloc_itr, const _SIZE_ vect_size_tot, const int line_size)
440{
441 static constexpr bool IS_ADD = (_TYPE_OP_ == TYPE_OPERATOR_SINGLE::ADD_), IS_SUB = (_TYPE_OP_ == TYPE_OPERATOR_SINGLE::SUB_),
442 IS_MULT = (_TYPE_OP_ == TYPE_OPERATOR_SINGLE::MULT_), IS_DIV = (_TYPE_OP_ == TYPE_OPERATOR_SINGLE::DIV_), IS_EGAL = (_TYPE_OP_ == TYPE_OPERATOR_SINGLE::EGAL_),
443 IS_NEGATE = (_TYPE_OP_ == TYPE_OPERATOR_SINGLE::NEGATE_), IS_INV = (_TYPE_OP_ == TYPE_OPERATOR_SINGLE::INV_), IS_ABS = (_TYPE_OP_ == TYPE_OPERATOR_SINGLE::ABS_),
444 IS_SQRT = (_TYPE_OP_ == TYPE_OPERATOR_SINGLE::SQRT_), IS_SQUARE = (_TYPE_OP_ == TYPE_OPERATOR_SINGLE::SQUARE_);
445
446 auto resu_view= resu.template view_rw<1, ExecSpace>().data();
447#ifdef TRUST_USE_GPU
448 if (nblocs_left>3) ToDo_Kokkos("nblocs_left too high, optimize by rewriting as local_operations_vect_bis_generic_kernel");
449#endif
450 for (; nblocs_left; nblocs_left--)
451 {
452 // Get index of next bloc start:
453 const _SIZE_ begin_bloc = (*(bloc_itr++)) * line_size;
454 const _SIZE_ end_bloc = (*(bloc_itr++)) * line_size;
455
456 assert(begin_bloc >= 0 && end_bloc <= vect_size_tot && end_bloc >= begin_bloc);
457 Kokkos::RangePolicy<ExecSpace> policy(begin_bloc, end_bloc);
458 if (statistics().get_use_gpu()) start_gpu_timer(__KERNEL_NAME__);
459 Kokkos::parallel_for(policy, KOKKOS_LAMBDA(const _SIZE_ i)
460 {
461 if (IS_SUB) resu_view[i] -= x;
462 if (IS_ADD) resu_view[i] += x;
463 if (IS_MULT) resu_view[i] *= x;
464 if (IS_EGAL) resu_view[i] = x;
465 if (IS_NEGATE) resu_view[i] = -resu_view[i];
466 if (IS_ABS) resu_view[i] = (_TYPE_) Kokkos::abs(resu_view[i]);
467 if (IS_SQRT) resu_view[i] = (_TYPE_) Kokkos::sqrt(resu_view[i]);
468 if (IS_SQUARE) resu_view[i] = resu_view[i]*resu_view[i];
469 if (IS_DIV) resu_view[i] /= x;
470 if (IS_INV) resu_view[i] = (_TYPE_) ((_TYPE_)1 /resu_view[i]);
471 });
472 if (statistics().get_use_gpu()) end_gpu_timer(__KERNEL_NAME__, is_default_exec_space<ExecSpace>);
473 }
474}
475}
476#endif
477
478
479template <typename _TYPE_, typename _SIZE_, TYPE_OPERATOR_SINGLE _TYPE_OP_ >
480void operator_vect_single_generic(TRUSTVect<_TYPE_,_SIZE_>& resu, const _TYPE_ x, Mp_vect_options opt)
481{
482#ifndef LATATOOLS
483 const int line_size = resu.line_size();
484 const _SIZE_ vect_size_tot = resu.size_totale();
485 const MD_Vector& md = resu.get_md_vector();
486 // Determine blocs of data to process, depending on " opt"
487 int nblocs_left;
488 Block_Iter<_SIZE_> bloc_itr = ::determine_blocks(opt, md, vect_size_tot, line_size, nblocs_left);
489 // Shortcut for empty arrays (avoid case line_size == 0)
490 if (bloc_itr.empty()) return;
491
492 bool kernelOnDevice = resu.checkDataOnDevice();
493
494 if (kernelOnDevice)
495 operator_vect_single_generic_kernel<Kokkos::DefaultExecutionSpace, _TYPE_, _SIZE_, _TYPE_OP_>(resu, x, nblocs_left, bloc_itr, vect_size_tot, line_size);
496 else
497 operator_vect_single_generic_kernel<Kokkos::DefaultHostExecutionSpace, _TYPE_, _SIZE_, _TYPE_OP_>(resu, x, nblocs_left, bloc_itr, vect_size_tot, line_size);
498
499 // In debug mode, put invalid values where data has not been computed
500#ifndef NDEBUG
501 invalidate_data(resu, opt);
502#endif
503 return;
504#else
505 Cerr << "Error! operator_vect_single_generic can't be called in your project!" << finl;
506 Process::exit();
507#endif
508}
509// Explicit instanciation for templates:
510template void operator_vect_single_generic<double, int, TYPE_OPERATOR_SINGLE::ADD_>(TRUSTVect<double, int>& resu, const double x, Mp_vect_options opt);
511template void operator_vect_single_generic<int, int, TYPE_OPERATOR_SINGLE::ADD_>(TRUSTVect<int, int>& resu, const int x, Mp_vect_options opt);
512template void operator_vect_single_generic<float, int, TYPE_OPERATOR_SINGLE::ADD_>(TRUSTVect<float, int>& resu, const float x, Mp_vect_options opt);
513template void operator_vect_single_generic<double, int, TYPE_OPERATOR_SINGLE::SUB_>(TRUSTVect<double, int>& resu, const double x, Mp_vect_options opt);
514template void operator_vect_single_generic<int, int, TYPE_OPERATOR_SINGLE::SUB_>(TRUSTVect<int, int>& resu, const int x, Mp_vect_options opt);
515template void operator_vect_single_generic<float, int, TYPE_OPERATOR_SINGLE::SUB_>(TRUSTVect<float, int>& resu, const float x, Mp_vect_options opt);
516template void operator_vect_single_generic<double, int, TYPE_OPERATOR_SINGLE::MULT_>(TRUSTVect<double, int>& resu, const double x, Mp_vect_options opt);
517template void operator_vect_single_generic<int, int, TYPE_OPERATOR_SINGLE::MULT_>(TRUSTVect<int, int>& resu, const int x, Mp_vect_options opt);
518template void operator_vect_single_generic<float, int, TYPE_OPERATOR_SINGLE::MULT_>(TRUSTVect<float, int>& resu, const float x, Mp_vect_options opt);
519template void operator_vect_single_generic<double, int, TYPE_OPERATOR_SINGLE::DIV_>(TRUSTVect<double, int>& resu, const double x, Mp_vect_options opt);
520template void operator_vect_single_generic<int, int, TYPE_OPERATOR_SINGLE::DIV_>(TRUSTVect<int, int>& resu, const int x, Mp_vect_options opt);
521template void operator_vect_single_generic<float, int, TYPE_OPERATOR_SINGLE::DIV_>(TRUSTVect<float, int>& resu, const float x, Mp_vect_options opt);
522template void operator_vect_single_generic<double, int, TYPE_OPERATOR_SINGLE::EGAL_>(TRUSTVect<double, int>& resu, const double x, Mp_vect_options opt);
523template void operator_vect_single_generic<int, int, TYPE_OPERATOR_SINGLE::EGAL_>(TRUSTVect<int, int>& resu, const int x, Mp_vect_options opt);
524template void operator_vect_single_generic<float, int, TYPE_OPERATOR_SINGLE::EGAL_>(TRUSTVect<float, int>& resu, const float x, Mp_vect_options opt);
525template void operator_vect_single_generic<double, int, TYPE_OPERATOR_SINGLE::NEGATE_>(TRUSTVect<double, int>& resu, const double x, Mp_vect_options opt);
526template void operator_vect_single_generic<int, int, TYPE_OPERATOR_SINGLE::NEGATE_>(TRUSTVect<int, int>& resu, const int x, Mp_vect_options opt);
527template void operator_vect_single_generic<float, int, TYPE_OPERATOR_SINGLE::NEGATE_>(TRUSTVect<float, int>& resu, const float x, Mp_vect_options opt);
528template void operator_vect_single_generic<double, int, TYPE_OPERATOR_SINGLE::INV_>(TRUSTVect<double, int>& resu, const double x, Mp_vect_options opt);
529template void operator_vect_single_generic<int, int, TYPE_OPERATOR_SINGLE::INV_>(TRUSTVect<int, int>& resu, const int x, Mp_vect_options opt);
530template void operator_vect_single_generic<float, int, TYPE_OPERATOR_SINGLE::INV_>(TRUSTVect<float, int>& resu, const float x, Mp_vect_options opt);
531template void operator_vect_single_generic<double, int, TYPE_OPERATOR_SINGLE::ABS_>(TRUSTVect<double, int>& resu, const double x, Mp_vect_options opt);
532template void operator_vect_single_generic<int, int, TYPE_OPERATOR_SINGLE::ABS_>(TRUSTVect<int, int>& resu, const int x, Mp_vect_options opt);
533template void operator_vect_single_generic<float, int, TYPE_OPERATOR_SINGLE::ABS_>(TRUSTVect<float, int>& resu, const float x, Mp_vect_options opt);
534template void operator_vect_single_generic<double, int, TYPE_OPERATOR_SINGLE::SQRT_>(TRUSTVect<double, int>& resu, const double x, Mp_vect_options opt);
535template void operator_vect_single_generic<int, int, TYPE_OPERATOR_SINGLE::SQRT_>(TRUSTVect<int, int>& resu, const int x, Mp_vect_options opt);
536template void operator_vect_single_generic<float, int, TYPE_OPERATOR_SINGLE::SQRT_>(TRUSTVect<float, int>& resu, const float x, Mp_vect_options opt);
537template void operator_vect_single_generic<double, int, TYPE_OPERATOR_SINGLE::SQUARE_>(TRUSTVect<double, int>& resu, const double x, Mp_vect_options opt);
538template void operator_vect_single_generic<int, int, TYPE_OPERATOR_SINGLE::SQUARE_>(TRUSTVect<int, int>& resu, const int x, Mp_vect_options opt);
539template void operator_vect_single_generic<float, int, TYPE_OPERATOR_SINGLE::SQUARE_>(TRUSTVect<float, int>& resu, const float x, Mp_vect_options opt);
540
541#if INT_is_64_ == 2
542template void operator_vect_single_generic<trustIdType, trustIdType, TYPE_OPERATOR_SINGLE::ADD_>(TRUSTVect<trustIdType, trustIdType>& resu, const trustIdType x, Mp_vect_options opt);
543template void operator_vect_single_generic<int, trustIdType, TYPE_OPERATOR_SINGLE::ADD_>(TRUSTVect<int, trustIdType>& resu, const int x, Mp_vect_options opt);
544template void operator_vect_single_generic<float, trustIdType, TYPE_OPERATOR_SINGLE::ADD_>(TRUSTVect<float, trustIdType>& resu, const float x, Mp_vect_options opt);
545template void operator_vect_single_generic<double, trustIdType, TYPE_OPERATOR_SINGLE::ADD_>(TRUSTVect<double, trustIdType>& resu, const double x, Mp_vect_options opt);
546
547template void operator_vect_single_generic<trustIdType, trustIdType, TYPE_OPERATOR_SINGLE::SUB_>(TRUSTVect<trustIdType, trustIdType>& resu, const trustIdType x, Mp_vect_options opt);
548template void operator_vect_single_generic<int, trustIdType, TYPE_OPERATOR_SINGLE::SUB_>(TRUSTVect<int, trustIdType>& resu, const int x, Mp_vect_options opt);
549template void operator_vect_single_generic<float, trustIdType, TYPE_OPERATOR_SINGLE::SUB_>(TRUSTVect<float, trustIdType>& resu, const float x, Mp_vect_options opt);
550template void operator_vect_single_generic<double, trustIdType, TYPE_OPERATOR_SINGLE::SUB_>(TRUSTVect<double, trustIdType>& resu, const double x, Mp_vect_options opt);
551
552template void operator_vect_single_generic<double, trustIdType, TYPE_OPERATOR_SINGLE::MULT_>(TRUSTVect<double, trustIdType>& resu, const double x, Mp_vect_options opt);
553template void operator_vect_single_generic<float, trustIdType, TYPE_OPERATOR_SINGLE::MULT_>(TRUSTVect<float, trustIdType>& resu, const float x, Mp_vect_options opt);
554
555template void operator_vect_single_generic<double, trustIdType, TYPE_OPERATOR_SINGLE::DIV_>(TRUSTVect<double, trustIdType>& resu, const double x, Mp_vect_options opt);
556
557
558#endif
559
560#ifndef LATATOOLS
561namespace
562{
563template<typename ExecSpace, typename _TYPE_, typename _SIZE_,typename _TYPE_RETURN_, TYPE_OPERATION_VECT _TYPE_OP_>
564void local_extrema_vect_generic_kernel(const TRUSTVect<_TYPE_,_SIZE_>& vx, int nblocs_left, Block_Iter<_SIZE_>& bloc_itr,
565 const _SIZE_ vect_size_tot, const int line_size, _TYPE_& min_max_val, int& i_min_max)
566{
567 // Shortcut for empty arrays (avoid case line_size == 0)
568 if (bloc_itr.empty()) return ;
569
570 static constexpr bool IS_IMAX = (_TYPE_OP_ == TYPE_OPERATION_VECT::IMAX_), IS_IMIN = (_TYPE_OP_ == TYPE_OPERATION_VECT::IMIN_), IS_MAX = (_TYPE_OP_ == TYPE_OPERATION_VECT::MAX_),
571 IS_MIN = (_TYPE_OP_ == TYPE_OPERATION_VECT::MIN_), IS_MAX_ABS = (_TYPE_OP_ == TYPE_OPERATION_VECT::MAX_ABS_), IS_MIN_ABS = (_TYPE_OP_ == TYPE_OPERATION_VECT::MIN_ABS_);
572
573 //For clearer code, is it max, is it a min, is it done on absolute values ?
574 static constexpr bool IS_MAXS = (IS_MAX || IS_MAX_ABS || IS_IMAX);
575 static constexpr bool IS_MINS = (IS_MIN || IS_MIN_ABS || IS_IMIN);
576 static constexpr bool IS_ABS = (IS_MAX_ABS || IS_MIN_ABS);
577
578 // Define the reducer, based on the reduction type
579 using reducer = typename std::conditional<IS_MAXS, Kokkos::MaxLoc<_TYPE_, int>, Kokkos::MinLoc<_TYPE_, int>>::type;
580 // Define the type of what the reducer will return ( a value + a index)
581 using reducer_value_type = typename reducer::value_type;
582
583 if (not(IS_MAXS || IS_MINS)) {Process::exit("Wrong operation type in local_extrema_vect_generic_kernel");}
584
585 auto vx_view= vx.template view_ro<1, ExecSpace>().data();
586#ifdef TRUST_USE_GPU
587 if (nblocs_left>3) ToDo_Kokkos("nblocs_left too high, optimize by rewriting as local_operations_vect_bis_generic_kernel");
588#endif
589 for (; nblocs_left; nblocs_left--)
590 {
591 // Get index of next bloc start:
592 const _SIZE_ begin_bloc = (*(bloc_itr++)) * line_size;
593 const _SIZE_ end_bloc = (*(bloc_itr++)) * line_size;
594
595 //Asserts
596 assert(begin_bloc >= 0 && end_bloc <= vect_size_tot && end_bloc >= begin_bloc);
597
598 //Define Policy
599 Kokkos::RangePolicy<ExecSpace> policy(begin_bloc, end_bloc);
600
601 // Define the object in which the reduction is saved
602 reducer_value_type bloc_min_max;
603
604 //Reduction
605 if (statistics().get_use_gpu()) start_gpu_timer(__KERNEL_NAME__);
606 Kokkos::parallel_reduce(policy,
607 KOKKOS_LAMBDA(const int i, reducer_value_type& local_min_max)
608 {
609 const _TYPE_ val = (IS_ABS) ? Kokkos::abs(vx_view[i]) : vx_view[i];
610
611 if ( (IS_MAXS && val>local_min_max.val) || (IS_MINS && val<local_min_max.val) )
612 {
613 local_min_max.val=val;
614 local_min_max.loc=i; // not begin_bloc + i ? This seems to be what was done before, although this is weird to me (dont we want the global index ?)
615 }
616 }
617 ,reducer(bloc_min_max)); //Reduce in bloc_min_max
618 if (statistics().get_use_gpu()) end_gpu_timer(__KERNEL_NAME__, is_default_exec_space<ExecSpace>);
619
620 //Bloc-level reduction
621 if ( (IS_MAXS && bloc_min_max.val > min_max_val) || (IS_MINS && bloc_min_max.val < min_max_val) )
622 {
623 min_max_val=bloc_min_max.val;
624 i_min_max= bloc_min_max.loc;
625 }
626 }
627}
628}
629#endif
630
631template <typename _TYPE_, typename _SIZE_, typename _TYPE_RETURN_, TYPE_OPERATION_VECT _TYPE_OP_ >
632_TYPE_RETURN_ local_extrema_vect_generic(const TRUSTVect<_TYPE_,_SIZE_>& vx, Mp_vect_options opt)
633{
634#ifndef LATATOOLS
635
636 //Array info
637 const int line_size = vx.line_size();
638 const _SIZE_ vect_size_tot = vx.size_totale();
639 const MD_Vector& md = vx.get_md_vector();
640
641 //Asserts
642 assert(vx.line_size() == line_size);
643 assert(vx.size_totale() == vect_size_tot); // this test is necessary if md is null
644 assert(vx.get_md_vector() == md);
645
646 // Determine blocs of data to process, depending on " opt"
647 int nblocs_left;
648 Block_Iter<_SIZE_> bloc_itr = ::determine_blocks(opt, md, vect_size_tot, line_size, nblocs_left);
649
650 //Initialize results
651 _TYPE_ min_max_val = neutral_value<_TYPE_,_TYPE_OP_>(); // _TYPE_ et pas _TYPE_RETURN_ desole ...
652 int i_min_max = -1 ; // seulement pour IMAX_ et IMIN_
653
654 //Localize data
655 bool kernelOnDevice = vx.checkDataOnDevice();
656
657 //Compute reduction
658 if (kernelOnDevice)
659 local_extrema_vect_generic_kernel<Kokkos::DefaultExecutionSpace, _TYPE_, _SIZE_, _TYPE_RETURN_, _TYPE_OP_>(vx, nblocs_left, bloc_itr, vect_size_tot, line_size, min_max_val, i_min_max);
660 else
661 local_extrema_vect_generic_kernel<Kokkos::DefaultHostExecutionSpace, _TYPE_, _SIZE_, _TYPE_RETURN_, _TYPE_OP_>(vx, nblocs_left, bloc_itr, vect_size_tot, line_size, min_max_val, i_min_max);
662
663 //Return index or value
664 static constexpr bool IS_IMAX = (_TYPE_OP_ == TYPE_OPERATION_VECT::IMAX_), IS_IMIN = (_TYPE_OP_ == TYPE_OPERATION_VECT::IMIN_);
665
666 return (IS_IMAX || IS_IMIN) ? (_TYPE_RETURN_)i_min_max : (_TYPE_RETURN_)min_max_val;
667
668#else
669 Cerr << "Error! local_extrema_vect_generic can't be called in your project!" << finl;
670 Process::exit();
671 return (_TYPE_RETURN_)0; // For compil in latatools
672#endif
673}
674// Explicit instanciation for templates:
675template double local_extrema_vect_generic<double, int, double, TYPE_OPERATION_VECT::IMAX_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
676template double local_extrema_vect_generic<double, int, double, TYPE_OPERATION_VECT::IMIN_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
677template double local_extrema_vect_generic<double, int, double, TYPE_OPERATION_VECT::MAX_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
678template double local_extrema_vect_generic<double, int, double, TYPE_OPERATION_VECT::MIN_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
679template double local_extrema_vect_generic<double, int, double, TYPE_OPERATION_VECT::MAX_ABS_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
680template double local_extrema_vect_generic<double, int, double, TYPE_OPERATION_VECT::MIN_ABS_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
681template int local_extrema_vect_generic<double, int, int, TYPE_OPERATION_VECT::IMAX_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
682template int local_extrema_vect_generic<double, int, int, TYPE_OPERATION_VECT::IMIN_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
683template int local_extrema_vect_generic<double, int, int, TYPE_OPERATION_VECT::MAX_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
684template int local_extrema_vect_generic<double, int, int, TYPE_OPERATION_VECT::MIN_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
685template int local_extrema_vect_generic<double, int, int, TYPE_OPERATION_VECT::MAX_ABS_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
686template int local_extrema_vect_generic<double, int, int, TYPE_OPERATION_VECT::MIN_ABS_>(const TRUSTVect<double, int>& vx, Mp_vect_options opt);
687template int local_extrema_vect_generic<int, int, int, TYPE_OPERATION_VECT::IMAX_>(const TRUSTVect<int, int>& vx, Mp_vect_options opt);
688template int local_extrema_vect_generic<int, int, int, TYPE_OPERATION_VECT::IMIN_>(const TRUSTVect<int, int>& vx, Mp_vect_options opt);
689template int local_extrema_vect_generic<int, int, int, TYPE_OPERATION_VECT::MAX_>(const TRUSTVect<int, int>& vx, Mp_vect_options opt);
690template int local_extrema_vect_generic<int, int, int, TYPE_OPERATION_VECT::MIN_>(const TRUSTVect<int, int>& vx, Mp_vect_options opt);
691template int local_extrema_vect_generic<int, int, int, TYPE_OPERATION_VECT::MAX_ABS_>(const TRUSTVect<int, int>& vx, Mp_vect_options opt);
692template int local_extrema_vect_generic<int, int, int, TYPE_OPERATION_VECT::MIN_ABS_>(const TRUSTVect<int, int>& vx, Mp_vect_options opt);
693
694#if INT_is_64_ == 2
695template double local_extrema_vect_generic<double, trustIdType, double, TYPE_OPERATION_VECT::MAX_ABS_>(const TRUSTVect<double, trustIdType>& vx, Mp_vect_options opt);
696template int local_extrema_vect_generic<int, trustIdType, int, TYPE_OPERATION_VECT::MAX_>(const TRUSTVect<int, trustIdType>& vx, Mp_vect_options opt);
697template trustIdType local_extrema_vect_generic<trustIdType, trustIdType, trustIdType, TYPE_OPERATION_VECT::MAX_>(const TRUSTVect<trustIdType, trustIdType>& vx, Mp_vect_options opt);
698#endif
699
700#ifndef LATATOOLS
701namespace
702{
703template<typename ExecSpace, typename _TYPE_, typename _SIZE_, TYPE_OPERATION_VECT_BIS _TYPE_OP_>
704void local_operations_vect_bis_generic_kernel(const TRUSTVect<_TYPE_,_SIZE_>& vx, int nblocs_left,
705 Block_Iter<_SIZE_>& bloc_itr, const _SIZE_ vect_size_tot, const int line_size, _TYPE_& sum)
706{
707 static constexpr bool IS_SQUARE = (_TYPE_OP_ == TYPE_OPERATION_VECT_BIS::SQUARE_), IS_SUM = (_TYPE_OP_ == TYPE_OPERATION_VECT_BIS::SOMME_);
708 // Performance important point for TRUSTArray dynamic kernel to have serial mode performance:
709 // Use pointer access into Kokkos loop with [] and getting raw pointer to view with .data() !
710 auto vx_view = vx.template view_ro<1, ExecSpace>().data();
711 if (nblocs_left>3)
712 {
713 // We use flattened items_blocs cause possible huge number in parallel of nblocs_left/kernel launch (e.g. during moyenne(Ps))
714 auto items = bloc_itr.items_->template view_ro<1, ExecSpace>().data();
715 // Reduction
716 if (statistics().get_use_gpu()) start_gpu_timer(__KERNEL_NAME__);
717 Kokkos::parallel_reduce(__KERNEL_NAME__,
718 Kokkos::RangePolicy<ExecSpace>(0, bloc_itr.items_->size_array()),
719 KOKKOS_LAMBDA(const int i, _TYPE_& local_sum)
720 {
721 _SIZE_ item = items[i] * line_size;
722 const _TYPE_ x = vx_view[item];
723 if (IS_SQUARE) local_sum += x * x;
724 if (IS_SUM) local_sum += x;
725 },sum);
726 if (statistics().get_use_gpu()) end_gpu_timer(__KERNEL_NAME__, is_default_exec_space<ExecSpace>);
727 }
728 else
729 {
730 for (; nblocs_left; nblocs_left--)
731 {
732 // Get index of next bloc start:
733 const _SIZE_ begin_bloc = (*(bloc_itr++)) * line_size;
734 const _SIZE_ end_bloc = (*(bloc_itr++)) * line_size;
735 //Asserts
736 assert(begin_bloc >= 0 && end_bloc <= vect_size_tot && end_bloc >= begin_bloc);
737 //Define Policy
738 Kokkos::RangePolicy <ExecSpace> policy(begin_bloc, end_bloc);
739 // Define the bloc sum
740 _TYPE_ bloc_sum = 0;
741 //Reduction
742 if (statistics().get_use_gpu()) start_gpu_timer(__KERNEL_NAME__);
743 Kokkos::parallel_reduce(policy, KOKKOS_LAMBDA(
744 const _SIZE_ i, _TYPE_
745 &local_sum)
746 {
747 const _TYPE_ x = vx_view[i];
748 if (IS_SQUARE) local_sum += x * x;
749 if (IS_SUM) local_sum += x;
750 }
751 ,bloc_sum); //Reduce in bloc_sum
752 if (statistics().get_use_gpu()) end_gpu_timer(__KERNEL_NAME__, is_default_exec_space<ExecSpace>);
753
754 //Bloc-level reduction
755 sum += bloc_sum;
756 }
757 }
758}
759}
760#endif
761
762template <typename _TYPE_, typename _SIZE_, TYPE_OPERATION_VECT_BIS _TYPE_OP_ >
763_TYPE_ local_operations_vect_bis_generic(const TRUSTVect<_TYPE_,_SIZE_>& vx,Mp_vect_options opt)
764{
765#ifndef LATATOOLS
766 _TYPE_ sum = 0;
767 // Master vect donne la structure de reference, les autres vecteurs doivent avoir la meme structure.
768 const TRUSTVect<_TYPE_,_SIZE_>& master_vect = vx;
769 const int line_size = master_vect.line_size();
770 const _SIZE_ vect_size_tot = master_vect.size_totale();
771 const MD_Vector& md = master_vect.get_md_vector();
772 assert(vx.line_size() == line_size);
773 assert(vx.size_totale() == vect_size_tot); // this test is necessary if md is null
774 assert(vx.get_md_vector() == md);
775 // Determine blocs of data to process, depending on " opt"
776 int nblocs_left;
777 Block_Iter<_SIZE_> bloc_itr = ::determine_blocks(opt, md, vect_size_tot, line_size, nblocs_left);
778 // Shortcut for empty arrays (avoid case line_size == 0)
779 if (bloc_itr.empty()) return sum;
780
781 bool kernelOnDevice = vx.checkDataOnDevice();
782
783 if (kernelOnDevice)
784 local_operations_vect_bis_generic_kernel<Kokkos::DefaultExecutionSpace, _TYPE_, _SIZE_, _TYPE_OP_>(vx, nblocs_left, bloc_itr, vect_size_tot, line_size, sum);
785 else
786 local_operations_vect_bis_generic_kernel<Kokkos::DefaultHostExecutionSpace, _TYPE_, _SIZE_, _TYPE_OP_>(vx, nblocs_left, bloc_itr, vect_size_tot, line_size, sum);
787
788 return sum;
789#else
790 Cerr << "Error! local_operations_vect_bis_generic can't be called in your project!" << finl;
791 Process::exit();
792 return (_TYPE_)0; // For compil in latatools
793#endif
794}
795// Explicit instanciation for templates:
796template double local_operations_vect_bis_generic<double, int, TYPE_OPERATION_VECT_BIS::SQUARE_>(const TRUSTVect<double, int>& vx,Mp_vect_options opt);
797template int local_operations_vect_bis_generic<int, int, TYPE_OPERATION_VECT_BIS::SQUARE_>(const TRUSTVect<int, int>& vx,Mp_vect_options opt);
798template float local_operations_vect_bis_generic<float, int, TYPE_OPERATION_VECT_BIS::SQUARE_>(const TRUSTVect<float, int>& vx,Mp_vect_options opt);
799template double local_operations_vect_bis_generic<double, int, TYPE_OPERATION_VECT_BIS::SOMME_>(const TRUSTVect<double, int>& vx,Mp_vect_options opt);
800template int local_operations_vect_bis_generic<int, int, TYPE_OPERATION_VECT_BIS::SOMME_>(const TRUSTVect<int, int>& vx,Mp_vect_options opt);
801template float local_operations_vect_bis_generic<float, int, TYPE_OPERATION_VECT_BIS::SOMME_>(const TRUSTVect<float, int>& vx,Mp_vect_options opt);
802
803#if INT_is_64_ == 2
804template double local_operations_vect_bis_generic<double, trustIdType, TYPE_OPERATION_VECT_BIS::SOMME_>(const TRUSTVect<double, trustIdType>& vx,Mp_vect_options opt);
805#endif
806
807// ==================================================================================================================================
808// BEGIN code for debug
809#ifndef NDEBUG
810// INVALID_SCALAR is used to fill arrays when values are not computed (virtual space might not be computed by operators).
811// The value below probably triggers errors on parallel test cases but does not prevent from doing "useless" computations with it.
812#ifndef LATATOOLS
813namespace
814{
815template<typename ExecSpace, typename _TYPE_, typename _SIZE_>
816void invalidate_data_kernel(TRUSTVect<_TYPE_,_SIZE_>& resu,
817 const ArrOfInt& items_blocs, const int line_size, const int blocs_size)
818{
819 _TYPE_ invalid = (_TYPE_)-987654321;
820 auto resu_view= resu.template view_rw<1, ExecSpace>().data();
821
822 int i = 0;
823 for (int blocs_idx = 0; blocs_idx < blocs_size; blocs_idx += 2) // process data until beginning of next bloc, or end of array
824 {
825 const int bloc_end = line_size * items_blocs[blocs_idx];
826 //Define Policy
827 Kokkos::RangePolicy<ExecSpace> policy(i, bloc_end);
828 //Loop
829 if (statistics().get_use_gpu()) start_gpu_timer(__KERNEL_NAME__);
830 Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int count)
831 {
832 resu_view[count]=invalid;
833 });
834 if (statistics().get_use_gpu()) end_gpu_timer(__KERNEL_NAME__, is_default_exec_space<ExecSpace>);
835 i = items_blocs[blocs_idx+1] * line_size;
836 }
837 const _SIZE_ bloc_end = resu.size_array(); // Process until end of vector
838 //Define Policy
839 Kokkos::RangePolicy<ExecSpace> policy(i, bloc_end);
840 //Loop
841 if (statistics().get_use_gpu()) start_gpu_timer(__KERNEL_NAME__);
842 Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int count)
843 {
844 resu_view[count]=invalid;
845 });
846 if (statistics().get_use_gpu()) end_gpu_timer(__KERNEL_NAME__, is_default_exec_space<ExecSpace>);
847}
848}
849#endif
850
851template <typename _TYPE_, typename _SIZE_>
852void invalidate_data(TRUSTVect<_TYPE_,_SIZE_>& resu, Mp_vect_options opt)
853{
854#ifndef LATATOOLS
855 if (Process::is_sequential()) return; // no invalid values in sequential
856
857 const MD_Vector& md = resu.get_md_vector();
858 const int line_size = resu.line_size();
859 if (opt == VECT_ALL_ITEMS || (!md)) return; // no invalid values
860 assert(opt == VECT_SEQUENTIAL_ITEMS || opt == VECT_REAL_ITEMS);
861 const ArrOfInt& items_blocs = (opt == VECT_SEQUENTIAL_ITEMS) ? md->get_blocs_items_to_sum() : md->get_blocs_items_to_compute();
862 const int blocs_size = items_blocs.size_array();
863
864 bool kernelOnDevice = resu.checkDataOnDevice();
865
866 if (kernelOnDevice)
867 invalidate_data_kernel<Kokkos::DefaultExecutionSpace, _TYPE_, _SIZE_>(resu, items_blocs, line_size, blocs_size);
868 else
869 invalidate_data_kernel<Kokkos::DefaultHostExecutionSpace, _TYPE_, _SIZE_>(resu, items_blocs, line_size, blocs_size);
870#else
871 Cerr << "Error! invalidate_data can't be called in your project!" << finl;
872 Process::exit();
873#endif
874}
875//END code for debug
876// ==================================================================================================================================
877// Explicit instanciation for templates:
878template void invalidate_data<double>(TRUSTVect<double, int>& resu, Mp_vect_options opt);
879template void invalidate_data<float>(TRUSTVect<float, int>& resu, Mp_vect_options opt);
880template void invalidate_data<int>(TRUSTVect<int, int>& resu, Mp_vect_options opt);
881#endif /* NDEBUG */
882
883#ifndef LATATOOLS
884namespace
885{
886template<typename ExecSpace, typename _TYPE_, typename _SIZE_>
887void local_prodscal_kernel(const TRUSTVect<_TYPE_,_SIZE_>& vx, const TRUSTVect<_TYPE_,_SIZE_>& vy, int nblocs_left,
888 Block_Iter<_SIZE_>& bloc_itr, const int vect_size_tot, const int line_size, _TYPE_& sum)
889{
890 auto vx_view= vx.template view_ro<1, ExecSpace>().data();
891 auto vy_view= vy.template view_ro<1, ExecSpace>().data();
892#ifdef TRUST_USE_GPU
893 if (nblocs_left>3) ToDo_Kokkos("nblocs_left too high, optimize by rewriting as local_operations_vect_bis_generic_kernel");
894#endif
895 for (; nblocs_left; nblocs_left--)
896 {
897 // Get index of next bloc start:
898 const _SIZE_ begin_bloc = (*(bloc_itr++)) * line_size;
899 const _SIZE_ end_bloc = (*(bloc_itr++)) * line_size;
900
901 //Asserts
902 assert(begin_bloc >= 0 && end_bloc <= vect_size_tot && end_bloc >= begin_bloc);
903
904 //Define Policy
905 Kokkos::RangePolicy<ExecSpace> policy(begin_bloc, end_bloc);
906
907 // Define the bloc sum
908 _TYPE_ bloc_sum=0;
909
910 //Reduction
911 if (statistics().get_use_gpu()) start_gpu_timer(__KERNEL_NAME__);
912 Kokkos::parallel_reduce(policy, KOKKOS_LAMBDA(const _SIZE_ i, _TYPE_& local_sum)
913 {
914 local_sum += vx_view[i]*vy_view[i];
915 }
916 , Kokkos::Sum<_TYPE_>(bloc_sum)); //Reduce in bloc_sum
917
918 //timer
919 if (statistics().get_use_gpu()) end_gpu_timer(__KERNEL_NAME__, is_default_exec_space<ExecSpace>);
920
921 //Bloc-level reduction
922 sum += bloc_sum;
923 }
924}
925}
926#endif
927
928template<typename _TYPE_, typename _SIZE_>
929_TYPE_ local_prodscal(const TRUSTVect<_TYPE_,_SIZE_>& vx, const TRUSTVect<_TYPE_,_SIZE_>& vy, Mp_vect_options opt)
930{
931#ifndef LATATOOLS
932 _TYPE_ sum = 0;
933
934 const int line_size = vx.line_size();
935 const _SIZE_ vect_size_tot = vx.size_totale();
936 const MD_Vector& md = vx.get_md_vector();
937
938 assert(vx.line_size() == line_size && vy.line_size() == line_size);
939 assert(vx.size_totale() == vect_size_tot && vy.size_totale() == vect_size_tot); // this test is necessary if md is null
940 assert(vx.get_md_vector() == md && vy.get_md_vector() == md);
941 // Determine blocs of data to process, depending on " opt"
942 int nblocs_left;
943 Block_Iter<_SIZE_> bloc_itr = ::determine_blocks(opt, md, vect_size_tot, line_size, nblocs_left);
944 // Shortcut for empty arrays (avoid case line_size == 0)
945 if (bloc_itr.empty()) return sum;
946
947 bool kernelOnDevice = const_cast<TRUSTVect<_TYPE_,_SIZE_>&>(vx).checkDataOnDevice(vy);
948
949 if (kernelOnDevice)
950 local_prodscal_kernel<Kokkos::DefaultExecutionSpace, _TYPE_, _SIZE_>(vx, vy, nblocs_left, bloc_itr, vect_size_tot,line_size, sum);
951 else
952 local_prodscal_kernel<Kokkos::DefaultHostExecutionSpace, _TYPE_, _SIZE_>(vx, vy, nblocs_left, bloc_itr, vect_size_tot,line_size, sum);
953
954 return sum;
955
956#else
957 Cerr << "Error! local_prodscal can't be called in your project!" << finl;
958 Process::exit();
959 return (_TYPE_)0; // For compil in latatools
960#endif
961}
962// Explicit instanciation for templates:
963template double local_prodscal(const TRUSTVect<double, int>& vx, const TRUSTVect<double, int>& vy, Mp_vect_options opt);
964template float local_prodscal(const TRUSTVect<float, int>& vx, const TRUSTVect<float, int>& vy, Mp_vect_options opt);
MD_Vector_base::get_blocs_items_to_sum
virtual const ArrOfInt & get_blocs_items_to_sum() const =0
MD_Vector_base::get_items_to_sum
virtual const ArrOfInt & get_items_to_sum() const =0
MD_Vector_base::use_blocks
virtual bool use_blocks() const =0
MD_Vector_base::get_items_to_compute
virtual const ArrOfInt & get_items_to_compute() const =0
MD_Vector_base::get_blocs_items_to_compute
virtual const ArrOfInt & get_blocs_items_to_compute() const =0
MD_Vector
: This class is an OWN_PTR but the pointed object is shared among multiple
Definition MD_Vector.h:48
Process::exit
static void exit(int exit_code=-1)
Exit routine for TRUST within a Kokkos region.
Definition Process.cpp:466
Process::is_sequential
static bool is_sequential()
Definition Process.cpp:113
TRUSTArray::size_array
_SIZE_ size_array() const
Definition TRUSTArray.tpp:187
TRUSTArray::addr
_TYPE_ * addr()
Definition TRUSTArray.tpp:159
TRUSTArray::data
_TYPE_ * data()
Definition TRUSTArray.tpp:173
TRUSTVect::size_totale
_SIZE_ size_totale() const
Definition TRUSTVect.tpp:61
TRUSTVect::line_size
int line_size() const
Definition TRUSTVect.tpp:67
TRUSTVect::get_md_vector
virtual const MD_Vector & get_md_vector() const
Definition TRUSTVect.h:123
Block_Iter
Definition TRUSTVect_tools.tpp:31
Block_Iter::empty
bool empty() const
Definition TRUSTVect_tools.tpp:53