Alien  1.3.0
User documentation
Loading...
Searching...
No Matches
SYCLBEllPackInternal.h
1// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
2//-----------------------------------------------------------------------------
3// Copyright 2000-2026 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
4// See the top-level COPYRIGHT file for details.
5// SPDX-License-Identifier: Apache-2.0
6//-----------------------------------------------------------------------------
7
8
9#pragma once
10
11#include <alien/kernels/sycl/SYCLPrecomp.h>
12
13#ifdef USE_SYCL2020
14#include <sycl/sycl.hpp>
15#else
16#include <CL/sycl.hpp>
17#endif
18
19#include <alien/kernels/sycl/data/SYCLSendRecvOp.h>
20#include <alien/kernels/sycl/data/SYCLLUSendRecvOp.h>
21
22#include <alien/kernels/sycl/data/BEllPackStructInfo.h>
23
24/*---------------------------------------------------------------------------*/
25
26namespace Alien::SYCLInternal
27{
28
29#ifndef USE_SYCL2020
30 using namespace cl ;
31#endif
32
33template <int EllPackSize, typename IndexT>
34struct ALIEN_EXPORT StructInfoInternal
35{
36 // clang-format off
37 static constexpr int ellpack_size = EllPackSize ;
38 using index_type = IndexT;
39 using IndexType = IndexT;
40 using index_buffer_type = sycl::buffer<index_type, 1>;
41 using IndexBufferType = sycl::buffer<index_type, 1>;
42 using MaskBufferType = sycl::buffer<uint8_t, 1>;
43 // clang-format on
44
45 StructInfoInternal(std::size_t nrows,
46 std::size_t nnz,
47 std::size_t block_nrows,
48 std::size_t block_nnz,
49 int const* h_kcol,
50 int const* h_cols,
51 int const* h_block_row_offset,
52 int const* h_local_row_size);
53
54 IndexBufferType& getBlockRowOffset() const { return m_block_row_offset; }
55
56 IndexBufferType& getBlockCols() const { return m_block_cols; }
57
58 IndexBufferType& getKCol() const { return m_kcol; }
59
60 int const* kcol() const
61 {
62 return m_h_kcol.data();
63 }
64
65 int const* cols() const
66 {
67 return m_h_cols.data();
68 }
69
70 int const* dcol() const
71 {
72 getUpperDiagOffset();
73 return m_h_dcol.data();
74 }
75
76 void getUpperDiagOffset() const;
77 void computeLowerUpperMask() const;
78
79 MaskBufferType& getLowerMask() const;
80 MaskBufferType& getUpperMask() const;
81
82 // clang-format off
83 std::size_t m_nrows = 0 ;
84 std::size_t m_nnz = 0 ;
85 std::size_t m_block_nrows = 0 ;
86 std::size_t m_block_nnz = 0 ;
87
88 std::vector<index_type> m_h_kcol ;
89 std::vector<index_type> m_h_cols ;
90 std::vector<index_type> m_h_block_cols ;
91
92 mutable IndexBufferType m_block_row_offset ;
93 mutable IndexBufferType m_block_cols ;
94 mutable IndexBufferType m_kcol ;
95
96 mutable bool m_lower_upper_mask_ready = false ;
97 mutable std::vector<index_type> m_h_dcol ;
98 mutable std::unique_ptr<MaskBufferType> m_lower_mask ;
99 mutable std::unique_ptr<MaskBufferType> m_upper_mask ;
100 // clang-format on
101};
102
103/*---------------------------------------------------------------------------*/
104template <int N> class MultNKernel;
105template <int N> class AddUMultNKernel;
106template <int N> class AddLMultNKernel;
107template <int N> class AddExtMultNKernel;
108
109template <int Nmin, typename F, int... Is>
110bool dispatchN(int n, F&& f, std::integer_sequence<int, Is...>)
111{
112 return ((n == Nmin + Is
113 ? (f(std::integral_constant<int, Nmin + Is>{}), true)
114 : false) || ...);
115}
116
117template <int Nmin, int Nmax, typename F>
118bool dispatchN(int n, F&& f)
119{
120 static_assert(Nmax >= Nmin, "plage invalide");
121 return dispatchN<Nmin>(n, std::forward<F>(f),
122 std::make_integer_sequence<int, Nmax - Nmin + 1>{});
123}
124
125template <typename ValueT, int EllPackSize>
126class MatrixInternal
127{
128 public:
129 // clang-format off
130 using ThisType = MatrixInternal<ValueT,EllPackSize>;
131
132 static constexpr int ellpack_size = EllPackSize ;
133
134 using ValueType = ValueT;
135 using value_type = ValueT;
136
137 using ProfileType = BEllPackStructInfo<EllPackSize,int>;
138 using InternalProfileType = typename ProfileType::InternalType;
139 using IndexType = typename InternalProfileType::IndexType;
140 using IndexBufferType = typename InternalProfileType::IndexBufferType;
141 using IndexBufferPtrType = std::unique_ptr<IndexBufferType>;
142
143 using value_buffer_type = sycl::buffer<value_type, 1>;
144 using ValueBufferType = sycl::buffer<value_type, 1>;
145 using ValueBufferPtrType = std::unique_ptr<ValueBufferType>;
146
147 using QueueType = sycl::queue;
148 // clang-format onMultNKernel
149
150 template<int N>
151 struct TileT
152 {
153 static constexpr int NxN = N*N ;
154 inline std::size_t ijk(std::size_t k, int i, int j) const
155 {
156 return (k*NxN + i*N + j)*ellpack_size;
157 }
158
159 inline std::size_t ij(std::size_t local_id,int i, int j) const
160 {
161 return local_id*NxN+ i*N + j;
162 }
163 };
164
165 struct Tile
166 {
167 static const int ellpack_size = EllPackSize ;
168 int m_N = 0 ;
169 int m_NxN = 0 ;
170
171 Tile(int N)
172 : m_N(N)
173 , m_NxN(N*N)
174 {}
175
176 inline std::size_t _ijk(std::size_t k, int i, int j) const
177 {
178 return (k*m_NxN + i*m_N + j)*ellpack_size;
179 }
180
181 inline std::size_t _ij(std::size_t local_id,int i, int j) const
182 {
183 return local_id*m_NxN+ i*m_N + j;
184 }
185
186 template<typename MatrixValueAccessorT,
187 typename MatrixColAccessorT,
188 typename VectorAccessorT>
189 ValueType mult(int ieq,
190 std::size_t local_id,
191 std::size_t k,
192 MatrixColAccessorT& cols,
193 MatrixValueAccessorT& matrix,
194 VectorAccessorT& x) const
195 {
196 ValueType value = 0. ;
197 auto x_offset = cols[k*ellpack_size+local_id]*m_N ;
198 if(x_offset>=0)
199 {
200 for(int j=0;j<m_N;++j)
201 {
202 auto mat_offset = _ijk(k,ieq,j)+local_id ;
203 value += matrix[mat_offset]*x[x_offset+j] ;
204 //printf("\n %d %d %d %d : %f += %f*%f ",ieq,j,int(k),int(mat_offset),value,matrix[mat_offset],x[x_offset+j]) ;
205 }
206 }
207 return value ;
208 }
209
210 template<typename MatrixValueAccessorT,
211 typename MatrixColAccessorT,
212 typename MaskAccessorT,
213 typename VectorAccessorT>
214 ValueType mult(int ieq,
215 std::size_t local_id,
216 std::size_t k,
217 MatrixColAccessorT& cols,
218 MaskAccessorT& mask,
219 MatrixValueAccessorT& matrix,
220 VectorAccessorT& x) const
221 {
222 ValueType value = 0. ;
223 auto x_offset = cols[k*ellpack_size+local_id]*m_N ;
224 auto ma = mask[k*ellpack_size+local_id] ;
225 if(x_offset>=0 && ma==1)
226 {
227 for(int j=0;j<m_N;++j)
228 {
229 auto mat_offset = _ijk(k,ieq,j)+local_id ;
230 value += matrix[mat_offset]*x[x_offset+j] ;
231 //printf("\n %d %d %d %d : %f += %f*%f ",ieq,j,int(k),int(mat_offset),value,matrix[mat_offset],x[x_offset+j]) ;
232 }
233 }
234 return value ;
235 }
236 };
237
238 template<typename MatrixAccT,
239 typename VectorAccT,
240 typename LUAccT>
241 struct LU
242 {
243 static constexpr int ellpack_size = EllPackSize ;
244 int m_N = 0 ;
245 int m_NxN = 0 ;
246 MatrixAccT m_matrix;
247
248 LU(int N, MatrixAccT& matrix)
249 : m_N(N)
250 , m_NxN(N*N)
251 , m_matrix(matrix)
252 {}
253
254 inline std::size_t _ijk(std::size_t k, int i, int j) const
255 {
256 return (k*m_NxN + i*m_N + j)*ellpack_size;
257 }
258
259 inline std::size_t _ij(std::size_t local_id,int i, int j) const
260 {
261 return local_id*m_NxN+ i*m_N + j;
262 }
263
264 void factorize([[maybe_unused]] std::size_t global_id,
265 std::size_t local_id,
266 std::size_t block_id,
267 std::size_t kcol,
268 LUAccT m_LU) const
269 {
270 // Copy Diag Matrix in A
271 for(int i=0;i<m_N;++i)
272 for(int j=0;j<m_N;++j)
273 m_LU[_ijk(block_id,i,j)+local_id] = m_matrix[_ijk(kcol,i,j)+local_id] ;
274
275 //Factorize A = LU
276 for (int k = 0; k < m_N; ++k)
277 {
278 //assert(m_LU[_ijk(block_id,k,k)+local_id] != 0);
279 m_LU[_ijk(block_id,k,k)+local_id] = 1 / m_LU[_ijk(block_id,k,k)+local_id];
280 for (int i = k + 1; i < m_N; ++i) {
281 m_LU[_ijk(block_id,i,k)+local_id] *= m_LU[_ijk(block_id,k,k)+local_id];
282 }
283 for (int i = k + 1; i < m_N; ++i) {
284 for (int j = k + 1; j < m_N; ++j) {
285 m_LU[_ijk(block_id,i,j)+local_id] -= m_LU[_ijk(block_id,i,k)+local_id] * m_LU[_ijk(block_id,k,j)+local_id];
286 }
287 }
288 }
289 }
290
291 void inverse(std::size_t global_id,
292 std::size_t local_id,
293 std::size_t block_id,
294 LUAccT m_LU,
295 VectorAccT m_y) const
296 {
297 // SET Y to Id
298 for(int i=0;i<m_N;++i)
299 for(int j=0;j<m_N;++j)
300 m_y[_ij(global_id,i,j)] = 0. ;
301 for(int i=0;i<m_N;++i)
302 m_y[_ij(global_id,i,i)] = 1. ;
303
304 // L solve
305 for (int i = 1; i < m_N; ++i)
306 {
307 for (int j = 0; j < i; ++j)
308 {
309 for(int k=0;k<m_N;++k)
310 m_y[_ij(global_id,i,k)] -= m_LU[_ijk(block_id,i,j)+local_id] * m_y[_ij(global_id,j,k)];
311 }
312 }
313
314 // U solve
315 for (int i = m_N - 1; i >= 0; --i)
316 {
317 for (int j = m_N - 1; j > i; --j)
318 {
319 for(int k=0;k<m_N;++k)
320 m_y[_ij(global_id,i,k)] -= m_LU[_ijk(block_id,i,j)+local_id] * m_y[_ij(global_id,j,k)];
321 }
322 for(int k=0;k<m_N;++k)
323 m_y[_ij(global_id,i,k)] *= m_LU[_ijk(block_id,i,i)+local_id];
324 }
325 }
326 };
327
328 public:
329 MatrixInternal(ProfileType const* profile, int blk_size=1);
330
331 ~MatrixInternal() {}
332
333 bool setMatrixValues(ValueType const* values, bool only_host);
334 bool setMatrixValuesFromHost();
335
336 bool setMatrixValues(ValueBufferType& values);
337 bool setMatrixValues(ValueBufferType& values,
338 ValueBufferType& ext_values);
339
340 bool copy(std::size_t nb_blocks,
341 Integer block_size,
342 ValueBufferType& rhs_values,
343 Integer rhs_block_size);
344
345 bool copy(std::size_t nb_blocks,
346 Integer block_size,
347 ValueBufferType& rhs_values,
348 ValueBufferType& rhs_ext_values,
349 Integer rhs_block_size);
350
351 bool needUpdate();
352 void notifyChanges();
353 void endUpdate();
354
355 template <int N>
356 void multN(ValueBufferType& x, ValueBufferType& y, QueueType& queue) const
357 {
358 auto device = queue.get_device();
359
360 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
361 // getting the maximum work group size per thread
362 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
363 // building the best number of global thread
364
365 // clang-format off
366 std::size_t pack_size = ellpack_size;
367 auto nrows = m_profile->getNRows();
368 auto nnz = m_profile->getNnz();
369
370 auto internal_profile = m_profile->internal();
371 auto& kcol = internal_profile->getKCol();
372 auto& block_row_offset = internal_profile->getBlockRowOffset();
373 auto& block_cols = internal_profile->getBlockCols();
374
375 auto blocks_needed = (nrows + pack_size - 1) / pack_size;
376 auto blocks_target = std::max(blocks_needed, num_groups * 4UL);
377 auto total_threads = blocks_target * pack_size;
378 std::size_t wg_size = std::min<std::size_t>(pack_size, 128);
379
380 queue.submit(
381 [&](sycl::handler& cgh)
382 {
383 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
384 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
385 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
386
387 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
388 auto access_y = y.template get_access<sycl::access::mode::discard_write>(cgh);
389
390 auto tile = TileT<N>() ;
391
392 //sycl::local_accessor<ValueType, 1> lds_x{pack_size*N, cgh};
393 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{wg_size}};
394 cgh.parallel_for<MultNKernel<N>>(r,
395 [=](sycl::nd_item<1> item_id)
396 {
397 //const std::size_t local_id = item_id.get_local_id(0);
398 const std::size_t global_id = item_id.get_global_id(0);
399 const std::size_t stride = item_id.get_global_range(0);
400
401 for (auto i = global_id; i < nrows; i += stride)
402 {
403 std::size_t block_id = i / pack_size ;
404 std::size_t local_id = i % pack_size;
405
406 int begin = access_block_row_offset[block_id] ;
407 int end = access_block_row_offset[block_id+1] ;
408
409 ValueType value[N];
410 #pragma unroll
411 for (int ieq = 0; ieq < N; ++ieq)
412 value[ieq] = 0;
413
414 for (int k = begin; k < end; ++k)
415 {
416 const int col = access_cols[k * pack_size + local_id];
417 if (col >= 0)
418 {
419 ValueType xl[N];
420 #pragma unroll
421 for (int ju = 0; ju < N; ++ju)
422 xl[ju] = access_x[col * N + ju];
423
424 #pragma unroll
425 for (int ieq = 0; ieq < N; ++ieq)
426 {
427 #pragma unroll
428 for (int ju = 0; ju < N; ++ju)
429 value[ieq] += access_values[tile.ijk(k, ieq, ju) + local_id] * xl[ju];
430 }
431 }
432 }
433
434 #pragma unroll
435 for (int ieq = 0; ieq < N; ++ieq)
436 access_y[i * N + ieq] = value[ieq];
437
438 }
439 });
440 });
441 }
442 void mult(ValueBufferType& x, ValueBufferType& y) const;
443 void mult(ValueBufferType& x, ValueBufferType& y, QueueType& queue) const;
444
445 template <int N>
446 void addExtMultN(ValueType alpha,
447 ValueBufferType& x,
448 ValueBufferType& y,
449 QueueType& queue) const
450 {
451 auto device = queue.get_device();
452
453 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
454 // getting the maximum work group size per thread
455 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
456
457 const std::size_t pack_size = ellpack_size;
458 const std::size_t interface_nrow = m_ext_profile->getNRows();
459
460 auto ext_profile = m_ext_profile->internal();
461 auto& block_row_offset = ext_profile->getBlockRowOffset();
462 auto& block_cols = ext_profile->getBlockCols();
463
464 const std::size_t blocks_needed = (interface_nrow + pack_size - 1) / pack_size;
465 const std::size_t blocks_target = std::max<std::size_t>(blocks_needed, num_groups * 4);
466 const std::size_t total_threads = blocks_target * pack_size;
467 const std::size_t wg_size = std::min<std::size_t>(pack_size, 128);
468
469 queue.submit([&](sycl::handler& cgh)
470 {
471 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
472 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
473 auto access_values = m_ext_values->template get_access<sycl::access::mode::read>(cgh);
474 auto access_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
475 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
476 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
477
478 const auto tile = TileT<N>();
479
480 cgh.parallel_for<AddExtMultNKernel<N>>(
481 sycl::nd_range<1>{ sycl::range<1>{total_threads}, sycl::range<1>{wg_size} },
482 [=](sycl::nd_item<1> item)
483 {
484 const std::size_t stride = item.get_global_range(0);
485 for (std::size_t i = item.get_global_id(0); i < interface_nrow; i += stride)
486 {
487 const std::size_t block_id = i / pack_size;
488 const std::size_t local_id = i % pack_size;
489
490 const int begin = access_block_row_offset[block_id];
491 const int end = access_block_row_offset[block_id + 1];
492
493 ValueT value[N];
494 #pragma unroll
495 for (int e = 0; e < N; ++e) value[e] = 0;
496
497 for (int k = begin; k < end; ++k)
498 {
499 const int col = access_cols[k * pack_size + local_id];
500 if (col >= 0) // ignore le padding
501 {
502 ValueT xl[N];
503 #pragma unroll
504 for (int ju = 0; ju < N; ++ju) xl[ju] = access_x[col * N + ju];
505 #pragma unroll
506 for (int ieq = 0; ieq < N; ++ieq)
507 #pragma unroll
508 for (int ju = 0; ju < N; ++ju)
509 value[ieq] += access_values[tile.ijk(k, ieq, ju) + local_id] * xl[ju];
510 }
511 }
512 const std::size_t row = access_row_ids[i];
513 #pragma unroll
514 for (int e = 0; e < N; ++e) access_y[row * N + e] += alpha * value[e];
515 }
516 });
517 });
518 }
519
520 void addExtMult(ValueBufferType& x, ValueBufferType& y) const;
521 void addExtMult(ValueBufferType& x, ValueBufferType& y, QueueType& queue) const;
522
523 template<int N>
524 void addLMultN(ValueType alpha, ValueBufferType& x, ValueBufferType& y, QueueType& queue) const
525 {
526 auto device = queue.get_device();
527
528 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
529 // getting the maximum work group size per thread
530 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
531 // building the best number of global thread
532
533 std::size_t pack_size = ellpack_size;
534 auto nrows = m_profile->getNRows();
535 auto nnz = m_profile->getNnz();
536
537 auto internal_profile = m_profile->internal();
538 auto& kcol = internal_profile->getKCol();
539 auto& block_row_offset = internal_profile->getBlockRowOffset();
540 auto& block_cols = internal_profile->getBlockCols();
541
542 auto& mask = internal_profile->getLowerMask();
543 // clang-format on
544 // clang-format off
545 queue.submit(
546 [&](sycl::handler& cgh)
547 {
548 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
549 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
550 auto access_mask = mask.template get_access<sycl::access::mode::read>(cgh);
551 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
552
553
554 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
555 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
556
557 auto blocks_needed = (nrows + pack_size - 1) / pack_size;
558 auto blocks_target = std::max(blocks_needed, num_groups * 4UL);
559 auto total_threads = blocks_target * pack_size;
560 std::size_t wg_size = std::min<std::size_t>(pack_size, 128);
561
562 auto tile = TileT<N>() ;
563 //sycl::local_accessor<ValueType, 1> lds_x{pack_size*N, cgh};
564 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{wg_size}};
565 cgh.parallel_for<AddLMultNKernel<N>>(r,
566 [=](sycl::nd_item<1> item_id)
567 {
568 const std::size_t global_id = item_id.get_global_id(0);
569 const std::size_t stride = item_id.get_global_range(0);
570
571 for (auto i = global_id; i < nrows; i += stride)
572 {
573 auto block_id = i / pack_size ;
574 auto local_id = i % pack_size ;
575
576 int begin = access_block_row_offset[block_id] ;
577 int end = access_block_row_offset[block_id+1] ;
578
579 ValueType value[N];
580 #pragma unroll
581 for (int ieq = 0; ieq < N; ++ieq)
582 value[ieq] = 0;
583
584 for (int k = begin; k < end; ++k)
585 {
586 const int col = access_cols[k * pack_size + local_id];
587 if (col >= 0 && access_mask[k * pack_size + local_id])
588 {
589 ValueType xl[N];
590 #pragma unroll
591 for (int ju = 0; ju < N; ++ju)
592 xl[ju] = access_x[col * N + ju];
593
594 #pragma unroll
595 for (int ieq = 0; ieq < N; ++ieq)
596 {
597 #pragma unroll
598 for (int ju = 0; ju < N; ++ju)
599 value[ieq] += access_values[tile.ijk(k, ieq, ju) + local_id] * xl[ju];
600 }
601 }
602 }
603
604 #pragma unroll
605 for (int ieq = 0; ieq < N; ++ieq)
606 access_y[i * N + ieq] += alpha * value[ieq];
607
608
609 /* OLD
610 #pragma unroll
611 for(int ieq=0;ieq<N;++ieq)
612 {
613 ValueType value = 0. ;
614 for(int k=begin;k<end;++k)
615 {
616 //auto k = block_row_offset+j*ellpack_size+local_id ;
617 const int col = access_cols[k * pack_size + local_id];
618 if(col>=0)
619 #pragma unroll
620 for(int ju=0;ju<N;++ju)
621 lds_x[N*local_id+ju] = access_x[N*col+ju];
622 item_id.barrier(sycl::access::fence_space::local_space);
623 if(access_mask[k * pack_size + local_id])
624 #pragma unroll
625 for(int ju=0;ju<N;++ju)
626 value += access_values[tile.ijk(k,ieq,ju) + local_id] * lds_x[N*local_id+ju] ;
627 item_id.barrier(sycl::access::fence_space::local_space);
628 }
629 access_y[i*N+ieq] += alpha*value ;
630 }
631 */
632 }
633 });
634 });
635 }
636
637 template<int N>
638 void addUMultN(ValueType alpha, ValueBufferType& x, ValueBufferType& y, QueueType& queue) const
639 {
640 auto device = queue.get_device();
641
642 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
643 // getting the maximum work group size per thread
644 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
645 // building the best number of global thread
646 //auto total_threads = num_groups * ellpack_size;
647
648 // clang-format off
649 std::size_t pack_size = ellpack_size;
650 auto nrows = m_profile->getNRows() ;
651 auto nnz = m_profile->getNnz() ;
652
653 auto blocks_needed = (nrows + pack_size - 1) / pack_size;
654 auto blocks_target = std::max(blocks_needed, num_groups * 4UL);
655 auto total_threads = blocks_target * pack_size;
656 std::size_t wg_size = std::min<std::size_t>(pack_size, 128);
657
658 auto internal_profile = m_profile->internal() ;
659 auto& kcol = internal_profile->getKCol() ;
660 auto& block_row_offset = internal_profile->getBlockRowOffset() ;
661 auto& block_cols = internal_profile->getBlockCols() ;
662 auto& mask = internal_profile->getUpperMask() ;
663 // COMPUTE VALUES
664 queue.submit(
665 [&](sycl::handler& cgh)
666 {
667 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
668 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
669 auto access_mask = mask.template get_access<sycl::access::mode::read>(cgh);
670 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
671
672
673 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
674 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
675
676 auto tile = TileT<N>() ;
677 //sycl::local_accessor<ValueType, 1> lds_x{pack_size*N, cgh};
678 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{wg_size}};
679 cgh.parallel_for<AddUMultNKernel<N>>(
680 r,
681 [=](sycl::nd_item<1> item_id)
682 {
683
684 const std::size_t local_id = item_id.get_local_id(0);
685 const std::size_t global_id = item_id.get_global_id(0);
686 const std::size_t stride = item_id.get_global_range(0);
687
688 for (auto i = global_id; i < nrows; i += stride)
689 {
690 auto block_id = i / pack_size ;
691
692 auto begin = access_block_row_offset[block_id] ;
693 auto end = access_block_row_offset[block_id+1] ;
694
695 ValueType value[N];
696 #pragma unroll
697 for (int ieq = 0; ieq < N; ++ieq)
698 value[ieq] = 0;
699
700 for (int k = begin; k < end; ++k)
701 {
702 const int col = access_cols[k * pack_size + local_id];
703 if (col >= 0 && access_mask[k * pack_size + local_id])
704 {
705 ValueType xl[N];
706 #pragma unroll
707 for (int ju = 0; ju < N; ++ju)
708 xl[ju] = access_x[col * N + ju];
709
710 #pragma unroll
711 for (int ieq = 0; ieq < N; ++ieq)
712 {
713 #pragma unroll
714 for (int ju = 0; ju < N; ++ju)
715 value[ieq] += access_values[tile.ijk(k, ieq, ju) + local_id] * xl[ju];
716 }
717 }
718 }
719
720 #pragma unroll
721 for (int ieq = 0; ieq < N; ++ieq)
722 access_y[i * N + ieq] += alpha * value[ieq];
723
724
725 /*
726 #pragma unroll
727 for(int ieq=0;ieq<N;++ieq)
728 {
729 ValueType value = 0. ;
730 for(int k=begin;k<end;++k)
731 {
732 const int col = access_cols[k * pack_size + local_id];
733 if(col>=0)
734 #pragma unroll
735 for(int ju=0;ju<N;++ju)
736 lds_x[local_id*N+ju] = access_x[col*N+ju];
737 item_id.barrier(sycl::access::fence_space::local_space);
738 if(access_mask[k * pack_size + local_id])
739 #pragma unroll
740 for(int ju=0;ju<N;++ju)
741 value += access_values[tile.ijk(k,ieq,ju) + local_id] * lds_x[local_id*N+ju] ;
742 item_id.barrier(sycl::access::fence_space::local_space);
743 }
744 access_y[i*N+ieq] += alpha*value ;
745 }*/
746 }
747 });
748 });
749 }
750
751 void addLMult(ValueType alpha, ValueBufferType& x, ValueBufferType& y) const;
752 void addUMult(ValueType alpha, ValueBufferType& x, ValueBufferType& y) const;
753
754 void addLMult(ValueType alpha, ValueBufferType& x, ValueBufferType& y, QueueType& queue) const;
755 void addUMult(ValueType alpha, ValueBufferType& x, ValueBufferType& y, QueueType& queue) const;
756
757
758 template<int N>
759 void addExtLMultN(ValueType alpha, ValueBufferType& x, ValueBufferType& y, QueueType& queue) const
760 {
761 auto device = queue.get_device();
762
763 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
764 // getting the maximum work group size per thread
765 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
766 // building the best number of global thread
767 //auto total_threads = num_groups * ellpack_size;
768
769 std::size_t pack_size = ellpack_size;
770
771 auto nrows = m_profile->getNRows();
772 auto nnz = m_profile->getNnz();
773
774 auto interface_nrow = m_ext_profile->getNRows();
775 auto first_upper_ghost = m_first_upper_ghost_index ;
776 auto ext_profile = m_ext_profile->internal();
777 auto& kcol = ext_profile->getKCol();
778 auto& block_row_offset = ext_profile->getBlockRowOffset();
779 auto& block_cols = ext_profile->getBlockCols();
780
781 auto blocks_needed = (interface_nrow + pack_size - 1) / pack_size;
782 auto blocks_target = std::max(blocks_needed, num_groups * 4UL);
783 auto total_threads = blocks_target * pack_size;
784 std::size_t wg_size = std::min<std::size_t>(pack_size, 128);
785
786 {
787 // clang-format off
788
789 queue.submit([&](sycl::handler& cgh)
790 {
791 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
792 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
793 auto access_values = m_ext_values->template get_access<sycl::access::mode::read>(cgh);
794
795 auto access_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
796
797 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
798 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
799
800 auto tile = TileT<N>() ;
801 //sycl::local_accessor<ValueType, 1> lds_x{pack_size*N, cgh};
802 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{wg_size}};
803 cgh.parallel_for<class compute_ext_lmult>(
804 r,
805 [=] (sycl::nd_item<1> item_id)
806 {
807 const std::size_t global_id = item_id.get_global_id(0);
808 const std::size_t stride = item_id.get_global_range(0);
809
810 for (auto i = global_id; i < interface_nrow; i += stride)
811 {
812 auto block_id = i / pack_size ;
813 auto local_id = i % pack_size ;
814
815 auto begin = access_block_row_offset[block_id] ;
816 auto end = access_block_row_offset[block_id+1] ;
817
818
819 ValueType value[N];
820 #pragma unroll
821 for (int ieq = 0; ieq < N; ++ieq)
822 value[ieq] = 0;
823
824
825 for (int k = begin; k < end; ++k)
826 {
827 const int col = access_cols[k * pack_size + local_id];
828 if (col >= 0 && col < first_upper_ghost)
829 {
830 ValueType xl[N];
831 #pragma unroll
832 for (int ju = 0; ju < N; ++ju)
833 xl[ju] = access_x[col * N + ju];
834
835 #pragma unroll
836 for (int ieq = 0; ieq < N; ++ieq)
837 {
838 #pragma unroll
839 for (int ju = 0; ju < N; ++ju)
840 value[ieq] += access_values[tile.ijk(k, ieq, ju) + local_id] * xl[ju];
841 }
842 }
843 }
844
845 #pragma unroll
846 for (int ieq = 0; ieq < N; ++ieq)
847 access_y[access_row_ids[i] * N + ieq] += alpha * value[ieq];
848
849 /*
850 int block_row_offset = access_block_row_offset[block_id] * pack_size ;
851 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
852 for(int ieq=0;ieq<N;++ieq)
853 {
854 ValueType value = 0. ;
855 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
856 {
857 // value += access_values[k]* access_x[access_cols[k]] ;
858 if(access_cols[k]<first_upper_ghost)
859 value += tile.mult(ieq,
860 local_id,
861 k,
862 access_cols,
863 access_values,
864 access_x) ;
865 }
866 access_y[access_row_ids[i]*N+ieq] += alpha * value;
867 }*/
868 }
869 });
870 });
871 // clang-format on
872 }
873 }
874 void addExtLMult(ValueType alpha, ValueBufferType& x, ValueBufferType& y) const;
875 void addExtLMult(ValueType alpha, ValueBufferType& x, ValueBufferType& y, QueueType& queue) const;
876
877 template<int N>
878 void addExtUMultN(ValueType alpha, ValueBufferType& x, ValueBufferType& y, QueueType& queue) const
879 {
880 auto device = queue.get_device();
881
882 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
883 // getting the maximum work group size per thread
884 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
885 // building the best number of global thread
886 //auto total_threads = num_groups * ellpack_size;
887
888 std::size_t pack_size = ellpack_size;
889
890
891 auto nrows = m_profile->getNRows();
892 auto nnz = m_profile->getNnz();
893
894 auto interface_nrow = m_ext_profile->getNRows();
895 auto first_upper_ghost = m_first_upper_ghost_index ;
896 auto ext_profile = m_ext_profile->internal();
897 auto& kcol = ext_profile->getKCol();
898 auto& block_row_offset = ext_profile->getBlockRowOffset();
899 auto& block_cols = ext_profile->getBlockCols();
900
901 auto blocks_needed = (interface_nrow + ellpack_size - 1) / ellpack_size;
902 auto blocks_target = std::max(blocks_needed, num_groups * 4UL);
903 auto total_threads = blocks_target * ellpack_size;
904 {
905 // clang-format off
906
907 queue.submit([&](sycl::handler& cgh)
908 {
909 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
910 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
911 auto access_values = m_ext_values->template get_access<sycl::access::mode::read>(cgh);
912
913 auto access_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
914
915 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
916 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
917
918 auto tile = TileT<N>() ;
919 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{pack_size}};
920 cgh.parallel_for<class compute_ext_lmult>(
921 r,
922 [=] (sycl::nd_item<1> item_id)
923 {
924 const std::size_t global_id = item_id.get_global_id(0);
925 const std::size_t stride = item_id.get_global_range(0);
926
927 for (auto i = global_id; i < interface_nrow; i += stride)
928 {
929 auto block_id = i / pack_size ;
930 auto local_id = i % pack_size ;
931
932
933 auto begin = access_block_row_offset[block_id] ;
934 auto end = access_block_row_offset[block_id+1] ;
935
936
937 ValueType value[N];
938 #pragma unroll
939 for (int ieq = 0; ieq < N; ++ieq)
940 value[ieq] = 0;
941
942
943 for (int k = begin; k < end; ++k)
944 {
945 const int col = access_cols[k * pack_size + local_id];
946 if (col >= 0 && col >= first_upper_ghost)
947 {
948 ValueType xl[N];
949 #pragma unroll
950 for (int ju = 0; ju < N; ++ju)
951 xl[ju] = access_x[col * N + ju];
952
953 #pragma unroll
954 for (int ieq = 0; ieq < N; ++ieq)
955 {
956 #pragma unroll
957 for (int ju = 0; ju < N; ++ju)
958 value[ieq] += access_values[tile.ijk(k, ieq, ju) + local_id] * xl[ju];
959 }
960 }
961 }
962
963 #pragma unroll
964 for (int ieq = 0; ieq < N; ++ieq)
965 access_y[access_row_ids[i] * N + ieq] += alpha * value[ieq];
966
967 /*
968
969 int block_row_offset = access_block_row_offset[block_id] * pack_size ;
970 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
971 for(int ieq=0;ieq<N;++ieq)
972 {
973 ValueType value = 0. ;
974 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
975 {
976 // value += access_values[k]* access_x[access_cols[k]] ;
977 if(access_cols[k]<first_upper_ghost)
978 value += tile.mult(ieq,
979 local_id,
980 k,
981 access_cols,
982 access_values,
983 access_x) ;
984 }
985 access_y[access_row_ids[i]*N+ieq] += alpha * value;
986 }*/
987 }
988 });
989 });
990 // clang-format on
991 }
992
993 }
994 void addExtUMult(ValueType alpha, ValueBufferType& x, ValueBufferType& y) const;
995 void addExtUMult(ValueType alpha, ValueBufferType& x, ValueBufferType& y, QueueType& queue) const;
996
997 void multDiag(ValueBufferType& x, ValueBufferType& y) const;
998 void multDiag(ValueBufferType& x, ValueBufferType& y, QueueType& queue) const;
999
1000 void multDiag(ValueBufferType& y) const;
1001 void multDiag(ValueBufferType& y, QueueType& queue) const;
1002
1003 void computeDiag(ValueBufferType& y) const;
1004 void computeDiag(ValueBufferType& y, QueueType& queue) const;
1005
1006 void computeBlockDiag(ValueBufferType& y) const;
1007 void computeBlockDiag(ValueBufferType& y, QueueType& queue) const;
1008
1009 void multInvDiag(ValueBufferType& y) const;
1010 void multInvDiag(ValueBufferType& y, QueueType& queue) const;
1011
1012 void computeInvDiag(ValueBufferType& y) const;
1013 void computeInvDiag(ValueBufferType& y, QueueType& queue) const;
1014
1015 void computeInvBlockDiag(ValueBufferType& y) const;
1016 void computeInvBlockDiag(ValueBufferType& y, QueueType& queue) const;
1017
1018 void scal(ValueBufferType& y);
1019
1020 void scal(ValueBufferType& y, QueueType& queue);
1021
1022 void copyDevicePointers(int local_offset,
1023 std::size_t nrows,
1024 std::size_t nnz,
1025 int* rows,
1026 int* ncols,
1027 int* cols,
1028 ValueT* values) const ;
1029
1030 ValueBufferType& getValues() { return m_values; }
1031
1032 ValueBufferType const getValues() const { return m_values; }
1033
1034 //ProfileType* getProfile() { return m_profile; }
1035
1036 ProfileType const* getProfile() const { return m_profile; }
1037
1038 ValueType const* getHCsrData() const
1039 {
1040 return m_h_csr_values.data();
1041 }
1042
1043 ValueType* getHCsrData()
1044 {
1045 return m_h_csr_values.data();
1046 }
1047
1048 IndexBufferType& getSendIds() const
1049 {
1050 return *m_send_ids;
1051 }
1052 IndexBufferType& getRecvIds() const
1053 {
1054 return *m_recv_ids;
1055 }
1056
1057 // clang-format off
1058 int m_N = 1;
1059 int m_NxN = 1;
1060 ProfileType const* m_profile = nullptr;
1061 ProfileType const* m_ext_profile = nullptr;
1062 int m_first_upper_ghost_index = 0 ;
1063
1064 std::vector<ValueType> m_h_csr_values ;
1065 std::vector<ValueType> m_h_values ;
1066 mutable ValueBufferType m_values ;
1067
1068 std::vector<ValueType> m_h_csr_ext_values ;
1069 std::vector<ValueType> m_h_ext_values ;
1070 mutable ValueBufferPtrType m_ext_values ;
1071 bool m_values_is_update = false ;
1072
1073 int const* m_h_interface_row_ids = nullptr;
1074 mutable IndexBufferPtrType m_interface_row_ids ;
1075 mutable IndexBufferPtrType m_send_ids ;
1076 mutable IndexBufferPtrType m_recv_ids ;
1077 mutable IndexBufferPtrType m_recv_uids ;
1078 // clang-format on
1079};
1080
1081
1082/*---------------------------------------------------------------------------*/
1083
1084} // namespace Alien::SYCLInternal
1085
1086/*---------------------------------------------------------------------------*/