Alien  1.3.0
Developer documentation
Loading...
Searching...
No Matches
SYCLBEllPackMatrix.cc
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 * Copyright 2020 IFPEN-CEA
9 *
10 * Licensed under the Apache License, Version 2.0 (the "License");
11 * you may not use this file except in compliance with the License.
12 * You may obtain a copy of the License at
13 *
14 * http://www.apache.org/licenses/LICENSE-2.0
15 *
16 * Unless required by applicable law or agreed to in writing, software
17 * distributed under the License is distributed on an "AS IS" BASIS,
18 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
19 * See the License for the specific language governing permissions and
20 * limitations under the License.
21 *
22 * SPDX-License-Identifier: Apache-2.0
23 */
24
25#include <cassert>
26
27#include <alien/kernels/sycl/data/SYCLEnvInternal.h>
28#include <alien/kernels/sycl/data/SYCLEnv.h>
29
30#include <alien/kernels/sycl/data/SYCLVectorInternal.h>
31#include <alien/kernels/sycl/data/SYCLVector.h>
32
33#include <alien/kernels/sycl/data/SYCLBEllPackInternal.h>
34#include <alien/kernels/sycl/data/SYCLBEllPackMatrix.h>
35
36/*---------------------------------------------------------------------------*/
37
38namespace Alien
39{
40
41SYCLEnv::SYCLEnv()
42{
43 SYCLInternal::EnvInternal::printPlatformInfo() ;
44 m_internal.reset(new SYCLInternal::EnvInternal{});
45}
46
47SYCLEnv::SYCLEnv(int device_id)
48{
49 SYCLInternal::EnvInternal::printPlatformInfo() ;
50 auto gpu_devices = sycl::device::get_devices(sycl::info::device_type::gpu);
51 if(device_id < gpu_devices.size())
52 {
53 sycl::device selected_device = gpu_devices[device_id];
54 m_internal.reset(new SYCLInternal::EnvInternal{selected_device,device_id});
55 }
56 else
57 m_internal.reset(new SYCLInternal::EnvInternal{});
58}
59
60
61SYCLEnv::~SYCLEnv()
62{
63}
64
65SYCLEnv* SYCLEnv::m_instance = nullptr;
66
67SYCLEnv* SYCLEnv::instance()
68{
69 if (!m_instance)
70 m_instance = new SYCLEnv{};
71 return m_instance;
72}
73
74SYCLEnv* SYCLEnv::instance(int device_id)
75{
76 if (!m_instance)
77 m_instance = new SYCLEnv{device_id};
78 if (m_instance->deviceId()!= device_id)
79 {
80 delete m_instance ;
81 m_instance = new SYCLEnv{device_id};
82 }
83 return m_instance;
84}
85
86int SYCLEnv::deviceId()
87{
88 return m_internal->deviceId() ;
89}
90
91std::size_t SYCLEnv::maxNumGroups()
92{
93 return m_internal->maxNumGroups();
94}
95
96std::size_t SYCLEnv::maxWorkGroupSize()
97{
98 return m_internal->maxWorkGroupSize();
99}
100
101std::size_t SYCLEnv::maxNumThreads()
102{
103 return m_internal->maxNumThreads();
104}
105
106template <int EllpackSize, typename IndexT>
107void BEllPackStructInfo<EllpackSize, IndexT>::computeBlockRowOffset(std::vector<int>& block_row_offset,
108 std::size_t nrows,
109 int const* kcol)
110{
111 std::size_t block_nrows = BEllPackStructInfo<EllpackSize, IndexT>::nbBlocks(nrows);
112 block_row_offset.resize(block_nrows + 1);
113 int offset = 0;
114 for (std::size_t ib = 0; ib < block_nrows; ++ib) {
115 block_row_offset[ib] = offset;
116 int max_row_size = 0;
117 for (int i = 0; i < std::min(EllpackSize, int(nrows - ib * EllpackSize)); ++i) {
118 auto irow = ib * EllpackSize + i;
119 int row_size = kcol[irow + 1] - kcol[irow];
120 max_row_size = std::max(max_row_size, row_size);
121 }
122 offset += max_row_size;
123 }
124 block_row_offset[block_nrows] = offset;
125}
126
127template <int EllpackSize, typename IndexT>
128SYCLInternal::StructInfoInternal<EllpackSize, IndexT>::
129StructInfoInternal(std::size_t nrows,
130 std::size_t nnz,
131 std::size_t block_nrows,
132 std::size_t block_nnz,
133 int const* h_kcol,
134 int const* h_cols,
135 int const* h_block_row_offset,
136 int const* h_local_row_size)
137: m_nrows(nrows)
138, m_nnz(nnz)
139, m_block_nrows(block_nrows)
140, m_block_nnz(block_nnz)
141, m_h_kcol(h_kcol, h_kcol + nrows + 1)
142, m_h_cols(h_cols, h_cols + nnz)
143, m_h_block_cols(block_nnz * ellpack_size)
144, m_block_row_offset(h_block_row_offset, sycl::range<1>(m_block_nrows + 1))
145, m_block_cols(m_h_block_cols.data(), sycl::range<1>(m_block_nnz * ellpack_size))
146, m_kcol(m_h_kcol.data(), sycl::range<1>(m_nrows + 1))
147{
148 auto env = SYCLEnv::instance();
149
150 auto& queue = env->internal()->queue();
151
152 auto num_groups = env->internal()->maxNumGroups();
153
154 auto max_work_group_size = env->internal()->maxWorkGroupSize();
155
156 // building the best number of global thread
157 auto total_threads = num_groups * ellpack_size;
158
159 // clang-format off
160 if(h_local_row_size==nullptr)
161 {
162 IndexBufferType cols_buffer(m_h_cols.data(), sycl::range<1>(m_nnz));
163
164 queue.submit([&](sycl::handler& cgh)
165 {
166 auto access_kcol_buffer = m_kcol.template get_access<sycl::access::mode::read>(cgh);
167 auto access_block_row_offset = m_block_row_offset.template get_access<sycl::access::mode::read>(cgh);
168
169 auto access_cols_buffer = cols_buffer.template get_access<sycl::access::mode::read>(cgh);
170 auto access_block_cols = m_block_cols.template get_access<sycl::access::mode::discard_write>(cgh);
171
172 cgh.parallel_for<class struct_info_sycl>(sycl::range<1>{total_threads},
173 [=] (sycl::item<1> itemId)
174 {
175 auto id = itemId.get_id(0);
176
177 for (auto i = id; i < nrows; i += itemId.get_range()[0])
178 {
179 auto block_id = i/ellpack_size ;
180 auto local_id = i%ellpack_size ;
181
182 int begin = access_kcol_buffer[i] ;
183 int end = access_kcol_buffer[i+1] ;
184 int row_size = end - begin ;
185
186 int block_row_offset = access_block_row_offset[block_id]*ellpack_size ;
187 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
188
189 //access_rx[i] = row_size;
190 /*
191 access_brs[i] = block_row_size ;
192 access_begin[i] = begin ;
193 access_end[i] = end ;
194 access_first_cols[i] = access_cols_buffer[begin] ;
195 */
196 for(int k=begin;k<end;++k)
197 {
198 access_block_cols[block_row_offset+(k-begin)*ellpack_size+local_id] = access_cols_buffer[k] ;
199 }
200 for(int k=row_size;k<block_row_size;++k)
201 {
202 access_block_cols[block_row_offset+k*ellpack_size+local_id] = -1 ;
203 }
204
205 }
206 });
207 });
208 }
209 else
210 {
211 IndexBufferType cols_buffer(m_h_cols.data(), sycl::range<1>(m_nnz));
212 IndexBufferType lrowsize_buffer(h_local_row_size, sycl::range<1>(nrows));
213 // clang-format off
214 queue.submit([&](sycl::handler& cgh)
215 {
216 auto access_kcol_buffer = m_kcol.template get_access<sycl::access::mode::read>(cgh);
217 auto access_block_row_offset = m_block_row_offset.template get_access<sycl::access::mode::read>(cgh);
218 auto access_lrowsize_buffer = lrowsize_buffer.template get_access<sycl::access::mode::read>(cgh);
219
220 auto access_cols_buffer = cols_buffer.template get_access<sycl::access::mode::read>(cgh);
221 auto access_block_cols = m_block_cols.template get_access<sycl::access::mode::discard_write>(cgh);
222
223 cgh.parallel_for<class struct_info_sycl2>(sycl::range<1>{total_threads},
224 [=] (sycl::item<1> itemId)
225 {
226 auto id = itemId.get_id(0);
227
228 for (auto i = id; i < nrows; i += itemId.get_range()[0])
229 {
230 auto block_id = i/ellpack_size ;
231 auto local_id = i%ellpack_size ;
232
233 int begin = access_kcol_buffer[i] ;
234 int lrow_size = access_lrowsize_buffer[i] ;
235
236 int block_row_offset = access_block_row_offset[block_id]*ellpack_size ;
237 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
238
239 for(int k=begin;k<begin+lrow_size;++k)
240 {
241 access_block_cols[block_row_offset+(k-begin)*ellpack_size+local_id] = access_cols_buffer[k] ;
242 }
243 for(int k=lrow_size;k<block_row_size;++k)
244 {
245 access_block_cols[block_row_offset+k*ellpack_size+local_id] = -1 ;
246 }
247 }
248 });
249 });
250 }
251 // clang-format on
252}
253
254template <int EllpackSize, typename IndexT>
255void SYCLInternal::StructInfoInternal<EllpackSize, IndexT>::getUpperDiagOffset() const
256{
257 if (m_h_dcol.size() == 0) {
258 m_h_dcol.resize(m_nrows);
259 for (std::size_t irow = 0; irow < m_nrows; ++irow) {
260 int index = m_h_kcol[irow];
261 for (int k = m_h_kcol[irow]; k < m_h_kcol[irow + 1]; ++k) {
262 if ((std::size_t)m_h_cols[k] < irow)
263 ++index;
264 else
265 break;
266 }
267 m_h_dcol[irow] = index;
268 }
269 }
270}
271
272template <int EllpackSize, typename IndexT>
273void SYCLInternal::StructInfoInternal<EllpackSize, IndexT>::computeLowerUpperMask() const
274{
275
276 if (not m_lower_upper_mask_ready) {
277 m_lower_mask.reset(new MaskBufferType(sycl::range<1>(m_block_nnz * ellpack_size)));
278 m_upper_mask.reset(new MaskBufferType(sycl::range<1>(m_block_nnz * ellpack_size)));
279
280 auto env = SYCLEnv::instance();
281
282 auto& queue = env->internal()->queue();
283
284 auto num_groups = env->internal()->maxNumGroups();
285
286 auto max_work_group_size = env->internal()->maxWorkGroupSize();
287
288 // building the best number of global thread
289 auto total_threads = num_groups * ellpack_size;
290
291 {
292 IndexBufferType dcol_buffer(m_h_dcol.data(), sycl::range<1>(m_nrows));
293
294 auto nrows = m_nrows;
295
296 // clang-format off
297 queue.submit([&](sycl::handler& cgh)
298 {
299 auto access_dcol_buffer = dcol_buffer.template get_access<sycl::access::mode::read>(cgh);
300
301 auto access_kcol_buffer = m_kcol.template get_access<sycl::access::mode::read>(cgh);
302 auto access_block_row_offset = m_block_row_offset.template get_access<sycl::access::mode::read>(cgh);
303
304 auto access_lower_mask = sycl::accessor { *m_lower_mask, cgh, sycl::write_only, sycl::property::no_init{}};
305 auto access_upper_mask = sycl::accessor { *m_upper_mask, cgh, sycl::write_only, sycl::property::no_init{}};
306
307
308 cgh.parallel_for<class lower_upper_mask>(sycl::range<1>{total_threads},
309 [=] (sycl::item<1> itemId)
310 {
311 auto id = itemId.get_id(0);
312
313 for (auto i = id; i < nrows; i += itemId.get_range()[0])
314 {
315 auto block_id = i/ellpack_size ;
316 auto local_id = i%ellpack_size ;
317
318 int begin = access_kcol_buffer[i] ;
319 int end = access_kcol_buffer[i+1] ;
320 int diag = access_dcol_buffer[i] ;
321 int row_size = end - begin ;
322
323 int block_row_offset = access_block_row_offset[block_id]*ellpack_size ;
324 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
325
326 access_lower_mask[block_row_offset+(diag-begin)*ellpack_size+local_id] = 0 ;
327 access_upper_mask[block_row_offset+(diag-begin)*ellpack_size+local_id] = 0 ;
328
329 for(int k=begin;k<diag;++k)
330 {
331 access_lower_mask[block_row_offset+(k-begin)*ellpack_size+local_id] = 1 ;
332 access_upper_mask[block_row_offset+(k-begin)*ellpack_size+local_id] = 0 ;
333 }
334 for(int k=diag+1;k<end;++k)
335 {
336 access_lower_mask[block_row_offset+(k-begin)*ellpack_size+local_id] = 0 ;
337 access_upper_mask[block_row_offset+(k-begin)*ellpack_size+local_id] = 1 ;
338 }
339 for(int k=row_size;k<block_row_size;++k)
340 {
341 access_lower_mask[block_row_offset+k*ellpack_size+local_id] = 0 ;
342 access_upper_mask[block_row_offset+k*ellpack_size+local_id] = 0 ;
343 }
344
345 }
346 });
347 });
348 // clang-format on
349 }
350 m_lower_upper_mask_ready = true;
351 }
352}
353
354template <int EllpackSize, typename IndexT>
355typename SYCLInternal::StructInfoInternal<EllpackSize, IndexT>::MaskBufferType&
356SYCLInternal::StructInfoInternal<EllpackSize, IndexT>::getLowerMask() const
357{
358 computeLowerUpperMask();
359 return *m_lower_mask;
360}
361
362template <int EllpackSize, typename IndexT>
363typename SYCLInternal::StructInfoInternal<EllpackSize, IndexT>::MaskBufferType&
364SYCLInternal::StructInfoInternal<EllpackSize, IndexT>::getUpperMask() const
365{
366 computeLowerUpperMask();
367 return *m_upper_mask;
368}
369
370template <int EllpackSize, typename IndexT>
371BEllPackStructInfo<EllpackSize, IndexT>::BEllPackStructInfo(std::size_t nrows,
372 int const* kcol,
373 int const* cols,
374 int const* h_block_row_offset,
375 int const* h_local_row_size)
376
377: BaseBEllPackStructInfo(nrows, kcol[nrows])
378, m_block_nrows(BEllPackStructInfo<EllpackSize, IndexT>::nbBlocks(nrows))
379, m_block_nnz(h_block_row_offset[m_block_nrows])
380, m_h_local_row_size(h_local_row_size)
381{
382 // clang-format off
383 alien_debug([&] {
384 cout() << "COMPUTE BELLPACK PROFILE";
385 cout() << " NROWS : "<<this->m_nrows;
386 cout() << " NNZ : "<<this->m_nnz;
387 cout() << " BLOCK NROWS : "<<m_block_nrows;
388 cout() <<" BLOCK NNZ : "<<m_block_nnz;
389 });
390 // clang-format on
391
392 //m_h_block_cols.assign(m_block_nnz*ellpack_size,-1) ;
393 m_internal = new InternalType{ this->m_nrows,
394 this->m_nnz,
395 m_block_nrows,
396 m_block_nnz,
397 kcol,
398 cols,
399 h_block_row_offset,
400 h_local_row_size };
401}
402
403template <int EllpackSize, typename IndexT>
404BEllPackStructInfo<EllpackSize, IndexT>::~BEllPackStructInfo()
405{
406 delete m_internal ;
407}
408
409template <int EllpackSize, typename IndexT>
410typename BEllPackStructInfo<EllpackSize, IndexT>::IndexType const*
411BEllPackStructInfo<EllpackSize, IndexT>::kcol() const
412{
413 return m_internal->kcol();
414}
415
416template <int EllpackSize, typename IndexT>
417typename BEllPackStructInfo<EllpackSize, IndexT>::IndexType const*
418BEllPackStructInfo<EllpackSize, IndexT>::cols() const
419{
420 return m_internal->cols();
421}
422
423template <int EllpackSize, typename IndexT>
424typename BEllPackStructInfo<EllpackSize, IndexT>::IndexType const*
425BEllPackStructInfo<EllpackSize, IndexT>::dcol() const
426{
427 return m_internal->dcol();
428}
429
430namespace SYCLInternal
431{
432 template <typename ValueT, int EllpackSize>
433 MatrixInternal<ValueT, EllpackSize>::MatrixInternal(ProfileType const* profile, int N)
434 : m_profile(profile)
435 , m_N(N)
436 , m_NxN(N*N)
437 , m_h_values(profile->getBlockNnz() * ellpack_size * N*N)
438 , m_values(m_h_values.data(), sycl::range<1>(profile->getBlockNnz() * ellpack_size * N*N))
439 {
440 //m_values.set_final_data(nullptr);
441 //alien_debug([&] { cout() << "SYCL InternalMATRIX" << profile->getBlockNnz() * ellpack_size<< "N="<<m_N; });
442 }
443
444 template <typename ValueT, int EllpackSize>
445 bool MatrixInternal<ValueT, EllpackSize>::setMatrixValuesFromHost()
446 {
447 //alien_debug([&] { cout() << "SYCLMatrix setMatrixValuesFromHost "; });
448 auto env = SYCLEnv::instance();
449 auto& queue = env->internal()->queue();
450 auto num_groups = env->internal()->maxNumGroups();
451 auto max_work_group_size = env->internal()->maxWorkGroupSize();
452 auto total_threads = num_groups * ellpack_size;
453
454 int N = this->m_N;
455 int NxN = this->m_NxN;
456 int pack_size = ellpack_size ;
457 auto nrows = m_profile->getNRows();
458 auto nnz = m_profile->getNnz();
459 auto block_nnz = m_profile->getBlockNnz();
460
461 auto internal_profile = m_profile->internal();
462 auto& kcol = internal_profile->getKCol();
463 auto& block_row_offset = internal_profile->getBlockRowOffset();
464 auto& block_cols = internal_profile->getBlockCols();
465 auto local_row_size = m_profile->localRowSize();
466 if(N==1)
467 {
468 if (local_row_size == nullptr)
469 {
470 assert(m_h_csr_values.size()>=nnz) ;
471 ValueBufferType values_buffer(m_h_csr_values.data(), sycl::range<1>(nnz));
472 // COMPUTE COLS
473 // clang-format off
474 queue.submit([&](sycl::handler& cgh)
475 {
476 auto access_kcol_buffer = kcol.template get_access<sycl::access::mode::read>(cgh);
477 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
478 auto access_values_buffer = values_buffer.template get_access<sycl::access::mode::read>(cgh);
479 auto access_block_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
480
481 //cgh.parallel_for<class vector_axpy>(sycl::nd_range<1>{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}},[=](sycl::nd_item<1> item_id)
482 cgh.parallel_for<class set_matrix_values>(
483 sycl::range<1>{total_threads},
484 [=] (sycl::item<1> item_id)
485 {
486 auto id = item_id.get_id(0);
487 //auto local_id = item_id.get_local_id(0);
488 //auto block_id = item_id.get_group(0) ;
489 //auto global_id = item_id.get_global_id(0);
490
491 for (auto i = id; i < nrows; i += item_id.get_range()[0])
492 {
493 auto block_id = i/ellpack_size ;
494 auto local_id = i%ellpack_size ;
495
496 auto begin = access_kcol_buffer[i] ;
497 auto end = access_kcol_buffer[i+1] ;
498 auto row_size = end - begin ;
499
500 int block_row_offset = access_block_row_offset[block_id];
501 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
502
503 for(int k=0;k<row_size;++k)
504 {
505 access_block_values[(block_row_offset+k)*ellpack_size+local_id] = access_values_buffer[begin+k] ;
506 }
507 for(int k=row_size;k<block_row_size;++k)
508 {
509 access_block_values[(block_row_offset+k)*ellpack_size+local_id] = 0 ;
510 }
511 }
512 });
513 }) ;
514 // clang-format on
515 m_values_is_update = true;
516 }
517 else
518 {
519 ValueBufferType values_buffer(m_h_csr_values.data(), sycl::range<1>(nnz));
520 IndexBufferType lrowsize_buffer(local_row_size, sycl::range<1>(nrows));
521 // COMPUTE COLS
522 // clang-format off
523 queue.submit([&](sycl::handler& cgh)
524 {
525 auto access_kcol_buffer = internal_profile->getKCol().template get_access<sycl::access::mode::read>(cgh);
526 auto access_lrowsize_buffer = lrowsize_buffer.template get_access<sycl::access::mode::read>(cgh);
527
528 auto access_block_row_offset = internal_profile->getBlockRowOffset().template get_access<sycl::access::mode::read>(cgh);
529 auto access_values_buffer = values_buffer.template get_access<sycl::access::mode::read>(cgh);
530 auto access_block_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
531
532 //cgh.parallel_for<class vector_axpy>(sycl::nd_range<1>{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}},[=](sycl::nd_item<1> item_id)
533 cgh.parallel_for<class set_matrix_values2>(
534 sycl::range<1>{total_threads},
535 [=] (sycl::item<1> item_id)
536 {
537 auto id = item_id.get_id(0);
538 //auto local_id = item_id.get_local_id(0);
539 //auto block_id = item_id.get_group(0) ;
540 //auto global_id = item_id.get_global_id(0);
541
542 for (auto i = id; i < nrows; i += item_id.get_range()[0])
543 {
544 auto block_id = i/pack_size ;
545 auto local_id = i%pack_size ;
546
547 auto begin = access_kcol_buffer[i] ;
548 auto lrow_size = access_lrowsize_buffer[i] ;
549 auto end = begin + lrow_size ;
550
551 int block_row_offset = access_block_row_offset[block_id];
552 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
553
554 for(int k=0;k<lrow_size;++k)
555 {
556 access_block_values[(block_row_offset+k)*pack_size+local_id] = access_values_buffer[begin+k] ;
557 }
558 for(int k=lrow_size;k<block_row_size;++k)
559 {
560 access_block_values[(block_row_offset+k)*pack_size+local_id] = 0. ;
561 }
562 }
563 });
564 }) ;
565
566 auto interface_nrows = m_ext_profile->getNRows();
567 auto ext_nnz = m_ext_profile->getNnz();
568 auto ext_block_nnz = m_ext_profile->getBlockNnz();
569
570 auto ext_internal_profile = m_ext_profile->internal();
571 auto& ext_kcol = ext_internal_profile->getKCol();
572 auto& ext_block_row_offset = ext_internal_profile->getBlockRowOffset();
573 //m_h_ext_values.resize(ext_block_nnz) ;
574 m_ext_values.reset(new ValueBufferType(ext_block_nnz * ellpack_size * NxN)) ;
575 {
576 ValueBufferType ext_csr_values_buffer(m_h_csr_ext_values.data(), sycl::range<1>(ext_nnz * NxN));
577 queue.submit([&](sycl::handler& cgh)
578 {
579 auto access_kcol_buffer = ext_kcol.template get_access<sycl::access::mode::read>(cgh);
580 auto access_block_row_offset = ext_block_row_offset.template get_access<sycl::access::mode::read>(cgh);
581 auto access_values_buffer = ext_csr_values_buffer.template get_access<sycl::access::mode::read>(cgh);
582 //auto access_block_values = m_ext_values->template get_access<sycl::access::mode::read_write>(cgh);
583 auto access_block_values = sycl::accessor { *m_ext_values, cgh, sycl::write_only, sycl::property::no_init{}};
584
585 cgh.parallel_for<class set_matrix_values3>(
586 sycl::range<1>{total_threads},
587 [=] (sycl::item<1> item_id)
588 {
589 auto id = item_id.get_id(0);
590
591 for (auto i = id; i < interface_nrows; i += item_id.get_range()[0])
592 {
593 auto block_id = i/ellpack_size ;
594 auto local_id = i%ellpack_size ;
595 auto begin = access_kcol_buffer[i] ;
596 auto end = access_kcol_buffer[i+1] ;
597 auto row_size = end - begin ;
598
599 int block_row_offset = access_block_row_offset[block_id]*ellpack_size*NxN ;
600 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
601
602 for(int k=begin*NxN;k<end*NxN;++k)
603 {
604 access_block_values[block_row_offset+(k-begin*NxN)*ellpack_size+local_id] = access_values_buffer[k] ;
605 }
606 for(int k=row_size*NxN;k<block_row_size*NxN;++k)
607 {
608 access_block_values[block_row_offset+k*ellpack_size+local_id] = 0 ;
609 }
610 }
611 });
612 }).wait() ;
613 // clang-format on
614 }
615 m_values_is_update = true;
616 }
617 }
618 else
619 {
620 if (local_row_size == nullptr)
621 {
622 assert(m_h_csr_values.size()>=nnz*NxN) ;
623 ValueBufferType values_buffer(m_h_csr_values.data(), sycl::range<1>(nnz*NxN));
624 // COMPUTE COLS
625 // clang-format off
626 queue.submit([&](sycl::handler& cgh)
627 {
628 auto access_kcol_buffer = kcol.template get_access<sycl::access::mode::read>(cgh);
629 auto access_cols_buffer = block_cols.template get_access<sycl::access::mode::read>(cgh);
630 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
631 auto access_values_buffer = values_buffer.template get_access<sycl::access::mode::read>(cgh);
632 auto access_block_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
633
634 //cgh.parallel_for<class vector_axpy>(sycl::nd_range<1>{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}},[=](sycl::nd_item<1> item_id)
635 cgh.parallel_for<class set_block_matrix_values>(
636 sycl::range<1>{total_threads},
637 [=] (sycl::item<1> item_id)
638 {
639 auto id = item_id.get_id(0);
640 //auto local_id = item_id.get_local_id(0);
641 //auto block_id = item_id.get_group(0) ;
642 //auto global_id = item_id.get_global_id(0);
643
644 for (std::size_t i = id; i < nrows; i += item_id.get_range()[0])
645 {
646 auto block_id = i/pack_size ;
647 auto local_id = i%pack_size ;
648
649 int begin = access_kcol_buffer[i] ;
650 int end = access_kcol_buffer[i+1] ;
651 int row_size = end - begin ;
652
653 int block_row_offset = access_block_row_offset[block_id] ;
654 int block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
655
656 for(int k=0;k<row_size;++k)
657 {
658 for(int i=0;i<N;++i)
659 for(int j=0;j<N;++j)
660 {
661 access_block_values[((block_row_offset+k)*NxN+i*N+j)*pack_size+local_id] = access_values_buffer[(begin+k)*NxN + i*N+j] ;
662 }
663 /*
664 printf("MATRIX[%zu,%d] (%d,%d) = (%f,%f,%f,%f)\n",i,k,block_row_offset,begin, access_values_buffer[(begin+k)*NxN], access_values_buffer[(begin+k)*NxN+1], access_values_buffer[(begin+k)*NxN+2], access_values_buffer[(begin+k)*NxN +3]);*/
665
666 }
667 for(int k=row_size;k<block_row_size;++k)
668 {
669 for(int i=0;i<N;++i)
670 for(int j=0;j<N;++j)
671 {
672 access_block_values[((block_row_offset+k)*NxN + i*N+j)*pack_size+local_id] = 0. ;
673 }
674 }
675 }
676 });
677 }).wait() ;
678 // clang-format on
679 m_values_is_update = true;
680 /*
681 {
682 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> h_acc(m_values);
683 sycl::host_accessor<int, 1, sycl::access::mode::read> kcol_acc(block_row_offset);
684 for(int irow=0;irow<nrows;++irow)
685 {
686 auto ib=irow/ellpack_size ;
687 auto il=irow%ellpack_size ;
688 std::cout<<"MAT["<<irow<<","<<ib<<","<<il<<"]:";
689 for(int k=kcol_acc[ib];k<kcol_acc[ib+1];++k)
690 {
691 for(int i=0;i<N;++i)
692 for(int j=0;j<N;++j)
693 std::cout<<h_acc[(k*NxN+i*N+j)*ellpack_size+il]<<",";
694 std::cout<<std::endl;
695 }
696 }
697 }*/
698 }
699 else
700 {
701 /*
702 alien_info([&] {
703 auto h_kcol = internal_profile->kcol() ;
704 auto cols = internal_profile->cols() ;
705 for(int irow=0;irow<nrows;++irow)
706 {
707 cout()<<"ROW["<<irow<<"]:";
708 for(int k=h_kcol[irow];k<h_kcol[irow]+local_row_size[irow];++k)
709 {
710 cout()<<"\t COL "<<cols[k];
711 for(int i=0;i<N;++i)
712 for(int j=0;j<N;++j)
713 cout()<<"\t\t MAT["<<i<<" "<<j<<"] "<<m_h_csr_values[k*NxN+i*N+j];
714 }
715 }
716 });*/
717 ValueBufferType values_buffer(m_h_csr_values.data(), sycl::range<1>(nnz*NxN));
718 IndexBufferType lrowsize_buffer(local_row_size, sycl::range<1>(nrows));
719 // COMPUTE COLS
720 // clang-format off
721 queue.submit([&](sycl::handler& cgh)
722 {
723 auto access_kcol_buffer = internal_profile->getKCol().template get_access<sycl::access::mode::read>(cgh);
724 auto access_lrowsize_buffer = lrowsize_buffer.template get_access<sycl::access::mode::read>(cgh);
725
726 auto access_block_row_offset = internal_profile->getBlockRowOffset().template get_access<sycl::access::mode::read>(cgh);
727 auto access_values_buffer = values_buffer.template get_access<sycl::access::mode::read>(cgh);
728 auto access_block_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
729
730 //cgh.parallel_for<class vector_axpy>(sycl::nd_range<1>{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}},[=](sycl::nd_item<1> item_id)
731 cgh.parallel_for<class set_block_matrix_values2>(
732 sycl::range<1>{total_threads},
733 [=] (sycl::item<1> item_id)
734 {
735 auto id = item_id.get_id(0);
736 //auto local_id = item_id.get_local_id(0);
737 //auto block_id = item_id.get_group(0) ;
738 //auto global_id = item_id.get_global_id(0);
739
740 for (auto i = id; i < nrows; i += item_id.get_range()[0])
741 {
742 auto block_id = i/pack_size ;
743 auto local_id = i%pack_size ;
744
745 auto begin = access_kcol_buffer[i] ;
746 auto lrow_size = access_lrowsize_buffer[i] ;
747 auto end = begin + lrow_size ;
748
749 auto block_row_offset = access_block_row_offset[block_id];
750 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
751
752 for(int k=0;k<lrow_size;++k)
753 {
754 for(int i=0;i<N;++i)
755 for(int j=0;j<N;++j)
756 {
757 access_block_values[((block_row_offset+k)*NxN+i*N+j)*pack_size+local_id] = access_values_buffer[(begin+k)*NxN+i*N+j] ;
758 }
759 }
760 for(int k=lrow_size;k<block_row_size;++k)
761 {
762 for(int i=0;i<N;++i)
763 for(int j=0;j<N;++j)
764 {
765 access_block_values[((block_row_offset+k)*NxN+i*N+j)*pack_size+local_id] = 0. ;
766 }
767 }
768 }
769 });
770 });
771
772
773 auto interface_nrows = m_ext_profile->getNRows();
774 auto ext_nnz = m_ext_profile->getNnz();
775 auto ext_block_nnz = m_ext_profile->getBlockNnz();
776
777 auto ext_internal_profile = m_ext_profile->internal();
778 auto& ext_kcol = ext_internal_profile->getKCol();
779 auto& ext_block_row_offset = ext_internal_profile->getBlockRowOffset();
780 //m_h_ext_values.resize(ext_block_nnz) ;
781 m_ext_values.reset(new ValueBufferType(ext_block_nnz * ellpack_size * NxN)) ;
782 {
783 ValueBufferType ext_csr_values_buffer(m_h_csr_ext_values.data(), sycl::range<1>(ext_nnz * NxN));
784 queue.submit([&](sycl::handler& cgh)
785 {
786 auto access_kcol_buffer = ext_kcol.template get_access<sycl::access::mode::read>(cgh);
787 auto access_block_row_offset = ext_block_row_offset.template get_access<sycl::access::mode::read>(cgh);
788 auto access_values_buffer = ext_csr_values_buffer.template get_access<sycl::access::mode::read>(cgh);
789 //auto access_block_values = m_ext_values->template get_access<sycl::access::mode::read_write>(cgh);
790 auto access_block_values = sycl::accessor { *m_ext_values, cgh, sycl::write_only, sycl::property::no_init{}};
791
792 cgh.parallel_for<class set_matrix_values3>(
793 sycl::range<1>{total_threads},
794 [=] (sycl::item<1> item_id)
795 {
796 auto id = item_id.get_id(0);
797
798 for (auto i = id; i < interface_nrows; i += item_id.get_range()[0])
799 {
800 auto block_id = i/pack_size ;
801 auto local_id = i%pack_size ;
802 auto begin = access_kcol_buffer[i] ;
803 auto end = access_kcol_buffer[i+1] ;
804 auto row_size = end - begin ;
805
806 int block_row_offset = access_block_row_offset[block_id] ;
807 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
808
809 for(int k=0;k<row_size;++k)
810 {
811 for(int i=0;i<N;++i)
812 for(int j=0;j<N;++j)
813 {
814 access_block_values[((block_row_offset+k)*NxN+i*N+j)*pack_size+local_id] = access_values_buffer[(begin+k)*NxN+i*N+j] ;
815 }
816 }
817 for(int k=row_size;k<block_row_size;++k)
818 {
819 for(int i=0;i<N;++i)
820 for(int j=0;j<N;++j)
821 {
822 access_block_values[((block_row_offset+k)*NxN+i*N+j)*pack_size+local_id] = 0 ;
823 }
824 }
825 }
826 });
827 }).wait() ;
828 // clang-format on
829 }
830 /*
831 {
832 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> h_acc(m_values);
833 sycl::host_accessor<int, 1, sycl::access::mode::read> kcol_acc(block_row_offset);
834 sycl::host_accessor<int, 1, sycl::access::mode::read> cols_acc(block_cols);
835 for(int irow=0;irow<nrows;++irow)
836 {
837 auto ib=irow/pack_size ;
838 auto il=irow%pack_size ;
839 alien_info([&] {
840 cout()<<"MAT["<<irow<<","<<ib<<","<<il<<"]:";
841 for(int k=kcol_acc[ib];k<kcol_acc[ib+1];++k)
842 {
843 cout()<<"\t K"<<k<<" "<<cols_acc[k*pack_size+il];
844 for(int i=0;i<N;++i)
845 for(int j=0;j<N;++j)
846 cout()<<"\t\t"<<h_acc[(k*NxN+i*N+j)*pack_size+il]<<",";
847 //std::cout<<std::endl;
848 }
849 });
850 }
851 }*/
852 m_values_is_update = true;
853 }
854 }
855 return true;
856 }
857
858 template <typename ValueT, int EllpackSize>
859 bool MatrixInternal<ValueT, EllpackSize>::setMatrixValues(ValueBufferType& values_buffer)
860 {
861 //alien_debug([&] { cout() << "SYCLMatrix setMatrixValues "; });
862 auto env = SYCLEnv::instance();
863 auto& queue = env->internal()->queue();
864 auto num_groups = env->internal()->maxNumGroups();
865 auto max_work_group_size = env->internal()->maxWorkGroupSize();
866 auto total_threads = num_groups * ellpack_size;
867
868 int NxN = this->m_NxN;
869
870 auto nrows = m_profile->getNRows();
871 auto nnz = m_profile->getNnz();
872 auto block_nnz = m_profile->getBlockNnz();
873
874 auto internal_profile = m_profile->internal();
875 auto& kcol = internal_profile->getKCol();
876 auto& block_row_offset = internal_profile->getBlockRowOffset();
877
878 auto local_row_size = m_profile->localRowSize();
879 if (local_row_size == nullptr)
880 {
881 //ValueBufferType values_buffer(m_h_csr_values.data(), sycl::range<1>(nnz));
882 // COMPUTE COLS
883 // clang-format off
884 queue.submit([&](sycl::handler& cgh)
885 {
886 auto access_kcol_buffer = internal_profile->getKCol().template get_access<sycl::access::mode::read>(cgh);
887 auto access_block_row_offset = internal_profile->getBlockRowOffset().template get_access<sycl::access::mode::read>(cgh);
888 auto access_values_buffer = values_buffer.template get_access<sycl::access::mode::read>(cgh);
889 auto access_block_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
890
891 //cgh.parallel_for<class vector_axpy>(sycl::nd_range<1>{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}},[=](sycl::nd_item<1> item_id)
892 cgh.parallel_for<class set_matrix_values4>(sycl::range<1>{total_threads},
893 [=] (sycl::item<1> item_id)
894 {
895 auto id = item_id.get_id(0);
896 //auto local_id = item_id.get_local_id(0);
897 //auto block_id = item_id.get_group(0) ;
898 //auto global_id = item_id.get_global_id(0);
899
900 for (auto i = id; i < nrows; i += item_id.get_range()[0])
901 {
902 auto block_id = i/ellpack_size ;
903 auto local_id = i%ellpack_size ;
904
905 auto begin = access_kcol_buffer[i] ;
906 auto end = access_kcol_buffer[i+1] ;
907 auto row_size = end - begin ;
908
909 int block_row_offset = access_block_row_offset[block_id]*ellpack_size*NxN ;
910 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
911
912 for(int k=begin*NxN;k<end*NxN;++k)
913 {
914 access_block_values[block_row_offset+(k-begin*NxN)*ellpack_size+local_id] = access_values_buffer[k] ;
915 }
916 for(int k=row_size*NxN;k<block_row_size*NxN;++k)
917 {
918 access_block_values[block_row_offset+k*ellpack_size+local_id] = 0 ;
919 }
920 }
921 });
922 }) ;
923 // clang-format on
924 m_values_is_update = true;
925 }
926 else
927 {
928 //ValueBufferType values_buffer(m_h_csr_values.data(), sycl::range<1>(nnz));
929 IndexBufferType lrowsize_buffer(local_row_size, sycl::range<1>(nrows));
930
931 // COMPUTE COLS
932 // clang-format off
933 queue.submit([&](sycl::handler& cgh)
934 {
935 auto access_kcol_buffer = internal_profile->getKCol().template get_access<sycl::access::mode::read>(cgh);
936 auto access_lrowsize_buffer = lrowsize_buffer.template get_access<sycl::access::mode::read>(cgh);
937
938 auto access_block_row_offset = internal_profile->getBlockRowOffset().template get_access<sycl::access::mode::read>(cgh);
939 auto access_values_buffer = values_buffer.template get_access<sycl::access::mode::read>(cgh);
940 auto access_block_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
941
942 //cgh.parallel_for<class vector_axpy>(sycl::nd_range<1>{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}},[=](sycl::nd_item<1> item_id)
943 cgh.parallel_for<class set_matrix_values5>(sycl::range<1>{total_threads},
944 [=] (sycl::item<1> item_id)
945 {
946 auto id = item_id.get_id(0);
947 //auto local_id = item_id.get_local_id(0);
948 //auto block_id = item_id.get_group(0) ;
949 //auto global_id = item_id.get_global_id(0);
950
951 for (auto i = id; i < nrows; i += item_id.get_range()[0])
952 {
953 auto block_id = i/ellpack_size ;
954 auto local_id = i%ellpack_size ;
955
956 auto begin = access_kcol_buffer[i] ;
957 auto lrow_size = access_lrowsize_buffer[i] ;
958 auto end = begin + lrow_size ;
959
960 int block_row_offset = access_block_row_offset[block_id]*ellpack_size*NxN ;
961 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
962
963 for(int k=begin*NxN;k<end*NxN;++k)
964 {
965 access_block_values[block_row_offset+(k-begin*NxN)*ellpack_size+local_id] = access_values_buffer[k] ;
966 }
967 for(int k=lrow_size*NxN;k<block_row_size*NxN;++k)
968 {
969 access_block_values[block_row_offset+k*ellpack_size+local_id] = 0 ;
970 }
971 }
972 });
973 }) ;
974 auto kcol = m_profile->kcol() ;
975 auto local_row_size = m_profile->localRowSize() ;
976
977 auto interface_nrows = m_ext_profile->getNRows();
978 auto ext_nnz = m_ext_profile->getNnz();
979 auto ext_block_nnz = m_ext_profile->getBlockNnz();
980
981 auto ext_internal_profile = m_ext_profile->internal();
982 auto& ext_kcol = ext_internal_profile->getKCol();
983 auto& ext_block_row_offset = ext_internal_profile->getBlockRowOffset();
984 //m_h_ext_values.resize(ext_block_nnz) ;
985 m_ext_values.reset(new ValueBufferType(ext_block_nnz * ellpack_size * NxN)) ;
986
987 // EXTRACT EXTERNAL PROFILE
988 {
989 ValueBufferType ext_csr_values_buffer(m_h_csr_ext_values.data(), sycl::range<1>(ext_nnz * NxN));
990
991 std::vector<Integer> h_local_row_offset(interface_nrows+1) ;
992 {
993 Integer offset = 0 ;
994 for(std::size_t i=0;i<interface_nrows;++i)
995 {
996 h_local_row_offset[i] = offset ;
997 offset += local_row_size[i] ;
998 }
999 h_local_row_offset[interface_nrows] = offset ;
1000 }
1001 IndexBufferType local_row_offset(h_local_row_offset.data(),sycl::range<1>(interface_nrows+1)) ;
1002
1003 queue.submit([&](sycl::handler& cgh)
1004 {
1005 auto access_kcol = internal_profile->getKCol().template get_access<sycl::access::mode::read>(cgh);
1006 auto access_interface_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
1007 auto access_local_row_size = lrowsize_buffer.template get_access<sycl::access::mode::read>(cgh);
1008 auto access_local_row_offset = local_row_offset.template get_access<sycl::access::mode::read>(cgh);
1009 auto access_csr_values = values_buffer.template get_access<sycl::access::mode::read>(cgh);
1010 auto access_ext_csr_values = ext_csr_values_buffer.template get_access<sycl::access::mode::read_write>(cgh);
1011 cgh.parallel_for<class set_matrix_values6>(sycl::range<1>{total_threads},
1012 [=] (sycl::item<1> item_id)
1013 {
1014 auto id = item_id.get_id(0);
1015 for (auto i = id; i < interface_nrows; i += item_id.get_range()[0])
1016 {
1017 Integer jcol = access_local_row_offset[i] ;
1018 for (int k = (access_kcol[i] + access_local_row_size[i])*NxN; k < access_kcol[i + 1] * NxN; ++k)
1019 access_ext_csr_values[jcol++] = access_csr_values[k];
1020 }
1021 });
1022 }).wait() ;
1023
1024 queue.submit([&](sycl::handler& cgh)
1025 {
1026 auto access_kcol_buffer = ext_kcol.template get_access<sycl::access::mode::read>(cgh);
1027 auto access_block_row_offset = ext_block_row_offset.template get_access<sycl::access::mode::read>(cgh);
1028 auto access_values_buffer = ext_csr_values_buffer.template get_access<sycl::access::mode::read>(cgh);
1029 //auto access_block_values = m_ext_values->template get_access<sycl::access::mode::read_write>(cgh);
1030 auto access_block_values = sycl::accessor { *m_ext_values, cgh, sycl::write_only, sycl::property::no_init{}};
1031
1032 cgh.parallel_for<class set_matrix_values7>(sycl::range<1>{total_threads},
1033 [=] (sycl::item<1> item_id)
1034 {
1035 auto id = item_id.get_id(0);
1036
1037 for (auto i = id; i < interface_nrows; i += item_id.get_range()[0])
1038 {
1039 auto block_id = i/ellpack_size ;
1040 auto local_id = i%ellpack_size ;
1041 auto begin = access_kcol_buffer[i] ;
1042 auto end = access_kcol_buffer[i+1] ;
1043 auto row_size = end - begin ;
1044
1045 int block_row_offset = access_block_row_offset[block_id]*ellpack_size*NxN ;
1046 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
1047
1048 for(int k=begin*NxN;k<end*NxN;++k)
1049 {
1050 access_block_values[block_row_offset+(k-begin*NxN)*ellpack_size+local_id] = access_values_buffer[k] ;
1051 }
1052 for(int k=row_size*NxN;k<block_row_size*NxN;++k)
1053 {
1054 access_block_values[block_row_offset+k*ellpack_size+local_id] = 0 ;
1055 }
1056 }
1057 });
1058 }).wait() ;
1059 // clang-format on
1060 }
1061 m_values_is_update = true;
1062 }
1063 return true;
1064 }
1065
1066 template <typename ValueT, int EllpackSize>
1067 bool MatrixInternal<ValueT, EllpackSize>::setMatrixValues(ValueBufferType& values_buffer,
1068 ValueBufferType& ext_values_buffer)
1069 {
1070 //alien_debug([&] { cout() << "SYCLMatrix setMatrixValues "; });
1071 auto env = SYCLEnv::instance();
1072 auto& queue = env->internal()->queue();
1073 auto num_groups = env->internal()->maxNumGroups();
1074 auto max_work_group_size = env->internal()->maxWorkGroupSize();
1075 auto total_threads = num_groups * ellpack_size;
1076
1077 int NxN = this->m_NxN;
1078
1079 auto nrows = m_profile->getNRows();
1080 auto nnz = m_profile->getNnz();
1081 auto block_nnz = m_profile->getBlockNnz();
1082
1083 auto internal_profile = m_profile->internal();
1084 auto& kcol = internal_profile->getKCol();
1085 auto& block_row_offset = internal_profile->getBlockRowOffset();
1086
1087 auto local_row_size = m_profile->localRowSize();
1088 if (local_row_size == nullptr)
1089 {
1090 //ValueBufferType values_buffer(m_h_csr_values.data(), sycl::range<1>(nnz));
1091 // COMPUTE COLS
1092 // clang-format off
1093 queue.submit([&](sycl::handler& cgh)
1094 {
1095 auto access_kcol_buffer = internal_profile->getKCol().template get_access<sycl::access::mode::read>(cgh);
1096 auto access_block_row_offset = internal_profile->getBlockRowOffset().template get_access<sycl::access::mode::read>(cgh);
1097 auto access_values_buffer = values_buffer.template get_access<sycl::access::mode::read>(cgh);
1098 auto access_block_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
1099
1100 //cgh.parallel_for<class vector_axpy>(sycl::nd_range<1>{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}},[=](sycl::nd_item<1> item_id)
1101 cgh.parallel_for<class set_matrix_values4>(sycl::range<1>{total_threads},
1102 [=] (sycl::item<1> item_id)
1103 {
1104 auto id = item_id.get_id(0);
1105 //auto local_id = item_id.get_local_id(0);
1106 //auto block_id = item_id.get_group(0) ;
1107 //auto global_id = item_id.get_global_id(0);
1108
1109 for (auto i = id; i < nrows; i += item_id.get_range()[0])
1110 {
1111 auto block_id = i/ellpack_size ;
1112 auto local_id = i%ellpack_size ;
1113
1114 auto begin = access_kcol_buffer[i] ;
1115 auto end = access_kcol_buffer[i+1] ;
1116 auto row_size = end - begin ;
1117
1118 int block_row_offset = access_block_row_offset[block_id]*ellpack_size*NxN ;
1119 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
1120
1121 for(int k=begin*NxN;k<end*NxN;++k)
1122 {
1123 access_block_values[block_row_offset+(k-begin*NxN)*ellpack_size+local_id] = access_values_buffer[k] ;
1124 }
1125 for(int k=row_size*NxN;k<block_row_size*NxN;++k)
1126 {
1127 access_block_values[block_row_offset+k*ellpack_size+local_id] = 0 ;
1128 }
1129 }
1130 });
1131 }) ;
1132 // clang-format on
1133 m_values_is_update = true;
1134 }
1135 else
1136 {
1137 //ValueBufferType values_buffer(m_h_csr_values.data(), sycl::range<1>(nnz));
1138 IndexBufferType lrowsize_buffer(local_row_size, sycl::range<1>(nrows));
1139
1140 // COMPUTE COLS
1141 // clang-format off
1142 queue.submit([&](sycl::handler& cgh)
1143 {
1144 auto access_kcol_buffer = internal_profile->getKCol().template get_access<sycl::access::mode::read>(cgh);
1145 auto access_lrowsize_buffer = lrowsize_buffer.template get_access<sycl::access::mode::read>(cgh);
1146
1147 auto access_block_row_offset = internal_profile->getBlockRowOffset().template get_access<sycl::access::mode::read>(cgh);
1148 auto access_values_buffer = values_buffer.template get_access<sycl::access::mode::read>(cgh);
1149 auto access_block_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
1150
1151 //cgh.parallel_for<class vector_axpy>(sycl::nd_range<1>{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}},[=](sycl::nd_item<1> item_id)
1152 cgh.parallel_for<class set_matrix_values5>(sycl::range<1>{total_threads},
1153 [=] (sycl::item<1> item_id)
1154 {
1155 auto id = item_id.get_id(0);
1156 //auto local_id = item_id.get_local_id(0);
1157 //auto block_id = item_id.get_group(0) ;
1158 //auto global_id = item_id.get_global_id(0);
1159
1160 for (auto i = id; i < nrows; i += item_id.get_range()[0])
1161 {
1162 auto block_id = i/ellpack_size ;
1163 auto local_id = i%ellpack_size ;
1164
1165 auto begin = access_kcol_buffer[i] ;
1166 auto lrow_size = access_lrowsize_buffer[i] ;
1167 auto end = begin + lrow_size ;
1168
1169 int block_row_offset = access_block_row_offset[block_id]*ellpack_size*NxN ;
1170 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
1171
1172 for(int k=begin*NxN;k<end*NxN;++k)
1173 {
1174 access_block_values[block_row_offset+(k-begin*NxN)*ellpack_size+local_id] = access_values_buffer[k] ;
1175 }
1176 for(int k=lrow_size*NxN;k<block_row_size*NxN;++k)
1177 {
1178 access_block_values[block_row_offset+k*ellpack_size+local_id] = 0 ;
1179 }
1180 }
1181 });
1182 }) ;
1183 auto kcol = m_profile->kcol() ;
1184 auto local_row_size = m_profile->localRowSize() ;
1185
1186 auto interface_nrows = m_ext_profile->getNRows();
1187 auto ext_nnz = m_ext_profile->getNnz();
1188 auto ext_block_nnz = m_ext_profile->getBlockNnz();
1189
1190 auto ext_internal_profile = m_ext_profile->internal();
1191 auto& ext_kcol = ext_internal_profile->getKCol();
1192 auto& ext_block_row_offset = ext_internal_profile->getBlockRowOffset();
1193 //m_h_ext_values.resize(ext_block_nnz) ;
1194 m_ext_values.reset(new ValueBufferType(ext_block_nnz * ellpack_size * NxN)) ;
1195
1196 // EXTRACT EXTERNAL PROFILE
1197 {
1198 ValueBufferType ext_csr_values_buffer(m_h_csr_ext_values.data(), sycl::range<1>(ext_nnz * NxN));
1199
1200 std::vector<Integer> h_local_row_offset(interface_nrows+1) ;
1201 {
1202 Integer offset = 0 ;
1203 for(std::size_t i=0;i<interface_nrows;++i)
1204 {
1205 h_local_row_offset[i] = offset ;
1206 offset += local_row_size[i] ;
1207 }
1208 h_local_row_offset[interface_nrows] = offset ;
1209 }
1210 IndexBufferType local_row_offset(h_local_row_offset.data(),sycl::range<1>(interface_nrows+1)) ;
1211
1212 queue.submit([&](sycl::handler& cgh)
1213 {
1214 auto access_kcol = internal_profile->getKCol().template get_access<sycl::access::mode::read>(cgh);
1215 auto access_interface_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
1216 auto access_local_row_size = lrowsize_buffer.template get_access<sycl::access::mode::read>(cgh);
1217 auto access_local_row_offset = local_row_offset.template get_access<sycl::access::mode::read>(cgh);
1218 auto access_csr_values = values_buffer.template get_access<sycl::access::mode::read>(cgh);
1219 auto access_ext_csr_values = ext_csr_values_buffer.template get_access<sycl::access::mode::read_write>(cgh);
1220 cgh.parallel_for<class set_matrix_values6>(sycl::range<1>{total_threads},
1221 [=] (sycl::item<1> item_id)
1222 {
1223 auto id = item_id.get_id(0);
1224 for (auto i = id; i < interface_nrows; i += item_id.get_range()[0])
1225 {
1226 Integer jcol = access_local_row_offset[i] ;
1227 for (int k = (access_kcol[i] + access_local_row_size[i])*NxN; k < access_kcol[i + 1] * NxN; ++k)
1228 access_ext_csr_values[jcol++] = access_csr_values[k];
1229 }
1230 });
1231 }).wait() ;
1232
1233 queue.submit([&](sycl::handler& cgh)
1234 {
1235 auto access_kcol_buffer = ext_kcol.template get_access<sycl::access::mode::read>(cgh);
1236 auto access_block_row_offset = ext_block_row_offset.template get_access<sycl::access::mode::read>(cgh);
1237 auto access_values_buffer = ext_csr_values_buffer.template get_access<sycl::access::mode::read>(cgh);
1238 //auto access_block_values = m_ext_values->template get_access<sycl::access::mode::read_write>(cgh);
1239 auto access_block_values = sycl::accessor { *m_ext_values, cgh, sycl::write_only, sycl::property::no_init{}};
1240
1241 cgh.parallel_for<class set_matrix_values7>(sycl::range<1>{total_threads},
1242 [=] (sycl::item<1> item_id)
1243 {
1244 auto id = item_id.get_id(0);
1245
1246 for (auto i = id; i < interface_nrows; i += item_id.get_range()[0])
1247 {
1248 auto block_id = i/ellpack_size ;
1249 auto local_id = i%ellpack_size ;
1250 auto begin = access_kcol_buffer[i] ;
1251 auto end = access_kcol_buffer[i+1] ;
1252 auto row_size = end - begin ;
1253
1254 int block_row_offset = access_block_row_offset[block_id]*ellpack_size*NxN ;
1255 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
1256
1257 for(int k=begin*NxN;k<end*NxN;++k)
1258 {
1259 access_block_values[block_row_offset+(k-begin*NxN)*ellpack_size+local_id] = access_values_buffer[k] ;
1260 }
1261 for(int k=row_size*NxN;k<block_row_size*NxN;++k)
1262 {
1263 access_block_values[block_row_offset+k*ellpack_size+local_id] = 0 ;
1264 }
1265 }
1266 });
1267 }).wait() ;
1268 // clang-format on
1269 }
1270 m_values_is_update = true;
1271 }
1272 return true;
1273 }
1274
1275 template <typename ValueT, int EllpackSize>
1276 bool MatrixInternal<ValueT, EllpackSize>::setMatrixValues(ValueT const* values, bool only_host)
1277 {
1278 //alien_debug([&] { cout() << "SYCLMatrix setMatrixValues " << only_host; });
1279
1280 int NxN = this->m_NxN;
1281
1282 auto nnz = m_profile->getNnz();
1283 m_h_csr_values.resize(nnz*NxN);
1284 std::copy(values, values + nnz*NxN, m_h_csr_values.begin());
1285 if (m_ext_profile) {
1286 // clang-format off
1287 auto kcol = m_profile->kcol() ;
1288 auto local_row_size = m_profile->localRowSize() ;
1289
1290 auto interface_nrows = m_ext_profile->getNRows() ;
1291 auto ext_nnz = m_ext_profile->getNnz();
1292 // clang-format on
1293
1294 m_h_csr_ext_values.resize(ext_nnz*NxN);
1295
1296 // EXTRACT EXTERNAL PROFILE
1297 {
1298 int jcol = 0;
1299 for (std::size_t i = 0; i < interface_nrows; ++i) {
1300 auto id = m_h_interface_row_ids[i];
1301 for (int k = (kcol[id] + local_row_size[id])*NxN; k < kcol[id + 1]*NxN; ++k) {
1302 m_h_csr_ext_values[jcol++] = values[k];
1303 }
1304 }
1305 }
1306 }
1307 if (only_host) {
1308 m_values_is_update = false;
1309 }
1310 else {
1311 setMatrixValuesFromHost();
1312 m_values_is_update = true;
1313 }
1314 return true;
1315 }
1316
1317 template <typename ValueT, int EllpackSize>
1318 bool MatrixInternal<ValueT, EllpackSize>::copy(std::size_t nb_blocks,
1319 int block_size,
1320 ValueBufferType& values_buffer,
1321 int rhs_block_size)
1322 {
1323 auto env = SYCLEnv::instance();
1324 auto& queue = env->internal()->queue();
1325 auto num_groups = env->internal()->maxNumGroups();
1326 auto max_work_group_size = env->internal()->maxWorkGroupSize();
1327 auto total_threads = num_groups * ellpack_size;
1328
1329 int pack_size = ellpack_size;
1330 int N1 = block_size ;
1331 int N1xN1 = N1*N1 ;
1332 int N2 = rhs_block_size ;
1333 int N2xN2 = N2*N2 ;
1334 auto nrows = m_profile->getNRows();
1335 auto nnz = m_profile->getNnz();
1336 auto block_nnz = m_profile->getBlockNnz();
1337
1338 auto internal_profile = m_profile->internal();
1339 auto& kcol = internal_profile->getKCol();
1340 auto& block_row_offset = internal_profile->getBlockRowOffset();
1341
1342 auto local_row_size = m_profile->localRowSize();
1343 if (local_row_size == nullptr)
1344 {
1345 //ValueBufferType values_buffer(m_h_csr_values.data(), sycl::range<1>(nnz));
1346 // COMPUTE COLS
1347 // clang-format off
1348 queue.submit([&](sycl::handler& cgh)
1349 {
1350 auto access_kcol_buffer = internal_profile->getKCol().template get_access<sycl::access::mode::read>(cgh);
1351 auto access_block_row_offset = internal_profile->getBlockRowOffset().template get_access<sycl::access::mode::read>(cgh);
1352 auto access_values_buffer = values_buffer.template get_access<sycl::access::mode::read>(cgh);
1353 auto access_block_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
1354
1355 //cgh.parallel_for<class vector_axpy>(sycl::nd_range<1>{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}},[=](sycl::nd_item<1> item_id)
1356 cgh.parallel_for<class set_matrix_values4>(
1357 sycl::range<1>{total_threads},
1358 [=] (sycl::item<1> item_id)
1359 {
1360 auto id = item_id.get_id(0);
1361 //auto local_id = item_id.get_local_id(0);
1362 //auto block_id = item_id.get_group(0) ;
1363 //auto global_id = item_id.get_global_id(0);
1364
1365 for (auto i = id; i < nrows; i += item_id.get_range()[0])
1366 {
1367 auto block_id = i/pack_size ;
1368 auto local_id = i%pack_size ;
1369
1370 auto begin = access_block_row_offset[block_id] ;
1371 auto end = access_block_row_offset[block_id+1] ;
1372
1373 for(int k=begin;k<end;++k)
1374 {
1375 for(int i=0;i<N1;++i)
1376 for(int j=0;j<N1;++j)
1377 access_block_values[(k*N1xN1 +i*N1+j)*pack_size+local_id] = access_values_buffer[(k*N2xN2+i*N2+j)*pack_size+local_id] ;
1378 }
1379 }
1380 });
1381 }) ;
1382 // clang-format on
1383 }
1384 else
1385 {
1386 //ValueBufferType values_buffer(m_h_csr_values.data(), sycl::range<1>(nnz));
1387 IndexBufferType lrowsize_buffer(local_row_size, sycl::range<1>(nrows));
1388
1389 // COMPUTE COLS
1390 // clang-format off
1391 queue.submit([&](sycl::handler& cgh)
1392 {
1393 auto access_kcol_buffer = internal_profile->getKCol().template get_access<sycl::access::mode::read>(cgh);
1394 auto access_lrowsize_buffer = lrowsize_buffer.template get_access<sycl::access::mode::read>(cgh);
1395
1396 auto access_block_row_offset = internal_profile->getBlockRowOffset().template get_access<sycl::access::mode::read>(cgh);
1397 auto access_values_buffer = values_buffer.template get_access<sycl::access::mode::read>(cgh);
1398 auto access_block_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
1399
1400 //cgh.parallel_for<class vector_axpy>(sycl::nd_range<1>{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}},[=](sycl::nd_item<1> item_id)
1401 cgh.parallel_for<class set_matrix_values5>(
1402 sycl::range<1>{total_threads},
1403 [=] (sycl::item<1> item_id)
1404 {
1405 auto id = item_id.get_id(0);
1406 //auto local_id = item_id.get_local_id(0);
1407 //auto block_id = item_id.get_group(0) ;
1408 //auto global_id = item_id.get_global_id(0);
1409
1410 for (auto i = id; i < nrows; i += item_id.get_range()[0])
1411 {
1412 auto block_id = i/pack_size ;
1413 auto local_id = i%pack_size ;
1414
1415
1416 int block_row_offset = access_block_row_offset[block_id] ;
1417
1418 auto begin = access_block_row_offset[block_id] ;
1419 auto end = access_block_row_offset[block_id+1] ;
1420 auto lrow_size = access_lrowsize_buffer[i] ;
1421
1422 for(int k=begin;k<begin+lrow_size;++k)
1423 {
1424 for(int i=0;i<N1;++i)
1425 for(int j=0;j<N1;++j)
1426 access_block_values[(k*N1xN1 +i*N1+j)*pack_size+local_id] = access_values_buffer[(k*N2xN2+i*N2+j)*pack_size+local_id] ;
1427 }
1428 for(int k=begin+lrow_size;k<end;++k)
1429 {
1430 for(int i=0;i<N1;++i)
1431 for(int j=0;j<N1;++j)
1432 access_block_values[(k*N1xN1 +i*N1+j)*pack_size+local_id] = 0 ;
1433 }
1434 }
1435 });
1436 }) ;
1437 // clang-format on
1438 }
1439 m_values_is_update = true;
1440 return true;
1441 }
1442
1443 template <typename ValueT, int EllpackSize>
1444 bool MatrixInternal<ValueT, EllpackSize>::copy(std::size_t nb_blocks,
1445 int block_size,
1446 ValueBufferType& values_buffer,
1447 ValueBufferType& ext_values_buffer,
1448 int rhs_block_size)
1449 {
1450 copy(nb_blocks,block_size,values_buffer,rhs_block_size) ;
1451
1452 auto env = SYCLEnv::instance();
1453 auto& queue = env->internal()->queue();
1454 auto num_groups = env->internal()->maxNumGroups();
1455 auto max_work_group_size = env->internal()->maxWorkGroupSize();
1456 auto total_threads = num_groups * ellpack_size;
1457
1458 int pack_size = ellpack_size;
1459 int N1 = block_size ;
1460 int N1xN1 = N1*N1 ;
1461 int N2 = rhs_block_size ;
1462 int N2xN2 = N2*N2 ;
1463 auto nrows = m_profile->getNRows();
1464 auto nnz = m_profile->getNnz();
1465 auto block_nnz = m_profile->getBlockNnz();
1466
1467 auto internal_profile = m_profile->internal();
1468 auto& kcol = internal_profile->getKCol();
1469 auto& block_row_offset = internal_profile->getBlockRowOffset();
1470
1471 assert(m_ext_profile) ;
1472 auto interface_nrows = m_ext_profile->getNRows();
1473 auto ext_nnz = m_ext_profile->getNnz();
1474 auto ext_block_nnz = m_ext_profile->getBlockNnz();
1475 auto ext_internal_profile = m_ext_profile->internal();
1476 auto& ext_kcol = ext_internal_profile->getKCol();
1477 auto& ext_block_row_offset = ext_internal_profile->getBlockRowOffset();
1478
1479 m_h_ext_values.resize(ext_block_nnz*N1xN1) ;
1480 m_ext_values.reset(new ValueBufferType(ext_block_nnz * N1xN1* ellpack_size)) ;
1481 {
1482 queue.submit([&](sycl::handler& cgh)
1483 {
1484 auto access_block_row_offset = ext_block_row_offset.template get_access<sycl::access::mode::read>(cgh);
1485 auto access_ext_values_buffer = ext_values_buffer.template get_access<sycl::access::mode::read>(cgh);
1486 auto access_ext_values = sycl::accessor { *m_ext_values, cgh, sycl::write_only, sycl::property::no_init{}};
1487
1488 cgh.parallel_for<class set_matrix_ext_values3>(
1489 sycl::range<1>{total_threads},
1490 [=] (sycl::item<1> item_id)
1491 {
1492 auto id = item_id.get_id(0);
1493 for (auto i = id; i < interface_nrows; i += item_id.get_range()[0])
1494 {
1495 auto block_id = i/pack_size ;
1496 auto local_id = i%pack_size ;
1497
1498 auto begin = access_block_row_offset[block_id] ;
1499 auto end = access_block_row_offset[block_id+1] ;
1500
1501 for(int k=begin;k<begin;++k)
1502 {
1503 for(int i=0;i<N1;++i)
1504 for(int j=0;j<N1;++j)
1505 access_ext_values[(k*N1xN1 +i*N1+j)*pack_size+local_id] = access_ext_values_buffer[(k*N2xN2+i*N2+j)*pack_size+local_id] ;
1506 }
1507 }
1508 });
1509 }) ;
1510 // clang-format on
1511 }
1512 m_values_is_update = true;
1513 return true;
1514 }
1515
1516
1517 template <typename ValueT, int EllpackSize>
1518 bool MatrixInternal<ValueT, EllpackSize>::needUpdate()
1519 {
1520 return m_values_is_update != true;
1521 }
1522
1523 template <typename ValueT, int EllpackSize>
1524 void MatrixInternal<ValueT, EllpackSize>::notifyChanges()
1525 {
1526 m_values_is_update = false;
1527 }
1528
1529 template <typename ValueT, int EllpackSize>
1530 void MatrixInternal<ValueT, EllpackSize>::endUpdate()
1531 {
1532 if (not m_values_is_update) {
1533 setMatrixValuesFromHost();
1534 }
1535 }
1536
1537 /*
1538 template<int N>
1539 template <typename ValueT, int EllPackSize>
1540 void MatrixInternal<ValueT, EllPackSize>::multN(ValueBufferType& x, ValueBufferType& y, sycl::queue& queue) const
1541 {
1542 auto device = queue.get_device();
1543
1544 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
1545 // getting the maximum work group size per thread
1546 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
1547 // building the best number of global thread
1548
1549 // clang-format off
1550 std::size_t pack_size = ellpack_size;
1551 auto nrows = m_profile->getNRows();
1552 auto nnz = m_profile->getNnz();
1553
1554 auto internal_profile = m_profile->internal();
1555 auto& kcol = internal_profile->getKCol();
1556 auto& block_row_offset = internal_profile->getBlockRowOffset();
1557 auto& block_cols = internal_profile->getBlockCols();
1558
1559 auto blocks_needed = (nrows + ellpack_size - 1) / ellpack_size;
1560 auto blocks_target = std::max(blocks_needed, num_groups * 4UL);
1561 auto total_threads = blocks_target * pack_size;
1562
1563 queue.submit(
1564 [&](sycl::handler& cgh)
1565 {
1566 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
1567 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
1568 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
1569
1570 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
1571 auto access_y = y.template get_access<sycl::access::mode::discard_write>(cgh);
1572
1573 auto tile = TileT<N>() ;
1574
1575 sycl::local_accessor<ValueType, 1> lds_x{pack_size*N, cgh};
1576 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{pack_size}};
1577 cgh.parallel_for<class compute_mult>(r,
1578 [=](sycl::nd_item<1> item_id)
1579 {
1580 auto local_id = item_id.get_local_id(0);
1581 auto global_id = item_id.get_global_id(0);
1582
1583 for (auto i = global_id; i < nrows; i += item_id.get_global_range()[0])
1584 {
1585 auto block_id = i/pack_size ;
1586
1587 int begin = access_block_row_offset[block_id] ;
1588 int end = access_block_row_offset[block_id+1] ;
1589 for(int ieq=0;ieq<N;++ieq)
1590 {
1591 ValueType value = 0. ;
1592 for(int k=begin;k<end;++k)
1593 {
1594 //auto k = block_row_offset+j*ellpack_size+local_id ;
1595 const int col = access_cols[k * pack_size + local_id];
1596 if(col>=0)
1597 for(int ju=0;ju<N;++ju)
1598 lds_x[N*local_id+ju] = access_x[col*N+ju];
1599 item_id.barrier(sycl::access::fence_space::local_space);
1600 if(col>=0)
1601 {
1602 for(int ju=0;ju<N;++ju)
1603 value += access_values[tile.ijk(k,ieq,ju) + local_id] * lds_x[local_id*N+ju] ;
1604 }
1605 item_id.barrier(sycl::access::fence_space::local_space);
1606 }
1607 access_y[i*N+ieq] = value ;
1608 }
1609 }
1610 });
1611 });
1612 }*/
1613
1614 template <typename ValueT, int EllPackSize>
1615 void MatrixInternal<ValueT, EllPackSize>::mult(ValueBufferType& x, ValueBufferType& y, sycl::queue& queue) const
1616 {
1617 // clang-format on
1618 if(this->m_N==1)
1619 {
1620 auto device = queue.get_device();
1621
1622 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
1623 // getting the maximum work group size per thread
1624 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
1625 // building the best number of global thread
1626 //auto total_threads = num_groups * ellpack_size;
1627
1628 // clang-format off
1629 //int NxN = this->m_NxN;
1630 std::size_t pack_size = ellpack_size;
1631 auto nrows = m_profile->getNRows();
1632 auto nnz = m_profile->getNnz();
1633
1634 auto internal_profile = m_profile->internal();
1635 auto& kcol = internal_profile->getKCol();
1636 auto& block_row_offset = internal_profile->getBlockRowOffset();
1637 auto& block_cols = internal_profile->getBlockCols();
1638
1639 auto blocks_needed = (nrows + ellpack_size - 1) / ellpack_size;
1640 auto blocks_target = std::max(blocks_needed, num_groups * 4UL);
1641 auto total_threads = blocks_target * ellpack_size;
1642
1643 // COMPUTE VALUES
1644 // clang-format off
1645 queue.submit([&](sycl::handler& cgh)
1646 {
1647 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
1648 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
1649 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
1650
1651
1652 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
1653 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
1654#ifdef OLD
1655 //sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}};
1656 //cgh.parallel_for<class compute_mult>(r, [&](sycl::nd_item<1> item_id)
1657 cgh.parallel_for<class compute_mult>(
1658 sycl::range<1>{total_threads},
1659 [=] (sycl::item<1> item_id)
1660 {
1661 auto id = item_id.get_id(0);
1662 //auto local_id = item_id.get_local_id(0);
1663 //auto block_id = item_id.get_group(0) ;
1664 //auto global_id = item_id.get_global_id(0);
1665
1666 for (auto i = id; i < nrows; i += item_id.get_range()[0])
1667 {
1668 auto block_id = i/ellpack_size ;
1669 auto local_id = i%ellpack_size ;
1670 int block_row_offset = access_block_row_offset[block_id]*ellpack_size ;
1671 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
1672
1673 ValueType value = 0. ;
1674 for(int j=0;j<block_row_size;++j)
1675 {
1676 auto k = block_row_offset+j*ellpack_size+local_id ;
1677 value += access_values[k]* access_x[access_cols[k]] ;
1678 }
1679 access_y[i] = value ;
1680 }
1681 });
1682
1683#endif
1684 sycl::local_accessor<ValueType, 1> lds_x{pack_size, cgh};
1685 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{pack_size}};
1686 cgh.parallel_for<class compute_mult>(r,
1687 [=](sycl::nd_item<1> item_id)
1688 {
1689 //auto id = item_id.get_id(0);
1690 auto local_id = item_id.get_local_id(0);
1691 //auto block_id = item_id.get_group(0) ;
1692 auto global_id = item_id.get_global_id(0);
1693
1694 //for (auto i = id; i < nrows; i += item_id.get_range()[0])
1695 for (auto i = global_id; i < nrows; i += item_id.get_global_range()[0])
1696 {
1697 auto block_id = i/pack_size ;
1698 //auto local_id = i%ellpack_size ;
1699
1700 int begin = access_block_row_offset[block_id] ;
1701 int end = access_block_row_offset[block_id+1] ;
1702 ValueType value = 0. ;
1703 for(int k=begin;k<end;++k)
1704 {
1705 //auto k = block_row_offset+j*ellpack_size+local_id ;
1706 const int col = access_cols[k * pack_size + local_id];
1707 if(col>=0)
1708 lds_x[local_id] = access_x[col];
1709 item_id.barrier(sycl::access::fence_space::local_space);
1710 if(col>=0)
1711 value += access_values[k * pack_size + local_id] * lds_x[local_id] ;
1712 item_id.barrier(sycl::access::fence_space::local_space);
1713 }
1714 access_y[i] = value ;
1715 }
1716 });
1717 });
1718 // clang-format on
1719 }
1720 else
1721 {
1722 const bool done = dispatchN<2, 20>(this->m_N, [&](auto N) {
1723 multN<decltype(N)::value>(x, y, queue);
1724 });
1725 if (!done)
1726 {
1727 auto device = queue.get_device();
1728
1729 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
1730 // getting the maximum work group size per thread
1731 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
1732 // building the best number of global thread
1733 auto total_threads = num_groups * ellpack_size;
1734
1735 // clang-format off
1736 int N = this->m_N;
1737 //int NxN = this->m_NxN;
1738 std::size_t pack_size = ellpack_size;
1739 auto nrows = m_profile->getNRows();
1740 auto nnz = m_profile->getNnz();
1741
1742 auto internal_profile = m_profile->internal();
1743 auto& kcol = internal_profile->getKCol();
1744 auto& block_row_offset = internal_profile->getBlockRowOffset();
1745 auto& block_cols = internal_profile->getBlockCols();
1746
1747 queue.submit(
1748 [&](sycl::handler& cgh)
1749 {
1750 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
1751 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
1752 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
1753
1754
1755 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
1756 auto access_y = y.template get_access<sycl::access::mode::discard_write>(cgh);
1757
1758 //sycl::local_accessor<ValueT,1> local_mem(sycl::range<1>(ellpack_size*N), cgh);
1759 //sycl::nd_range<1> range{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}};
1760 sycl::range<1> range{total_threads} ;
1761 auto tile = Tile(N) ;
1762 //cgh.parallel_for<class compute_mult>(r, [&](sycl::nd_item<1> item_id)
1763 cgh.parallel_for<class compute_block_mult>(range,
1764 [=] (sycl::item<1> item_id /*sycl::nd_item<1> item_id*/)
1765 {
1766 auto id = item_id.get_id(0);
1767 //auto local_id = item_id.get_local_id(0);
1768 //auto block_id = item_id.get_group(0) ;
1769 //auto global_id = item_id.get_global_id(0);
1770
1771
1772
1773 for (auto i = id; i < nrows; i += item_id.get_range()[0] /*item_id.get_local_range(0)*/)
1774 {
1775 auto block_id = i/pack_size ;
1776 auto local_id = i%pack_size ;
1777
1778 std::size_t block_row_offset = access_block_row_offset[block_id] ;
1779 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
1780
1781 for(int ieq=0;ieq<N;++ieq)
1782 {
1783 ValueType value = 0. ;
1784 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
1785 {
1786 // value += access_values[k]* access_x[access_cols[k]] ;
1787 value += tile.mult(ieq,
1788 local_id,
1789 k,
1790 access_cols,
1791 access_values,
1792 access_x) ;
1793 }
1794 access_y[i*N+ieq] = value;
1795 }
1796 }
1797 });
1798 });
1799 }
1800
1801#ifdef OLD
1802 switch(this->m_N)
1803 {
1804 case 2:
1805 multN<2>(x,y,queue) ;
1806 break ;
1807 case 3:
1808 multN<3>(x,y,queue) ;
1809 break;
1810 case 4:
1811 multN<4>(x,y,queue) ;
1812 break;
1813 default:
1814 auto device = queue.get_device();
1815
1816 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
1817 // getting the maximum work group size per thread
1818 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
1819 // building the best number of global thread
1820 auto total_threads = num_groups * ellpack_size;
1821
1822 // clang-format off
1823 int N = this->m_N;
1824 //int NxN = this->m_NxN;
1825 std::size_t pack_size = ellpack_size;
1826 auto nrows = m_profile->getNRows();
1827 auto nnz = m_profile->getNnz();
1828
1829 auto internal_profile = m_profile->internal();
1830 auto& kcol = internal_profile->getKCol();
1831 auto& block_row_offset = internal_profile->getBlockRowOffset();
1832 auto& block_cols = internal_profile->getBlockCols();
1833
1834 queue.submit(
1835 [&](sycl::handler& cgh)
1836 {
1837 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
1838 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
1839 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
1840
1841
1842 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
1843 auto access_y = y.template get_access<sycl::access::mode::discard_write>(cgh);
1844
1845 //sycl::local_accessor<ValueT,1> local_mem(sycl::range<1>(ellpack_size*N), cgh);
1846 //sycl::nd_range<1> range{sycl::range<1>{total_threads},sycl::range<1>{ellpack_size}};
1847 sycl::range<1> range{total_threads} ;
1848 auto tile = Tile(N) ;
1849 //cgh.parallel_for<class compute_mult>(r, [&](sycl::nd_item<1> item_id)
1850 cgh.parallel_for<class compute_block_mult>(range,
1851 [=] (sycl::item<1> item_id /*sycl::nd_item<1> item_id*/)
1852 {
1853 auto id = item_id.get_id(0);
1854 //auto local_id = item_id.get_local_id(0);
1855 //auto block_id = item_id.get_group(0) ;
1856 //auto global_id = item_id.get_global_id(0);
1857
1858
1859
1860 for (auto i = id; i < nrows; i += item_id.get_range()[0] /*item_id.get_local_range(0)*/)
1861 {
1862 auto block_id = i/pack_size ;
1863 auto local_id = i%pack_size ;
1864
1865 std::size_t block_row_offset = access_block_row_offset[block_id] ;
1866 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
1867
1868 for(int ieq=0;ieq<N;++ieq)
1869 {
1870 ValueType value = 0. ;
1871 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
1872 {
1873 // value += access_values[k]* access_x[access_cols[k]] ;
1874 value += tile.mult(ieq,
1875 local_id,
1876 k,
1877 access_cols,
1878 access_values,
1879 access_x) ;
1880 }
1881 access_y[i*N+ieq] = value;
1882 }
1883 }
1884 });
1885 });
1886 }
1887#endif
1888 /*
1889 {
1890 int N = this->m_N;
1891 //int NxN = this->m_NxN;
1892 std::size_t pack_size = ellpack_size;
1893 auto nrows = m_profile->getNRows();
1894 auto nnz = m_profile->getNnz();
1895
1896 auto internal_profile = m_profile->internal();
1897 auto& kcol = internal_profile->getKCol();
1898 auto& block_row_offset = internal_profile->getBlockRowOffset();
1899 auto& block_cols = internal_profile->getBlockCols();
1900
1901 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> x_acc(x);
1902 for(int il=0;il<nrows;++il)
1903 {
1904 alien_info([&] { cout()<<"X["<<il<<"]:"; });
1905 for(int i=0;i<N;++i)
1906 alien_info([&] { cout()<<x_acc[il*N+i]<<","; });
1907
1908 }
1909 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> y_acc(y);
1910 for(int il=0;il<nrows;++il)
1911 {
1912 alien_info([&] { cout()<<"Y["<<il<<"]:";});
1913 for(int i=0;i<N;++i)
1914 alien_info([&] { cout()<<y_acc[il*N+i]<<","; });
1915 }
1916 Universe().traceMng()->flush() ;
1917
1918 auto tile = Tile(N) ;
1919 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> matrix_acc(m_values);
1920 sycl::host_accessor<int, 1, sycl::access::mode::read> cols_acc(block_cols);
1921 sycl::host_accessor<int, 1, sycl::access::mode::read> brow_offset_acc(block_row_offset);
1922 for(std::size_t irow=0;irow<nrows;++irow)
1923 {
1924 auto ib = irow/ellpack_size ;
1925 auto il = irow%ellpack_size ;
1926 for(int ieq=0;ieq<N;++ieq)
1927 {
1928 std::cout<<"LINE["<<il<<","<<ieq<<"] : ";
1929 ValueType value = 0. ;
1930 for(std::size_t k=brow_offset_acc[ib];k<brow_offset_acc[ib+1];++k)
1931 {
1932 // value += access_values[k]* access_x[access_cols[k]] ;
1933 value += tile.mult(ieq,
1934 il,
1935 k,
1936 cols_acc,
1937 matrix_acc,
1938 x_acc) ;
1939 }
1940 std::cout<<"\nY_CPU["<<irow<<","<<ieq<<"]:"<<value<<std::endl;
1941 }
1942 }
1943 }*/
1944 // clang-format on
1945 }
1946 }
1947
1948 template <typename ValueT, int EllPackSize>
1949 void MatrixInternal<ValueT, EllPackSize>::mult(ValueBufferType& x, ValueBufferType& y) const
1950 {
1951 this->mult(x, y, SYCLEnv::instance()->internal()->queue());
1952 }
1953
1954 template <typename ValueT, int EllPackSize>
1955 void MatrixInternal<ValueT, EllPackSize>::addExtMult(ValueBufferType& x,
1956 ValueBufferType& y,
1957 sycl::queue& queue) const
1958 {
1959 auto device = queue.get_device();
1960
1961 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
1962 // getting the maximum work group size per thread
1963 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
1964 // building the best number of global thread
1965 auto total_threads = num_groups * ellpack_size;
1966
1967 //int N = this->m_N;
1968 //int NxN = this->m_NxN;
1969 std::size_t pack_size = ellpack_size;
1970
1971 // clang-format off
1972 auto nrows = m_profile->getNRows();
1973 auto nnz = m_profile->getNnz();
1974
1975 auto interface_nrow = m_ext_profile->getNRows();
1976 auto ext_profile = m_ext_profile->internal();
1977 auto& kcol = ext_profile->getKCol();
1978 auto& block_row_offset = ext_profile->getBlockRowOffset();
1979 auto& block_cols = ext_profile->getBlockCols();
1980 // clang-format on
1981 if(this->m_N==1)
1982 {
1983 // clang-format off
1984 queue.submit([&](sycl::handler& cgh)
1985 {
1986 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
1987 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
1988 auto access_values = m_ext_values->template get_access<sycl::access::mode::read>(cgh);
1989
1990 auto access_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
1991
1992 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
1993 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
1994
1995 cgh.parallel_for<class compute_ext_mult>(
1996 sycl::range<1>{total_threads},
1997 [=] (sycl::item<1> item_id)
1998 {
1999 auto id = item_id.get_id(0);
2000
2001 for (auto i = id; i < interface_nrow; i += item_id.get_range()[0])
2002 {
2003 auto block_id = i / pack_size ;
2004 auto local_id = i % pack_size ;
2005
2006 int block_row_offset = access_block_row_offset[block_id]*ellpack_size ;
2007 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
2008
2009 ValueType value = 0. ;
2010 for(int j=0;j<block_row_size;++j)
2011 {
2012 auto k = block_row_offset+j*ellpack_size+local_id ;
2013 value += access_values[k]* access_x[access_cols[k]] ;
2014 }
2015 access_y[access_row_ids[i]] += value ;
2016 }
2017 });
2018 });
2019 // clang-format on
2020 }
2021 else
2022 {
2023 const bool done = dispatchN<2, 20>(
2024 this->m_N,
2025 [&](auto N)
2026 {
2027 addExtMultN<decltype(N)::value>(1., x, y, queue);
2028 });
2029 if (!done)
2030 {
2031 // clang-format off
2032
2033 int N = this->m_N;
2034 queue.submit([&](sycl::handler& cgh)
2035 {
2036 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2037 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2038 auto access_values = m_ext_values->template get_access<sycl::access::mode::read>(cgh);
2039
2040 auto access_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
2041
2042 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2043 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2044
2045 auto tile = Tile(N) ;
2046 cgh.parallel_for<class compute_ext_mult>(
2047 sycl::range<1>{total_threads},
2048 [=] (sycl::item<1> item_id)
2049 {
2050 auto id = item_id.get_id(0);
2051
2052 for (auto i = id; i < interface_nrow; i += item_id.get_range()[0])
2053 {
2054 auto block_id = i / pack_size ;
2055 auto local_id = i % pack_size ;
2056
2057 int block_row_offset = access_block_row_offset[block_id] * pack_size ;
2058 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
2059
2060 for(int ieq=0;ieq<N;++ieq)
2061 {
2062 ValueType value = 0. ;
2063 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
2064 {
2065 // value += access_values[k]* access_x[access_cols[k]] ;
2066 value += tile.mult(ieq,
2067 local_id,
2068 k,
2069 access_cols,
2070 access_values,
2071 access_x) ;
2072 }
2073 access_y[access_row_ids[i]*N+ieq] += value;
2074 }
2075 }
2076 });
2077 });
2078 // clang-format on
2079 }
2080 }
2081 }
2082
2083 template <typename ValueT, int EllPackSize>
2084 void MatrixInternal<ValueT, EllPackSize>::addExtMult(ValueBufferType& x,
2085 ValueBufferType& y) const
2086 {
2087 this->addExtMult(x, y, SYCLEnv::instance()->internal()->queue());
2088 }
2089
2090 template <typename ValueT, int EllPackSize>
2091 void MatrixInternal<ValueT, EllPackSize>::addLMult(ValueType alpha,
2092 ValueBufferType& x,
2093 ValueBufferType& y,
2094 sycl::queue& queue) const
2095 {
2096
2097 if(this->m_N==1)
2098 {
2099 auto device = queue.get_device();
2100
2101 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
2102 // getting the maximum work group size per thread
2103 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
2104
2105 std::size_t pack_size = ellpack_size;
2106 auto nrows = m_profile->getNRows();
2107 auto nnz = m_profile->getNnz();
2108
2109 auto internal_profile = m_profile->internal();
2110 auto& kcol = internal_profile->getKCol();
2111 auto& block_row_offset = internal_profile->getBlockRowOffset();
2112 auto& block_cols = internal_profile->getBlockCols();
2113
2114 auto& mask = internal_profile->getLowerMask();
2115
2116 // clang-format off
2117 queue.submit(
2118 [&](sycl::handler& cgh)
2119 {
2120 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2121 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2122 auto access_mask = mask.template get_access<sycl::access::mode::read>(cgh);
2123 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
2124
2125
2126 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2127 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2128
2129 auto blocks_needed = (nrows + ellpack_size - 1) / ellpack_size;
2130 auto blocks_target = std::max(blocks_needed, num_groups * 4UL);
2131 auto total_threads = blocks_target * ellpack_size;
2132
2133 sycl::local_accessor<ValueType, 1> lds_x{pack_size, cgh};
2134 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{pack_size}};
2135 //cgh.parallel_for<class compute_lmult>(sycl::range<1>{total_threads},[=] (sycl::item<1> item_id)
2136 cgh.parallel_for<class compute_lmult>(r,
2137 [=](sycl::nd_item<1> item_id)
2138 {
2139 //auto id = item_id.get_id(0);
2140 auto local_id = item_id.get_local_id(0);
2141 //auto block_id = item_id.get_group(0) ;
2142 auto global_id = item_id.get_global_id(0);
2143
2144 //for (auto i = id; i < nrows; i += item_id.get_range()[0])
2145 for (auto i = global_id; i < nrows; i += item_id.get_global_range()[0])
2146 {
2147 auto block_id = i/pack_size ;
2148 //auto local_id = i%ellpack_size ;
2149
2150 int begin = access_block_row_offset[block_id] ;
2151 int end = access_block_row_offset[block_id+1] ;
2152
2153 ValueType value = access_y[i] ;
2154 for(int k=begin;k<end;++k)
2155 {
2156 //auto k = block_row_offset+j*ellpack_size+local_id ;
2157 const int col = access_cols[k * pack_size + local_id];
2158 if(col>=0)
2159 lds_x[local_id] = access_x[col];
2160 item_id.barrier(sycl::access::fence_space::local_space);
2161 if(access_mask[k * pack_size + local_id])
2162 value += alpha * access_values[k * pack_size + local_id] * lds_x[local_id] ;
2163 item_id.barrier(sycl::access::fence_space::local_space);
2164 }
2165 access_y[i] = value ;
2166 }
2167 });
2168 });
2169 }
2170 else
2171 {
2172 const bool done = dispatchN<2, 20>(
2173 this->m_N,
2174 [&](auto N)
2175 {
2176 addLMultN<decltype(N)::value>(alpha,x, y, queue);
2177 });
2178
2179 if (!done)
2180 {
2181 auto device = queue.get_device();
2182
2183 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
2184 // getting the maximum work group size per thread
2185 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
2186 auto total_threads = num_groups * ellpack_size;
2187
2188 int N = this->m_N;
2189 int NxN = this->m_NxN;
2190 std::size_t pack_size = ellpack_size;
2191 auto nrows = m_profile->getNRows();
2192 auto nnz = m_profile->getNnz();
2193
2194 auto internal_profile = m_profile->internal();
2195 auto& kcol = internal_profile->getKCol();
2196 auto& block_row_offset = internal_profile->getBlockRowOffset();
2197 auto& block_cols = internal_profile->getBlockCols();
2198
2199 auto& mask = internal_profile->getLowerMask();
2200
2201 // clang-format off
2202 queue.submit(
2203 [&](sycl::handler& cgh)
2204 {
2205 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2206 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2207 auto access_mask = mask.template get_access<sycl::access::mode::read>(cgh);
2208 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
2209
2210 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2211 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2212 sycl::range<1> range{total_threads} ;
2213 auto tile = Tile(N) ;
2214 cgh.parallel_for<class compute_block_lmult>(
2215 sycl::range<1>{total_threads},
2216 [=] (sycl::item<1> item_id)
2217 {
2218 auto id = item_id.get_id(0);
2219 //auto local_id = item_id.get_local_id(0);
2220 //auto block_id = item_id.get_group(0) ;
2221 //auto global_id = item_id.get_global_id(0);
2222
2223 for (auto i = id; i < nrows; i += item_id.get_range()[0])
2224 {
2225 auto block_id = i/ellpack_size ;
2226 auto local_id = i%ellpack_size ;
2227
2228 int block_row_offset = access_block_row_offset[block_id] ;
2229 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
2230
2231 for(int ieq=0;ieq<N;++ieq)
2232 {
2233 ValueType value = 0. ;
2234 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
2235 {
2236 // value += access_values[k]* access_x[access_cols[k]] ;
2237 value += tile.mult(ieq,
2238 local_id,
2239 k,
2240 access_cols,
2241 access_mask,
2242 access_values,
2243 access_x) ;
2244 }
2245 access_y[i*N+ieq] += alpha*value;
2246 }
2247 }
2248 });
2249 });
2250 }
2251#ifdef OLD
2252 switch(this->m_N)
2253 {
2254 case 2:
2255 addLMultN<2>(alpha,x,y,queue);
2256 break;
2257 case 3:
2258 addLMultN<3>(alpha,x,y,queue);
2259 break;
2260 case 4:
2261 addLMultN<4>(alpha,x,y,queue);
2262 break;
2263 default:
2264 {
2265 auto device = queue.get_device();
2266
2267 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
2268 // getting the maximum work group size per thread
2269 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
2270 auto total_threads = num_groups * ellpack_size;
2271
2272 int N = this->m_N;
2273 int NxN = this->m_NxN;
2274 std::size_t pack_size = ellpack_size;
2275 auto nrows = m_profile->getNRows();
2276 auto nnz = m_profile->getNnz();
2277
2278 auto internal_profile = m_profile->internal();
2279 auto& kcol = internal_profile->getKCol();
2280 auto& block_row_offset = internal_profile->getBlockRowOffset();
2281 auto& block_cols = internal_profile->getBlockCols();
2282
2283 auto& mask = internal_profile->getLowerMask();
2284
2285 // clang-format off
2286 queue.submit(
2287 [&](sycl::handler& cgh)
2288 {
2289 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2290 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2291 auto access_mask = mask.template get_access<sycl::access::mode::read>(cgh);
2292 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
2293
2294 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2295 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2296 sycl::range<1> range{total_threads} ;
2297 auto tile = Tile(N) ;
2298 cgh.parallel_for<class compute_block_lmult>(
2299 sycl::range<1>{total_threads},
2300 [=] (sycl::item<1> item_id)
2301 {
2302 auto id = item_id.get_id(0);
2303 //auto local_id = item_id.get_local_id(0);
2304 //auto block_id = item_id.get_group(0) ;
2305 //auto global_id = item_id.get_global_id(0);
2306
2307 for (auto i = id; i < nrows; i += item_id.get_range()[0])
2308 {
2309 auto block_id = i/ellpack_size ;
2310 auto local_id = i%ellpack_size ;
2311
2312 int block_row_offset = access_block_row_offset[block_id] ;
2313 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
2314
2315 for(int ieq=0;ieq<N;++ieq)
2316 {
2317 ValueType value = 0. ;
2318 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
2319 {
2320 // value += access_values[k]* access_x[access_cols[k]] ;
2321 value += tile.mult(ieq,
2322 local_id,
2323 k,
2324 access_cols,
2325 access_mask,
2326 access_values,
2327 access_x) ;
2328 }
2329 access_y[i*N+ieq] += alpha*value;
2330 }
2331 }
2332 });
2333 });
2334 }
2335 break;
2336 }
2337#endif
2338 }
2339
2340#ifdef PRINT_DEBUG_INFO
2341 {
2342 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> h_acc(m_values);
2343 sycl::host_accessor<int, 1, sycl::access::mode::read> kcol_acc(block_row_offset);
2344 sycl::host_accessor<uint8_t, 1, sycl::access::mode::read> mask_acc(mask) ;
2345 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> x_acc(x);
2346 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> y_acc(y);
2347 for(int irow=0;irow<nrows;++irow)
2348 {
2349 std::cout<<"XK["<<irow<<"]=\n";
2350 for(int i=0;i<N;++i)
2351 std::cout<<x_acc[irow*N+i]<<"\n";
2352 }
2353 for(int irow=0;irow<nrows;++irow)
2354 {
2355 auto ib=irow/ellpack_size ;
2356 auto il=irow%ellpack_size ;
2357 std::cout<<"LSOLVE MAT["<<irow<<"]:"<<std::endl;
2358 for(int k=kcol_acc[ib];k<kcol_acc[ib+1];++k)
2359 {
2360 if(mask_acc[k*ellpack_size+il])
2361 {
2362 for(int i=0;i<N;++i)
2363 {
2364 for(int j=0;j<N;++j)
2365 std::cout<<h_acc[(k*NxN+i*N+j)*ellpack_size+il]<<" ";
2366 std::cout<<std::endl;
2367 }
2368 }
2369 }
2370 std::cout<<"Y ="<<std::endl;
2371 for(int i=0;i<N;++i)
2372 std::cout<<y_acc[irow*N+i]<<"\n";
2373 }
2374 }
2375#endif
2376 }
2377
2378 template <typename ValueT, int EllPackSize>
2379 void MatrixInternal<ValueT, EllPackSize>::addLMult(ValueType alpha, ValueBufferType& x, ValueBufferType& y) const
2380 {
2381 this->addLMult(alpha, x, y, SYCLEnv::instance()->internal()->queue());
2382 }
2383
2384 template <typename ValueT, int EllPackSize>
2385 void MatrixInternal<ValueT, EllPackSize>::addExtLMult(ValueType alpha,
2386 ValueBufferType& x,
2387 ValueBufferType& y,
2388 sycl::queue& queue) const
2389 {
2390 auto device = queue.get_device();
2391
2392 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
2393 // getting the maximum work group size per thread
2394 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
2395 // building the best number of global thread
2396 auto total_threads = num_groups * ellpack_size;
2397
2398 //int N = this->m_N;
2399 //int NxN = this->m_NxN;
2400 std::size_t pack_size = ellpack_size;
2401
2402 // clang-format off
2403 auto nrows = m_profile->getNRows();
2404 auto nnz = m_profile->getNnz();
2405
2406 auto interface_nrow = m_ext_profile->getNRows();
2407 auto first_upper_ghost = m_first_upper_ghost_index ;
2408 auto ext_profile = m_ext_profile->internal();
2409 auto& kcol = ext_profile->getKCol();
2410 auto& block_row_offset = ext_profile->getBlockRowOffset();
2411 auto& block_cols = ext_profile->getBlockCols();
2412
2413 if(this->m_N==1)
2414 {
2415 // clang-format off
2416 queue.submit([&](sycl::handler& cgh)
2417 {
2418 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2419 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2420 auto access_values = m_ext_values->template get_access<sycl::access::mode::read>(cgh);
2421
2422 auto access_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
2423
2424 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2425 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2426
2427 cgh.parallel_for<class compute_ext_lmult>(
2428 sycl::range<1>{total_threads},
2429 [=] (sycl::item<1> item_id)
2430 {
2431 auto id = item_id.get_id(0);
2432
2433 for (auto i = id; i < interface_nrow; i += item_id.get_range()[0])
2434 {
2435 auto block_id = i / pack_size ;
2436 auto local_id = i % pack_size ;
2437
2438 int block_row_offset = access_block_row_offset[block_id] * pack_size ;
2439 auto block_row_size = access_block_row_offset[block_id+1] - access_block_row_offset[block_id] ;
2440
2441 ValueType value = 0. ;
2442 for(int j=0;j<block_row_size;++j)
2443 {
2444 auto k = block_row_offset+j*ellpack_size+local_id ;
2445 if(access_cols[k]<first_upper_ghost)
2446 value += access_values[k]* access_x[access_cols[k]] ;
2447 }
2448 access_y[access_row_ids[i]] += alpha * value ;
2449 }
2450 });
2451 });
2452 // clang-format on
2453 }
2454 else
2455 {
2456 const bool done = dispatchN<2, 20>(
2457 this->m_N,
2458 [&](auto N)
2459 {
2460 addExtLMultN<decltype(N)::value>(alpha, x, y, queue);
2461 });
2462 if (!done)
2463 {
2464 // clang-format off
2465
2466 int N = this->m_N;
2467 queue.submit([&](sycl::handler& cgh)
2468 {
2469 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2470 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2471 auto access_values = m_ext_values->template get_access<sycl::access::mode::read>(cgh);
2472
2473 auto access_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
2474
2475 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2476 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2477
2478 auto tile = Tile(N) ;
2479 cgh.parallel_for<class compute_ext_lmult>(
2480 sycl::range<1>{total_threads},
2481 [=] (sycl::item<1> item_id)
2482 {
2483 auto id = item_id.get_id(0);
2484
2485 for (auto i = id; i < interface_nrow; i += item_id.get_range()[0])
2486 {
2487 auto block_id = i / pack_size ;
2488 auto local_id = i % pack_size ;
2489
2490 int block_row_offset = access_block_row_offset[block_id] * pack_size ;
2491 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
2492
2493 for(int ieq=0;ieq<N;++ieq)
2494 {
2495 ValueType value = 0. ;
2496 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
2497 {
2498 // value += access_values[k]* access_x[access_cols[k]] ;
2499 if(access_cols[k]<first_upper_ghost)
2500 value += tile.mult(ieq,
2501 local_id,
2502 k,
2503 access_cols,
2504 access_values,
2505 access_x) ;
2506 }
2507 access_y[access_row_ids[i]*N+ieq] += alpha * value;
2508 }
2509 }
2510 });
2511 });
2512 // clang-format on
2513 }
2514 }
2515 }
2516
2517 template <typename ValueT, int EllPackSize>
2518 void MatrixInternal<ValueT, EllPackSize>::addExtLMult(ValueType alpha,
2519 ValueBufferType& x,
2520 ValueBufferType& y) const
2521 {
2522 this->addExtLMult(alpha, x, y, SYCLEnv::instance()->internal()->queue());
2523 }
2524
2525 template <typename ValueT, int EllPackSize>
2526 void MatrixInternal<ValueT, EllPackSize>::addUMult(ValueType alpha,
2527 ValueBufferType& x,
2528 ValueBufferType& y,
2529 sycl::queue& queue) const
2530 {
2531 if(this->m_N==1)
2532 {
2533 auto device = queue.get_device();
2534
2535 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
2536 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
2537
2538 // clang-format off
2539 int N = this->m_N;
2540 int NxN = this->m_NxN;
2541 std::size_t pack_size = ellpack_size;
2542 auto nrows = m_profile->getNRows() ;
2543 auto nnz = m_profile->getNnz() ;
2544
2545 auto blocks_needed = (nrows + ellpack_size - 1) / ellpack_size;
2546 auto blocks_target = std::max(blocks_needed, num_groups * 4UL);
2547 auto total_threads = blocks_target * ellpack_size;
2548
2549 auto internal_profile = m_profile->internal() ;
2550 auto& kcol = internal_profile->getKCol() ;
2551 auto& block_row_offset = internal_profile->getBlockRowOffset() ;
2552 auto& block_cols = internal_profile->getBlockCols() ;
2553 auto& mask = internal_profile->getUpperMask() ;
2554
2555 // COMPUTE VALUES
2556 queue.submit(
2557 [&](sycl::handler& cgh)
2558 {
2559 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2560 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2561 auto access_mask = mask.template get_access<sycl::access::mode::read>(cgh);
2562 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
2563
2564
2565 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2566 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2567
2568 sycl::local_accessor<ValueType, 1> lds_x{pack_size, cgh};
2569 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{pack_size}};
2570 //cgh.parallel_for<class compute_umult>(sycl::range<1>{total_threads},[=] (sycl::item<1> item_id)
2571 cgh.parallel_for<class compute_umult>(r,
2572 [=](sycl::nd_item<1> item_id)
2573 {
2574 //auto id = item_id.get_id(0);
2575 auto local_id = item_id.get_local_id(0);
2576 //auto block_id = item_id.get_group(0) ;
2577 auto global_id = item_id.get_global_id(0);
2578
2579 for (auto i = global_id; i < nrows; i += item_id.get_global_range()[0])
2580 {
2581 auto block_id = i/pack_size ;
2582 //auto local_id = i%ellpack_size ;
2583
2584 auto begin = access_block_row_offset[block_id] ;
2585 auto end = access_block_row_offset[block_id+1] ;
2586 //auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
2587
2588 ValueType value = access_y[i] ;
2589 for(int k=begin;k<end;++k)
2590 {
2591 //auto k = block_row_offset+j*ellpack_size+local_id ;
2592 //value += alpha * access_mask[k] * access_values[k]* access_x[access_cols[k]] ;
2593 const int col = access_cols[k * pack_size + local_id];
2594 if(col>=0)
2595 lds_x[local_id] = access_x[col];
2596 item_id.barrier(sycl::access::fence_space::local_space);
2597 if(access_mask[k * pack_size + local_id])
2598 value += alpha * access_values[k * pack_size + local_id] * lds_x[local_id] ;
2599 item_id.barrier(sycl::access::fence_space::local_space);
2600 }
2601 access_y[i] = value ;
2602 }
2603 });
2604 });
2605 }
2606 else
2607 {
2608 const bool done = dispatchN<2, 20>(
2609 this->m_N,
2610 [&](auto N)
2611 {
2612 addUMultN<decltype(N)::value>(alpha, x, y, queue);
2613 });
2614
2615 if (!done)
2616 {
2617 auto device = queue.get_device();
2618
2619 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
2620 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
2621 auto total_threads = num_groups * ellpack_size;
2622
2623 // clang-format off
2624 int N = this->m_N;
2625 int NxN = this->m_NxN;
2626 std::size_t pack_size = ellpack_size;
2627 auto nrows = m_profile->getNRows() ;
2628 auto nnz = m_profile->getNnz() ;
2629
2630 auto internal_profile = m_profile->internal() ;
2631 auto& kcol = internal_profile->getKCol() ;
2632 auto& block_row_offset = internal_profile->getBlockRowOffset() ;
2633 auto& block_cols = internal_profile->getBlockCols() ;
2634 auto& mask = internal_profile->getUpperMask() ;
2635
2636 // COMPUTE VALUES
2637 queue.submit(
2638 [&](sycl::handler& cgh)
2639 {
2640 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2641 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2642 auto access_mask = mask.template get_access<sycl::access::mode::read>(cgh);
2643 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
2644
2645
2646 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2647 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2648
2649 auto tile = Tile(N) ;
2650 cgh.parallel_for<class compute_block_umult>(
2651 sycl::range<1>{total_threads},
2652 [=] (sycl::item<1> item_id)
2653 {
2654 auto id = item_id.get_id(0);
2655 //auto local_id = item_id.get_local_id(0);
2656 //auto block_id = item_id.get_group(0) ;
2657 //auto global_id = item_id.get_global_id(0);
2658
2659 for (auto i = id; i < nrows; i += item_id.get_range()[0])
2660 {
2661 auto block_id = i/ellpack_size ;
2662 auto local_id = i%ellpack_size ;
2663
2664 std::size_t block_row_offset = access_block_row_offset[block_id] ;
2665 std::size_t block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
2666 for(int ieq=0;ieq<N;++ieq)
2667 {
2668 ValueType value = 0. ;
2669 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
2670 {
2671 // value += access_values[k]* access_x[access_cols[k]] ;
2672 value += tile.mult(ieq,
2673 local_id,
2674 k,
2675 access_cols,
2676 access_mask,
2677 access_values,
2678 access_x) ;
2679 }
2680 access_y[i*N+ieq] += alpha*value;
2681 }
2682 }
2683 });
2684 });
2685 }
2686#ifdef OLD
2687 switch(this->m_N)
2688 {
2689 case 2:
2690 addUMultN<2>(alpha,x,y,queue);
2691 break;
2692 case 3:
2693 addUMultN<3>(alpha,x,y,queue);
2694 break;
2695 case 4:
2696 addUMultN<4>(alpha,x,y,queue);
2697 break;
2698 default:
2699 {
2700 auto device = queue.get_device();
2701
2702 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
2703 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
2704 auto total_threads = num_groups * ellpack_size;
2705
2706 // clang-format off
2707 int N = this->m_N;
2708 int NxN = this->m_NxN;
2709 std::size_t pack_size = ellpack_size;
2710 auto nrows = m_profile->getNRows() ;
2711 auto nnz = m_profile->getNnz() ;
2712
2713 auto internal_profile = m_profile->internal() ;
2714 auto& kcol = internal_profile->getKCol() ;
2715 auto& block_row_offset = internal_profile->getBlockRowOffset() ;
2716 auto& block_cols = internal_profile->getBlockCols() ;
2717 auto& mask = internal_profile->getUpperMask() ;
2718
2719 // COMPUTE VALUES
2720 queue.submit(
2721 [&](sycl::handler& cgh)
2722 {
2723 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2724 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2725 auto access_mask = mask.template get_access<sycl::access::mode::read>(cgh);
2726 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
2727
2728
2729 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2730 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2731
2732 auto tile = Tile(N) ;
2733 cgh.parallel_for<class compute_block_umult>(
2734 sycl::range<1>{total_threads},
2735 [=] (sycl::item<1> item_id)
2736 {
2737 auto id = item_id.get_id(0);
2738 //auto local_id = item_id.get_local_id(0);
2739 //auto block_id = item_id.get_group(0) ;
2740 //auto global_id = item_id.get_global_id(0);
2741
2742 for (auto i = id; i < nrows; i += item_id.get_range()[0])
2743 {
2744 auto block_id = i/ellpack_size ;
2745 auto local_id = i%ellpack_size ;
2746
2747 int block_row_offset = access_block_row_offset[block_id] ;
2748 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
2749 for(int ieq=0;ieq<N;++ieq)
2750 {
2751 ValueType value = 0. ;
2752 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
2753 {
2754 // value += access_values[k]* access_x[access_cols[k]] ;
2755 value += tile.mult(ieq,
2756 local_id,
2757 k,
2758 access_cols,
2759 access_mask,
2760 access_values,
2761 access_x) ;
2762 }
2763 access_y[i*N+ieq] += alpha*value;
2764 }
2765 }
2766 });
2767 });
2768 }
2769 break;
2770 }
2771#endif
2772 }
2773 // clang-format on
2774
2775#ifdef PRINT_DEBUG_INFO
2776 {
2777 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> h_acc(m_values);
2778 sycl::host_accessor<int, 1, sycl::access::mode::read> kcol_acc(block_row_offset);
2779 sycl::host_accessor<uint8_t, 1, sycl::access::mode::read> mask_acc(mask) ;
2780 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> x_acc(x);
2781 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> y_acc(y);
2782 for(int irow=0;irow<nrows;++irow)
2783 {
2784 std::cout<<"XK["<<irow<<"]=\n";
2785 for(int i=0;i<N;++i)
2786 std::cout<<x_acc[irow*N+i]<<"\n";
2787 }
2788 for(int irow=0;irow<nrows;++irow)
2789 {
2790 auto ib=irow/ellpack_size ;
2791 auto il=irow%ellpack_size ;
2792 std::cout<<"USOLVE MAT["<<irow<<"]:"<<std::endl;
2793 for(int k=kcol_acc[ib];k<kcol_acc[ib+1];++k)
2794 {
2795 if(mask_acc[k*ellpack_size+il])
2796 {
2797 for(int i=0;i<N;++i)
2798 {
2799 for(int j=0;j<N;++j)
2800 std::cout<<h_acc[(k*NxN+i*N+j)*ellpack_size+il]<<" ";
2801 std::cout<<std::endl;
2802 }
2803 }
2804 }
2805 std::cout<<"Y ="<<std::endl;
2806 for(int i=0;i<N;++i)
2807 std::cout<<y_acc[irow*N+i]<<"\n";
2808 }
2809 }
2810#endif
2811 }
2812
2813 template <typename ValueT, int EllPackSize>
2814 void MatrixInternal<ValueT, EllPackSize>::addUMult(ValueType alpha, ValueBufferType& x, ValueBufferType& y) const
2815 {
2816 this->addUMult(alpha, x, y, SYCLEnv::instance()->internal()->queue());
2817 }
2818
2819 template <typename ValueT, int EllPackSize>
2820 void MatrixInternal<ValueT, EllPackSize>::addExtUMult(ValueType alpha,
2821 ValueBufferType& x,
2822 ValueBufferType& y,
2823 sycl::queue& queue) const
2824 {
2825 auto device = queue.get_device();
2826
2827 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
2828 // getting the maximum work group size per thread
2829 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
2830 // building the best number of global thread
2831 auto total_threads = num_groups * ellpack_size;
2832
2833 //int N = this->m_N;
2834 //int NxN = this->m_NxN;
2835 std::size_t pack_size = ellpack_size;
2836
2837 // clang-format off
2838 auto nrows = m_profile->getNRows();
2839 auto nnz = m_profile->getNnz();
2840
2841 auto interface_nrow = m_ext_profile->getNRows();
2842 auto first_upper_ghost = m_first_upper_ghost_index ;
2843 auto ext_profile = m_ext_profile->internal();
2844 auto& kcol = ext_profile->getKCol();
2845 auto& block_row_offset = ext_profile->getBlockRowOffset();
2846 auto& block_cols = ext_profile->getBlockCols();
2847 // clang-format on
2848 if(this->m_N==1)
2849 {
2850 // clang-format off
2851 queue.submit([&](sycl::handler& cgh)
2852 {
2853 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2854 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2855 auto access_values = m_ext_values->template get_access<sycl::access::mode::read>(cgh);
2856
2857 auto access_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
2858
2859 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2860 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2861
2862 cgh.parallel_for<class compute_ext_umult>(
2863 sycl::range<1>{total_threads},
2864 [=] (sycl::item<1> item_id)
2865 {
2866 auto id = item_id.get_id(0);
2867
2868 for (auto i = id; i < interface_nrow; i += item_id.get_range()[0])
2869 {
2870 auto block_id = i / pack_size ;
2871 auto local_id = i % pack_size ;
2872
2873 int block_row_offset = access_block_row_offset[block_id] * pack_size ;
2874 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
2875
2876 ValueType value = 0. ;
2877 for(int j=0;j<block_row_size;++j)
2878 {
2879 auto k = block_row_offset+j*pack_size+local_id ;
2880 if(access_cols[k]>=first_upper_ghost)
2881 value += access_values[k]* access_x[access_cols[k]] ;
2882 }
2883 access_y[access_row_ids[i]] += alpha * value ;
2884 }
2885 });
2886 });
2887 // clang-format on
2888 }
2889 else
2890 {
2891 const bool done = dispatchN<2, 20>(
2892 this->m_N,
2893 [&](auto N)
2894 {
2895 addExtUMultN<decltype(N)::value>(alpha, x, y, queue);
2896 });
2897 if (!done)
2898 {
2899 // clang-format off
2900
2901 int N = this->m_N;
2902 queue.submit([&](sycl::handler& cgh)
2903 {
2904 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
2905 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
2906 auto access_values = m_ext_values->template get_access<sycl::access::mode::read>(cgh);
2907
2908 auto access_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
2909
2910 auto access_x = x.template get_access<sycl::access::mode::read>(cgh);
2911 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
2912
2913 auto tile = Tile(N) ;
2914 cgh.parallel_for<class compute_ext_umult>(
2915 sycl::range<1>{total_threads},
2916 [=] (sycl::item<1> item_id)
2917 {
2918 auto id = item_id.get_id(0);
2919
2920 for (auto i = id; i < interface_nrow; i += item_id.get_range()[0])
2921 {
2922 auto block_id = i / pack_size ;
2923 auto local_id = i % pack_size ;
2924
2925 int block_row_offset = access_block_row_offset[block_id] * pack_size ;
2926 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
2927
2928 for(int ieq=0;ieq<N;++ieq)
2929 {
2930 ValueType value = 0. ;
2931 for(std::size_t k=block_row_offset;k<block_row_offset+block_row_size;++k)
2932 {
2933 // value += access_values[k]* access_x[access_cols[k]] ;
2934 if(access_cols[k]>=first_upper_ghost)
2935 value += tile.mult(ieq,
2936 local_id,
2937 k,
2938 access_cols,
2939 access_values,
2940 access_x) ;
2941 }
2942 access_y[access_row_ids[i]*N+ieq] += alpha * value;
2943 }
2944 }
2945 });
2946 });
2947 // clang-format on
2948 }
2949 }
2950 }
2951
2952 template <typename ValueT, int EllPackSize>
2953 void MatrixInternal<ValueT, EllPackSize>::addExtUMult(ValueType alpha, ValueBufferType& x,
2954 ValueBufferType& y) const
2955 {
2956 this->addExtUMult(alpha, x, y, SYCLEnv::instance()->internal()->queue());
2957 }
2958
2959 template <typename ValueT, int EllPackSize>
2960 void MatrixInternal<ValueT, EllPackSize>::multDiag(ValueBufferType& x, ValueBufferType& y, sycl::queue& queue) const
2961 {
2962 }
2963
2964 template <typename ValueT, int EllPackSize>
2965 void MatrixInternal<ValueT, EllPackSize>::multDiag(ValueBufferType& x, ValueBufferType& y) const
2966 {
2967 this->multDiag(x, y, SYCLEnv::instance()->internal()->queue());
2968 }
2969
2970 template <typename ValueT, int EllPackSize>
2971 void MatrixInternal<ValueT, EllPackSize>::multDiag(ValueBufferType& y, sycl::queue& queue) const
2972 {
2973 throw NotImplementedException(
2974 A_FUNCINFO, "MatrixInternal<ValueT, EllPackSize>::multDiag not implemented");
2975 }
2976
2977 template <typename ValueT, int EllPackSize>
2978 void MatrixInternal<ValueT, EllPackSize>::multDiag(ValueBufferType& y) const
2979 {
2980 this->multDiag(y, SYCLEnv::instance()->internal()->queue());
2981 }
2982
2983 template <typename ValueT, int EllPackSize>
2984 void MatrixInternal<ValueT, EllPackSize>::computeDiag(ValueBufferType& y, sycl::queue& queue) const
2985 {
2986
2987 auto device = queue.get_device();
2988
2989 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
2990 // getting the maximum work group size per thread
2991 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
2992 // building the best number of global thread
2993 auto total_threads = num_groups * ellpack_size;
2994
2995 // clang-format off
2996 int N = this->m_N;
2997 int NxN = this->m_NxN;
2998 int pack_size = ellpack_size ;
2999 auto nrows = m_profile->getNRows() ;
3000 auto nnz = m_profile->getNnz() ;
3001
3002 auto internal_profile = m_profile->internal() ;
3003 auto& kcol = internal_profile->getKCol() ;
3004 auto& block_row_offset = internal_profile->getBlockRowOffset() ;
3005 auto& block_cols = internal_profile->getBlockCols() ;
3006
3007 // COMPUTE VALUES
3008 queue.submit([&](sycl::handler& cgh)
3009 {
3010 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
3011 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
3012 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
3013 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
3014 if(N==1)
3015 {
3016 cgh.parallel_for<class compute_diag>(
3017 sycl::range<1>{total_threads},
3018 [=] (sycl::item<1> item_id)
3019 {
3020 auto id = item_id.get_id(0);
3021 //auto local_id = item_id.get_local_id(0);
3022 //auto block_id = item_id.get_group(0) ;
3023 //auto global_id = item_id.get_global_id(0);
3024
3025 for (auto i = id; i < nrows; i += item_id.get_range()[0])
3026 {
3027 auto block_id = i/pack_size ;
3028 auto local_id = i%pack_size ;
3029
3030 auto begin = access_block_row_offset[block_id] ;
3031 auto end = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
3032 for(int k=begin;k<end;++k)
3033 {
3034 auto kcol = k*pack_size+local_id ;
3035 if((access_cols[kcol])==int(i))
3036 access_y[i] = access_values[kcol] ;
3037 }
3038 }
3039 });
3040 }
3041 else
3042 {
3043 cgh.parallel_for<class compute_diagN>(
3044 sycl::range<1>{total_threads},
3045 [=] (sycl::item<1> item_id)
3046 {
3047 auto id = item_id.get_id(0);
3048 //auto local_id = item_id.get_local_id(0);
3049 //auto block_id = item_id.get_group(0) ;
3050 //auto global_id = item_id.get_global_id(0);
3051
3052 for (std::size_t i = id; i < nrows; i += item_id.get_range()[0])
3053 {
3054 auto block_id = i/pack_size ;
3055 auto local_id = i%pack_size ;
3056
3057 int block_row_offset = access_block_row_offset[block_id] ;
3058 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
3059 for(int k=0;k<block_row_size;++k)
3060 {
3061 auto kcol = (block_row_offset+k)*pack_size+local_id ;
3062 std::size_t col = access_cols[kcol] ;
3063 if(col==i)
3064 {
3065 for(int ieq=0;ieq<N;++ieq)
3066 {
3067 auto kval = (block_row_offset+k)*NxN + ieq*N + ieq ;
3068 access_y[i*N+ieq] = access_values[kval*pack_size+local_id] ;
3069 }
3070 }
3071 }
3072 }
3073 });
3074 }
3075 });
3076 /*
3077 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> h_acc(y);
3078 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> mat_acc(m_values);
3079 sycl::host_accessor<int, 1, sycl::access::mode::read> kcol_acc(block_row_offset);
3080 sycl::host_accessor<int, 1, sycl::access::mode::read> cols_acc(block_cols);
3081 for(int irow=0;irow<nrows;++irow)
3082 {
3083 auto ib=irow/ellpack_size ;
3084 auto il=irow%ellpack_size ;
3085 std::cout<<"MAT["<<irow<<","<<ib<<"]:";
3086 for(int k=kcol_acc[ib];k<kcol_acc[ib+1];++k)
3087 {
3088 std::cout<<"("<<cols_acc[k*ellpack_size+il]<<",";
3089 for(int i=0;i<N;++i)
3090 for(int j=0;j<N;++j)
3091 std::cout<<mat_acc[(k*NxN+i*N+j)*ellpack_size+il]<<",";
3092 std::cout<<")"<<std::endl;
3093 }
3094 }
3095 for(int irow=0;irow<nrows;++irow)
3096 {
3097 auto ib=irow/ellpack_size ;
3098 auto il=irow%ellpack_size ;
3099 std::cout<<"DIAG["<<irow<<","<<il<<"]:";
3100 for(int i=0;i<N;++i)
3101 std::cout<<h_acc[il*N+i]<<",";
3102 std::cout<<std::endl;
3103 }*/
3104 // clang-format on
3105 }
3106
3107 template <typename ValueT, int EllPackSize>
3108 void MatrixInternal<ValueT, EllPackSize>::computeBlockDiag(ValueBufferType& y, sycl::queue& queue) const
3109 {
3110
3111 auto device = queue.get_device();
3112
3113 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
3114 // getting the maximum work group size per thread
3115 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
3116 // building the best number of global thread
3117 auto total_threads = num_groups * ellpack_size;
3118
3119 // clang-format off
3120 int N = this->m_N;
3121 int NxN = this->m_NxN;
3122 int pack_size = ellpack_size ;
3123 auto nrows = m_profile->getNRows() ;
3124 auto nnz = m_profile->getNnz() ;
3125
3126 assert(y.size()>=NxN*nrows) ;
3127
3128 auto internal_profile = m_profile->internal() ;
3129 auto& kcol = internal_profile->getKCol() ;
3130 auto& block_row_offset = internal_profile->getBlockRowOffset() ;
3131 auto& block_cols = internal_profile->getBlockCols() ;
3132 {
3133 ValueBufferType buffer (sycl::range<1>(nrows*NxN)) ;
3134 // COMPUTE VALUES
3135 queue.submit([&](sycl::handler& cgh)
3136 {
3137 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
3138 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
3139 auto matrix = m_values.template get_access<sycl::access::mode::read>(cgh);
3140 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
3141 auto A = buffer.template get_access<sycl::access::mode::read_write>(cgh);
3142
3143 auto tile = Tile(N) ;
3144 cgh.parallel_for<class compute_block_diag>(
3145 sycl::range<1>{total_threads},
3146 [=] (sycl::item<1> item_id)
3147 {
3148 auto id = item_id.get_id(0);
3149 //auto local_id = item_id.get_local_id(0);
3150 //auto block_id = item_id.get_group(0) ;
3151 //auto global_id = item_id.get_global_id(0);
3152
3153 for (std::size_t global_id = id; global_id < nrows; global_id += item_id.get_range()[0])
3154 {
3155 auto block_id = global_id/pack_size ;
3156 auto local_id = global_id%pack_size ;
3157
3158 int block_row_offset = access_block_row_offset[block_id] ;
3159 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
3160 for(int k=0;k<block_row_size;++k)
3161 {
3162 auto kcol = (block_row_offset+k)*ellpack_size+local_id ;
3163 std::size_t col = access_cols[kcol] ;
3164 if(col==global_id)
3165 {
3166 // Copy Diag Matrix in A
3167 for(int i=0;i<N;++i)
3168 for(int j=0;j<N;++j)
3169 A[tile._ij(local_id,i,j)] = matrix[tile._ijk(k,i,j)+local_id] ;
3170 }
3171 }
3172 }
3173 });
3174 });
3175 }
3176 // clang-format on
3177 }
3178
3179 template <typename ValueT, int EllPackSize>
3180 void MatrixInternal<ValueT, EllPackSize>::computeDiag(ValueBufferType& y) const
3181 {
3182 this->computeDiag(y, SYCLEnv::instance()->internal()->queue());
3183 }
3184
3185 template <typename ValueT, int EllPackSize>
3186 void MatrixInternal<ValueT, EllPackSize>::computeBlockDiag(ValueBufferType& y) const
3187 {
3188 this->computeBlockDiag(y, SYCLEnv::instance()->internal()->queue());
3189 }
3190
3191
3192 template <typename ValueT, int EllPackSize>
3193 void MatrixInternal<ValueT, EllPackSize>::multInvDiag(ValueBufferType& y, sycl::queue& queue) const
3194 {
3195 throw NotImplementedException(
3196 A_FUNCINFO, "MatrixInternal<ValueT, EllPackSize>::multInvDiag not implemented");
3197 }
3198
3199 template <typename ValueT, int EllPackSize>
3200 void MatrixInternal<ValueT, EllPackSize>::multInvDiag(ValueBufferType& y) const
3201 {
3202 this->multInvDiag(y, SYCLEnv::instance()->internal()->queue());
3203 }
3204
3205 template <typename ValueT, int EllPackSize>
3206 void MatrixInternal<ValueT, EllPackSize>::computeInvDiag(ValueBufferType& y, sycl::queue& queue) const
3207 {
3208
3209 auto device = queue.get_device();
3210
3211 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
3212 // getting the maximum work group size per thread
3213 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
3214
3215 // clang-format off
3216 int N = this->m_N;
3217 int NxN = this->m_NxN;
3218 std::size_t pack_size = ellpack_size ;
3219 auto nrows = m_profile->getNRows() ;
3220 auto nnz = m_profile->getNnz() ;
3221
3222 auto internal_profile = m_profile->internal() ;
3223 auto& kcol = internal_profile->getKCol() ;
3224 auto& block_row_offset = internal_profile->getBlockRowOffset() ;
3225 auto& block_cols = internal_profile->getBlockCols() ;
3226
3227 auto blocks_needed = (nrows + ellpack_size - 1) / ellpack_size;
3228 auto blocks_target = std::max(blocks_needed, num_groups * 4UL);
3229 auto total_threads = blocks_target * ellpack_size;
3230
3231 // COMPUTE VALUES
3232 queue.submit([&](sycl::handler& cgh)
3233 {
3234 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
3235 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
3236 auto access_values = m_values.template get_access<sycl::access::mode::read>(cgh);
3237 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
3238 if(N==1)
3239 {
3240 //sycl::local_accessor<ValueType, 1> lds_x{pack_size, cgh};
3241 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{pack_size}};
3242
3243 cgh.parallel_for<class compute_inv_diag>(r,
3244 [=](sycl::nd_item<1> item_id)
3245 {
3246 auto local_id = item_id.get_local_id(0);
3247 auto global_id = item_id.get_global_id(0);
3248 for (auto i = global_id; i < nrows; i += item_id.get_global_range()[0])
3249 {
3250 auto block_id = i/pack_size ;
3251 //auto local_id = i%ellpack_size ;
3252
3253 int begin = access_block_row_offset[block_id];
3254 int end = access_block_row_offset[block_id+1];
3255 for(int k=begin;k<end;++k)
3256 {
3257 auto kcol = k*pack_size+local_id ;
3258 if((access_cols[kcol])==int(i) && (access_values[kcol]!=0) )
3259 access_y[i] = 1./access_values[kcol] ;
3260 }
3261 }
3262 });
3263
3264 /*
3265 cgh.parallel_for<class compute_inv_diag>(
3266 sycl::range<1>{total_threads},
3267 [=] (sycl::item<1> item_id)
3268 {
3269 auto id = item_id.get_id(0);
3270 //auto local_id = item_id.get_local_id(0);
3271 //auto block_id = item_id.get_group(0) ;
3272 //auto global_id = item_id.get_global_id(0);
3273
3274 for (auto i = id; i < nrows; i += item_id.get_range()[0])
3275 {
3276 auto block_id = i/ellpack_size ;
3277 auto local_id = i%ellpack_size ;*/
3278
3279 }
3280 else
3281 {
3282 //sycl::local_accessor<ValueType, 1> lds_x{pack_size, cgh};
3283 sycl::nd_range<1> r{sycl::range<1>{total_threads},sycl::range<1>{pack_size}};
3284 cgh.parallel_for<class compute_inv_diagN>(r,
3285 [=](sycl::nd_item<1> item_id)
3286 {
3287 auto local_id = item_id.get_local_id(0);
3288 auto global_id = item_id.get_global_id(0);
3289 for (auto i = global_id; i < nrows; i += item_id.get_global_range()[0])
3290 {
3291 auto block_id = i/pack_size ;
3292
3293 int begin = access_block_row_offset[block_id];
3294 int end = access_block_row_offset[block_id+1];
3295
3296 for(int k=begin;k<end;++k)
3297 {
3298 auto kcol = k*pack_size+local_id ;
3299 std::size_t col = access_cols[kcol] ;
3300 if(col==i)
3301 {
3302 for(int ieq=0;ieq<N;++ieq)
3303 {
3304 auto kval = k*NxN + ieq*N + ieq ;
3305 auto diag = access_values[kval*pack_size+local_id] ;
3306 if(diag!=0)
3307 access_y[i*N+ieq] = 1./diag;
3308 else
3309 access_y[i*N+ieq] = 1.;
3310 }
3311 }
3312 }
3313 }
3314 });
3315
3316 /*
3317 cgh.parallel_for<class compute_inv_diagN>(
3318 sycl::range<1>{total_threads},
3319 [=] (sycl::item<1> item_id)
3320 {
3321 auto id = item_id.get_id(0);
3322 //auto local_id = item_id.get_local_id(0);
3323 //auto block_id = item_id.get_group(0) ;
3324 //auto global_id = item_id.get_global_id(0);
3325
3326 for (std::size_t i = id; i < nrows; i += item_id.get_range()[0])
3327 {
3328 auto block_id = i/pack_size ;
3329 auto local_id = i%pack_size ;
3330*/
3331 }
3332 });
3333 /*
3334 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> h_acc(y);
3335 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> mat_acc(m_values);
3336 sycl::host_accessor<int, 1, sycl::access::mode::read> kcol_acc(block_row_offset);
3337 sycl::host_accessor<int, 1, sycl::access::mode::read> cols_acc(block_cols);
3338 for(int irow=0;irow<nrows;++irow)
3339 {
3340 auto ib=irow/ellpack_size ;
3341 auto il=irow%ellpack_size ;
3342 std::cout<<"MAT["<<irow<<","<<ib<<"]:";
3343 for(int k=kcol_acc[ib];k<kcol_acc[ib+1];++k)
3344 {
3345 std::cout<<"("<<cols_acc[k*ellpack_size+il]<<",";
3346 for(int i=0;i<N;++i)
3347 for(int j=0;j<N;++j)
3348 std::cout<<mat_acc[(k*NxN+i*N+j)*ellpack_size+il]<<",";
3349 std::cout<<")"<<std::endl;
3350 }
3351 }
3352 for(int irow=0;irow<nrows;++irow)
3353 {
3354 auto ib=irow/ellpack_size ;
3355 auto il=irow%ellpack_size ;
3356 std::cout<<"DIAG["<<irow<<","<<il<<"]:";
3357 for(int i=0;i<N;++i)
3358 std::cout<<h_acc[il*N+i]<<",";
3359 std::cout<<std::endl;
3360 }*/
3361 // clang-format on
3362 }
3363
3364 template <typename ValueT, int EllPackSize>
3365 void MatrixInternal<ValueT, EllPackSize>::computeInvBlockDiag(ValueBufferType& y, sycl::queue& queue) const
3366 {
3367
3368 auto device = queue.get_device();
3369
3370 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
3371 // getting the maximum work group size per thread
3372 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
3373 // building the best number of global thread
3374 auto total_threads = num_groups * ellpack_size;
3375
3376 // clang-format off
3377 int N = this->m_N;
3378 int NxN = this->m_NxN;
3379 int pack_size = ellpack_size ;
3380 auto nrows = m_profile->getNRows() ;
3381 auto nnz = m_profile->getNnz() ;
3382
3383 assert(y.size()>=NxN*nrows) ;
3384
3385 auto internal_profile = m_profile->internal() ;
3386 auto& kcol = internal_profile->getKCol() ;
3387 auto& block_row_offset = internal_profile->getBlockRowOffset() ;
3388 auto& block_cols = internal_profile->getBlockCols() ;
3389 {
3390 ValueBufferType buffer (sycl::range<1>(internal_profile->m_block_nrows*NxN*pack_size)) ;
3391 // COMPUTE VALUES
3392 queue.submit([&](sycl::handler& cgh)
3393 {
3394 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
3395 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
3396 auto matrix = m_values.template get_access<sycl::access::mode::read>(cgh);
3397 auto access_y = y.template get_access<sycl::access::mode::read_write>(cgh);
3398 auto A = buffer.template get_access<sycl::access::mode::read_write>(cgh);
3399
3400 auto lu = LU<decltype(matrix),decltype(access_y),decltype(A)>{N,matrix} ;
3401 //sycl::local_accessor<ValueT,1> local_mem(sycl::range<1>(pack_size*NxN), cgh);
3402
3403 cgh.parallel_for<class compute_inv_block_diag>(
3404 sycl::range<1>{total_threads},
3405 [=] (sycl::item<1> item_id)
3406 {
3407 auto id = item_id.get_id(0);
3408 //auto local_id = item_id.get_local_id(0);
3409 //auto block_id = item_id.get_group(0) ;
3410 //auto global_id = item_id.get_global_id(0);
3411
3412 for (std::size_t global_id = id; global_id < nrows; global_id += item_id.get_range()[0])
3413 {
3414 auto block_id = global_id/pack_size ;
3415 auto local_id = global_id%pack_size ;
3416
3417 int begin = access_block_row_offset[block_id];
3418 int end = access_block_row_offset[block_id+1];
3419 for(int k=begin;k<end;++k)
3420 {
3421 auto kcol = k*pack_size+local_id ;
3422 std::size_t col = access_cols[kcol] ;
3423 if(col==global_id)
3424 {
3425 lu.factorize(global_id,
3426 local_id,
3427 block_id,
3428 k,
3429 A);
3430
3431 lu.inverse(global_id,
3432 local_id,
3433 block_id,
3434 A,
3435 access_y);
3436 /*
3437 {
3438 // Copy Diag Matrix in A
3439 for(int i=0;i<N;++i)
3440 for(int j=0;j<N;++j)
3441 A[tile._ij(local_id,i,j)] = matrix[tile._ijk(k,i,j)+local_id] ;
3442
3443 //Factorize A = LU
3444 for (int k = 0; k < N; ++k)
3445 {
3446 assert(A[tile._ij(local_id,k,k)] != 0);
3447 A[tile._ij(global_id,k,k)] = 1 / A[tile._ij(global_id,k,k)];
3448 for (int i = k + 1; i < N; ++i) {
3449 A[tile._ij(global_id,i,k)] *= A[tile._ij(global_id,k,k)];
3450 }
3451 for (int i = k + 1; i < N; ++i) {
3452 for (int j = k + 1; j < N; ++j) {
3453 A[tile._ij(global_id,i,j)] -= A[tile._ij(global_id,i,k)] * A[tile._ij(global_id,k,j)];
3454 }
3455 }
3456 }
3457
3458 // SET Y to Id
3459 for(int i=0;i<N;++i)
3460 for(int j=0;j<N;++j)
3461 access_y[tile._ij(global_id,i,j)] = 0. ;
3462 for(int i=0;i<m_N;++i)
3463 access_y[tile._ij(global_id,i,i)] = 1. ;
3464
3465 // L solve
3466 for (int i = 1; i < N; ++i)
3467 {
3468 for (int j = 0; j < i; ++j)
3469 {
3470 for(int k=0;k<N;++k)
3471 access_y[tile._ij(global_id,i,k)] -= A[tile._ij(global_id,i,j)] * access_y[tile._ij(global_id,j,k)];
3472 }
3473 }
3474
3475 // U solve
3476 for (int i = N - 1; i >= 0; --i)
3477 {
3478 for (int j = N - 1; j > i; --j)
3479 {
3480 for(int k=0;k<N;++k)
3481 access_y[tile._ij(global_id,i,k)] -= A[tile._ij(global_id,i,j)] * access_y[tile._ij(global_id,j,k)];
3482 }
3483 for(int k=0;k<N;++k)
3484 access_y[tile._ij(global_id,i,k)] *= A[tile._ij(global_id,i,i)];
3485 }
3486 }*/
3487 }
3488 }
3489 }
3490 });
3491 });
3492 }
3493#ifdef PRINT_DEBUG_INFO
3494 {
3495 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> h_acc(y);
3496 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> mat_acc(m_values);
3497 sycl::host_accessor<int, 1, sycl::access::mode::read> kcol_acc(block_row_offset);
3498 sycl::host_accessor<int, 1, sycl::access::mode::read> cols_acc(block_cols);
3499 std::vector<ValueT> test_y(nrows*NxN) ;
3500 std::vector<ValueT> test_A(internal_profile->m_block_nrows*NxN*pack_size) ;
3501 auto lu = LU<decltype(mat_acc),decltype(test_y.data()),decltype(test_A.data())>{N,mat_acc} ;
3502 for(int irow=0;irow<nrows;++irow)
3503 {
3504 auto ib=irow/ellpack_size ;
3505 auto il=irow%ellpack_size ;
3506 alien_info([&] {
3507 cout() <<"MAT["<<irow<<","<<ib<<"]:";
3508 });
3509 for(int k=kcol_acc[ib];k<kcol_acc[ib+1];++k)
3510 {
3511 alien_info([&] {
3512 cout() <<"\tCOLS["<<k<<"] = "<<cols_acc[k*pack_size+il];
3513 });
3514 /*
3515 for(int i=0;i<N;++i)
3516 {
3517 for(int j=0;j<N;++j)
3518 std::cout<<mat_acc[(k*NxN+i*N+j)*pack_size+il]<<" ";
3519 std::cout<<std::endl;
3520 }*/
3521 if(cols_acc[k*pack_size+il]==irow)
3522 {
3523
3524 alien_info([&] {
3525 cout() <<"\tCOMPUTE INVERSE["<<cols_acc[k*pack_size+il]<<"]" ;
3526 for(int i=0;i<N;++i)
3527 {
3528 for(int j=0;j<N;++j)
3529 cout()<<mat_acc[(k*NxN+i*N+j)*pack_size+il]<<" ";
3530 }
3531 });
3532 lu.factorize(irow,
3533 il,
3534 ib,
3535 k,
3536 test_A.data());
3537
3538 lu.inverse(irow,
3539 il,
3540 ib,
3541 test_A.data(),
3542 test_y.data());
3543
3544 alien_info([&] {
3545 for(int i=0;i<N;++i)
3546 {
3547 for(int j=0;j<N;++j)
3548 std::cout<<test_y[irow*NxN+i*N+j]<<",";
3549 std::cout<<std::endl;
3550 }
3551 });
3552 }
3553 }
3554 }
3555 for(int irow=0;irow<nrows;++irow)
3556 {
3557 auto ib=irow/ellpack_size ;
3558 auto il=irow%ellpack_size ;
3559 alien_info([&] {
3560 cout()<<"INV BLOCK DIAG["<<irow<<","<<il<<"]:\n";
3561 for(int i=0;i<N;++i)
3562 {
3563 for(int j=0;j<N;++j)
3564 cout()<<h_acc[irow*NxN+i*N+j]<<",";
3565 //std::cout<<std::endl;
3566 }
3567 });
3568 }
3569 }
3570#endif
3571 // clang-format on
3572 }
3573
3574 template <typename ValueT, int EllPackSize>
3575 void MatrixInternal<ValueT, EllPackSize>::computeInvDiag(ValueBufferType& y) const
3576 {
3577 this->computeInvDiag(y, SYCLEnv::instance()->internal()->queue());
3578 }
3579
3580 template <typename ValueT, int EllPackSize>
3581 void MatrixInternal<ValueT, EllPackSize>::computeInvBlockDiag(ValueBufferType& y) const
3582 {
3583 this->computeInvBlockDiag(y, SYCLEnv::instance()->internal()->queue());
3584 }
3585
3586
3587 template <typename ValueT, int EllPackSize>
3588 void MatrixInternal<ValueT, EllPackSize>::scal(ValueBufferType& y, sycl::queue& queue)
3589 {
3590
3591 auto device = queue.get_device();
3592
3593 auto num_groups = queue.get_device().get_info<sycl::info::device::max_compute_units>();
3594 // getting the maximum work group size per thread
3595 auto max_work_group_size = queue.get_device().get_info<sycl::info::device::max_work_group_size>();
3596 // building the best number of global thread
3597 auto total_threads = num_groups * ellpack_size;
3598
3599 // clang-format off
3600 auto nrows = m_profile->getNRows() ;
3601 auto nnz = m_profile->getNnz() ;
3602
3603 auto internal_profile = m_profile->internal() ;
3604 auto& kcol = internal_profile->getKCol() ;
3605 auto& block_row_offset = internal_profile->getBlockRowOffset() ;
3606 auto& block_cols = internal_profile->getBlockCols() ;
3607 {
3608 // COMPUTE VALUES
3609 queue.submit([&](sycl::handler& cgh)
3610 {
3611 auto access_block_row_offset = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
3612 auto access_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
3613 auto access_values = m_values.template get_access<sycl::access::mode::read_write>(cgh);
3614 auto access_y = y.template get_access<sycl::access::mode::read>(cgh);
3615
3616
3617 cgh.parallel_for<class scal_matrix>(sycl::range<1>{total_threads},
3618 [=] (sycl::item<1> item_id)
3619 {
3620 auto id = item_id.get_id(0);
3621 //auto local_id = item_id.get_local_id(0);
3622 //auto block_id = item_id.get_group(0) ;
3623 //auto global_id = item_id.get_global_id(0);
3624
3625 for (auto i = id; i < nrows; i += item_id.get_range()[0])
3626 {
3627 auto block_id = i/ellpack_size ;
3628 auto local_id = i%ellpack_size ;
3629
3630 int block_row_offset = access_block_row_offset[block_id]*ellpack_size ;
3631 auto block_row_size = access_block_row_offset[block_id+1]-access_block_row_offset[block_id] ;
3632 for(int j=0;j<block_row_size;++j)
3633 {
3634 auto k = block_row_offset+j*ellpack_size+local_id ;
3635 access_values[k] *= access_y[i] ;
3636 }
3637 }
3638 });
3639 });
3640 }
3641
3642 if(m_ext_profile)
3643 {
3644
3645 auto interface_nrows = m_ext_profile->getNRows();
3646 auto ext_nnz = m_ext_profile->getNnz();
3647 auto ext_block_nnz = m_ext_profile->getBlockNnz();
3648
3649 auto ext_internal_profile = m_ext_profile->internal();
3650 auto& ext_block_row_offset = ext_internal_profile->getBlockRowOffset();
3651
3652 {
3653 queue.submit([&](sycl::handler& cgh)
3654 {
3655 auto access_block_row_offset = ext_block_row_offset.template get_access<sycl::access::mode::read>(cgh);
3656 auto access_block_values = m_ext_values->template get_access<sycl::access::mode::read_write>(cgh);
3657 auto access_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
3658 auto access_y = y.template get_access<sycl::access::mode::read>(cgh);
3659
3660 cgh.parallel_for<class set_matrix_values7>(sycl::range<1>{total_threads},
3661 [=] (sycl::item<1> item_id)
3662 {
3663 auto id = item_id.get_id(0);
3664
3665 for (auto i = id; i < interface_nrows; i += item_id.get_range()[0])
3666 {
3667 auto block_id = i/ellpack_size ;
3668 auto local_id = i%ellpack_size ;
3669 auto begin = access_block_row_offset[block_id] ;
3670 auto end = access_block_row_offset[block_id+1] ;
3671
3672 for(int k=begin;k<end;++k)
3673 {
3674 access_block_values[k*ellpack_size+local_id] *= access_y[access_row_ids[i]] ;
3675 }
3676 }
3677 });
3678 });
3679 }
3680 }
3681 }
3682
3683 template <typename ValueT, int EllPackSize>
3684 void MatrixInternal<ValueT, EllPackSize>::scal(ValueBufferType& y)
3685 {
3686 this->scal(y, SYCLEnv::instance()->internal()->queue());
3687 }
3688
3689 template <typename ValueT, int EllPackSize>
3690 void MatrixInternal<ValueT, EllPackSize>::
3691 copyDevicePointers(int local_offset,
3692 std::size_t nrows,
3693 std::size_t nnz,
3694 int* rows,
3695 int* ncols,
3696 int* cols,
3697 ValueT* values) const
3698 {
3699#ifdef PRINT_DEBUG_INFO
3700 alien_info([&] {
3701 cout()<<"SYCLMATRIX COPY DEVICE POINTERS : "<<nrows<<" "<<nnz;
3702 });
3703#endif
3704
3705 auto env = SYCLEnv::instance() ;
3706 auto max_num_threads = env->maxNumThreads() ;
3707 auto& queue = env->internal()->queue() ;
3708 auto pack_size = ellpack_size;
3709
3710 auto internal_profile = m_profile->internal() ;
3711 auto& kcol = internal_profile->getKCol() ;
3712 auto& block_row_offset = internal_profile->getBlockRowOffset() ;
3713 auto& block_cols = internal_profile->getBlockCols() ;
3714
3715 auto local_row_size = m_profile->localRowSize();
3716 /*
3717 {
3718 std::cout<<"CHECK KCOL : "<<std::endl ;
3719 sycl::host_accessor<int, 1, sycl::access::mode::read> h_kcol(kcol);
3720 sycl::host_accessor<int, 1, sycl::access::mode::read> h_bkcol(block_row_offset);
3721 sycl::host_accessor<int, 1, sycl::access::mode::read> h_bcols(block_cols);
3722 sycl::host_accessor<ValueT, 1, sycl::access::mode::read> h_values(m_values);
3723 for(int irow=0;irow<nrows;++irow)
3724 {
3725 auto block_id = irow/pack_size ;
3726 auto local_id = irow%pack_size ;
3727 std::cout<<"KCOL["<<irow<<"] = "<<h_kcol[irow]<<" ROW SIZE = "<<h_kcol[irow+1]-h_kcol[irow]<<std::endl ;
3728
3729 int begin = h_bkcol[block_id] ;
3730 int end = h_bkcol[block_id+1] ;
3731 int row_size = h_kcol[irow+1] - h_kcol[irow];
3732
3733 std::cout<<"ELLPACK MAT["<<irow<<","<<block_id<<","<<local_id<<"]:";
3734 for(int k=begin;k<end;++k)
3735 {
3736 std::cout<<h_bcols[k*pack_size+local_id]<<" "<<h_values[k*pack_size+local_id]<<" ";
3737 }
3738 std::cout<<std::endl ;
3739 }
3740 }*/
3741 if(local_row_size == nullptr)
3742 {
3743
3744 queue.submit( [&](sycl::handler& cgh)
3745 {
3746 auto acc_val = m_values.template get_access<sycl::access::mode::read>(cgh);
3747 auto acc_bkcol = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
3748 auto acc_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
3749 auto acc_kcol = kcol.template get_access<sycl::access::mode::read>(cgh);
3750 cgh.parallel_for<class copy_device_ptr>(
3751 sycl::range<1>{max_num_threads},
3752 [=] (sycl::item<1> itemId)
3753 {
3754 auto id = itemId.get_id(0);
3755 for (auto i = id; i < nrows; i += itemId.get_range()[0])
3756 {
3757 auto block_id = i/pack_size ;
3758 auto local_id = i%pack_size ;
3759
3760 int begin = acc_bkcol[block_id] ;
3761 int end = acc_bkcol[block_id+1] ;
3762 int row_size = acc_kcol[i+1] - acc_kcol[i];
3763
3764 rows[i] = (int) (local_offset + i) ;
3765 ncols[i] = row_size ;
3766
3767 auto offset = acc_kcol[i];
3768
3769 for(int k=begin;k<end;++k)
3770 {
3771 auto j = k-begin ;
3772 if(j<row_size)
3773 {
3774 values[offset+j] = acc_val[k*pack_size+local_id] ;
3775 cols[offset+j] = acc_cols[k*pack_size+local_id] ;
3776 }
3777 }
3778 }
3779 });
3780 }).wait();
3781#ifdef PRINT_DEBUG_INFO
3782 {
3783 std::vector<ValueT> h_data(nnz) ;
3784 queue.memcpy(h_data.data(), values, nnz * sizeof(ValueT)).wait();
3785 std::vector<int> h_cols(nnz) ;
3786 queue.memcpy(h_cols.data(), cols, nnz * sizeof(int)).wait();
3787 sycl::host_accessor<int, 1, sycl::access::mode::read> h_kcol(kcol);
3788 alien_info([&]{
3789 cout()<<"SYCL DEVICE COPY : "<<nrows<<" "<<nnz;
3790 for(int irow=0;irow<nrows;++irow)
3791 {
3792 cout()<<"COPY MAT["<<irow<<"]:"<<h_kcol[irow+1]-h_kcol[irow];
3793 for(int k=h_kcol[irow];k<h_kcol[irow+1];++k)
3794 {
3795 std::cout<<h_cols[k]<<" "<<h_data[k]<<" ";
3796 }
3797 std::cout<<std::endl ;
3798 }
3799 });
3800 }
3801#endif
3802 }
3803 else
3804 {
3805
3806 IndexBufferType lrowsize(local_row_size, sycl::range<1>(nrows));
3807 queue.submit( [&](sycl::handler& cgh)
3808 {
3809 auto acc_val = m_values.template get_access<sycl::access::mode::read>(cgh);
3810 auto acc_bkcol = block_row_offset.template get_access<sycl::access::mode::read>(cgh);
3811 auto acc_cols = block_cols.template get_access<sycl::access::mode::read>(cgh);
3812 auto acc_kcol = kcol.template get_access<sycl::access::mode::read>(cgh);
3813 auto acc_lrowsize = lrowsize.template get_access<sycl::access::mode::read>(cgh);
3814 cgh.parallel_for<class copy_device_ptr>(
3815 sycl::range<1>{max_num_threads},
3816 [=] (sycl::item<1> itemId)
3817 {
3818 auto id = itemId.get_id(0);
3819 for (auto i = id; i < nrows; i += itemId.get_range()[0])
3820 {
3821 auto block_id = i/pack_size ;
3822 auto local_id = i%pack_size ;
3823
3824 int begin = acc_bkcol[block_id] ;
3825 int end = acc_bkcol[block_id+1] ;
3826 int row_size = acc_kcol[i+1] - acc_kcol[i];
3827 int lrow_size = acc_lrowsize[i] ;
3828 rows[i] = (int) (local_offset + i) ;
3829 ncols[i] = row_size ;
3830
3831 auto offset = acc_kcol[i];
3832
3833 for(int k=begin;k<end;++k)
3834 {
3835 auto j = k-begin ;
3836 if(j<lrow_size)
3837 {
3838 values[offset+j] = acc_val[k*pack_size+local_id];
3839 cols[offset+j] = local_offset + acc_cols[k*pack_size+local_id];
3840 }
3841 }
3842 }
3843 });
3844 });
3845
3846 assert(m_ext_profile) ;
3847
3848 auto interface_nrows = m_ext_profile->getNRows();
3849 auto ext_nnz = m_ext_profile->getNnz();
3850 auto ext_block_nnz = m_ext_profile->getBlockNnz();
3851
3852 auto ext_internal_profile = m_ext_profile->internal();
3853 auto& ext_block_row_offset = ext_internal_profile->getBlockRowOffset();
3854 auto& ext_block_cols = ext_internal_profile->getBlockCols() ;
3855
3856 queue.submit([&](sycl::handler& cgh)
3857 {
3858 auto acc_kcol = kcol.template get_access<sycl::access::mode::read>(cgh);
3859 auto acc_ext_bkcol = ext_block_row_offset.template get_access<sycl::access::mode::read>(cgh);
3860 auto acc_ext_cols = ext_block_cols.template get_access<sycl::access::mode::read>(cgh);
3861 auto acc_ext_values = m_ext_values->template get_access<sycl::access::mode::read_write>(cgh);
3862 auto acc_row_ids = m_interface_row_ids->template get_access<sycl::access::mode::read>(cgh);
3863 auto acc_lrowsize = lrowsize.template get_access<sycl::access::mode::read>(cgh);
3864 auto acc_recv_uids = m_recv_uids->template get_access<sycl::access::mode::read>(cgh);
3865
3866 cgh.parallel_for<class ext_copy_device_ptr>(
3867 sycl::range<1>{max_num_threads},
3868 [=] (sycl::item<1> item_id)
3869 {
3870 auto id = item_id.get_id(0);
3871
3872 for (auto i = id; i < interface_nrows; i += item_id.get_range()[0])
3873 {
3874 auto block_id = i/pack_size ;
3875 auto local_id = i%pack_size ;
3876 auto begin = acc_ext_bkcol[block_id] ;
3877 auto end = acc_ext_bkcol[block_id+1] ;
3878
3879 auto row_id = acc_row_ids[i];
3880 int row_size = acc_kcol[row_id+1] - acc_kcol[row_id];
3881 int lrow_size = acc_lrowsize[row_id] ;
3882 auto offset = acc_kcol[row_id] + lrow_size;
3883
3884 for(int k=begin;k<end;++k)
3885 {
3886 auto j = k-begin ;
3887 if(j<row_size-lrow_size)
3888 {
3889 values[offset+j] = acc_ext_values[k*pack_size+local_id] ;
3890 cols[offset+j] = acc_recv_uids[acc_ext_cols[k*pack_size+local_id]] ;
3891 }
3892 }
3893 }
3894 });
3895 });
3896#ifdef PRINT_DEBUG_INFO
3897 {
3898 std::vector<ValueT> h_data(nnz) ;
3899 queue.memcpy(h_data.data(), values, nnz * sizeof(ValueT)).wait();
3900 std::vector<int> h_cols(nnz) ;
3901 queue.memcpy(h_cols.data(), cols, nnz * sizeof(int)).wait();
3902 sycl::host_accessor<int, 1, sycl::access::mode::read> h_kcol(kcol);
3903 alien_info([&] {
3904 cout()<<"SYCL DEVICE COPY : "<<nrows<<" "<<nnz<<" "<<ext_nnz<<" "<<interface_nrows;
3905 for(int irow=0;irow<nrows;++irow)
3906 {
3907 std::cout<<"COPY MAT["<<irow<<"]:"<<local_row_size[irow]<<" "<<h_kcol[irow+1]-h_kcol[irow];
3908 for(int k=h_kcol[irow];k<h_kcol[irow+1];++k)
3909 {
3910 std::cout<<h_cols[k]<<" "<<h_data[k]<<" ";
3911 }
3912 std::cout<<std::endl ;
3913 }
3914 });
3915 }
3916#endif
3917 }
3918 }
3919} // namespace SYCLInternal
3920
3921
3922
3923
3924template <typename ValueT>
3925SYCLBEllPackMatrix<ValueT>::SYCLBEllPackMatrix()
3926: IMatrixImpl(nullptr, AlgebraTraits<BackEnd::tag::sycl>::name())
3927, m_send_policy(SimpleCSRInternal::CommProperty::ASynch)
3928, m_recv_policy(SimpleCSRInternal::CommProperty::ASynch)
3929{}
3930
3931template <typename ValueT>
3932SYCLBEllPackMatrix<ValueT>::SYCLBEllPackMatrix(const MultiMatrixImpl* multi_impl)
3933: IMatrixImpl(multi_impl, AlgebraTraits<BackEnd::tag::sycl>::name())
3934, m_send_policy(SimpleCSRInternal::CommProperty::ASynch)
3935, m_recv_policy(SimpleCSRInternal::CommProperty::ASynch)
3936{}
3937
3938template <typename ValueT>
3939SYCLBEllPackMatrix<ValueT>::~SYCLBEllPackMatrix()
3940{
3941#ifdef ALIEN_USE_PERF_TIMER
3942 m_timer.printInfo("SYCLBELLPACK-MATRIX");
3943#endif
3944}
3945
3946template <typename ValueT>
3947bool SYCLBEllPackMatrix<ValueT>::
3948initMatrix(Arccore::MessagePassing::IMessagePassingMng* parallel_mng,
3949 Integer local_offset,
3950 Integer global_size,
3951 std::size_t nrows,
3952 int const* kcol,
3953 int const* cols,
3954 SimpleCSRInternal::DistStructInfo const& matrix_dist_info,
3955 int block_size)
3956{
3957#ifdef PRINT_DEBUG_INFO
3958 alien_debug([&] {
3959 cout()<<"SYCLBellPackMatrix::initMatrix : "<<local_offset<<" "<<nrows<<" "<<global_size;
3960 cout()<<"ParalleMng : "<<parallel_mng;
3961 });
3962#endif
3963
3964 m_nproc = 1;
3965 m_myrank = 0;
3966 m_parallel_mng = parallel_mng;
3967 if (parallel_mng) {
3968 m_nproc = m_parallel_mng->commSize();
3969 m_myrank = m_parallel_mng->commRank();
3970 }
3971 m_is_parallel = (m_nproc > 1);
3972
3973
3974 m_matrix_dist_info.copy(matrix_dist_info);
3975
3976 m_ellpack_size = PKSIZE;
3977 if (m_nproc > 1)
3978 {
3979 m_local_size = nrows;
3980 m_local_offset = local_offset;
3981 m_global_size = global_size;
3982
3983 m_ghost_size = m_matrix_dist_info.m_ghost_nrow;
3984
3985 auto& local_row_size = m_matrix_dist_info.m_local_row_size;
3986 auto sorted_cols = m_matrix_dist_info.m_cols.data();
3987
3988 std::size_t block_nrows = ProfileInternalType::nbBlocks(nrows);
3989
3990 ProfileInternalType::computeBlockRowOffset(m_block_row_offset, nrows, kcol);
3991
3992 // clang-format off
3993 alien_debug([&] {
3994 cout() << "OFFSET = "<<m_local_offset ;
3995 cout() << "NROWS = "<<nrows;
3996 cout() << "NNZ = "<<kcol[nrows];
3997 cout() << "BNROWS = "<<block_nrows;
3998 cout() << "BNNZ = "<<m_block_row_offset[block_nrows];
3999 });
4000 // clang-format on
4001 //Universe().traceMng()->flush();
4002 m_profilePKSIZE.reset( new ProfileInternalType{ nrows,
4003 kcol,
4004 sorted_cols,
4005 m_block_row_offset.data(),
4006 local_row_size.data() });
4007
4008 m_matrixPKSIZE.reset(new MatrixInternalType{ m_profilePKSIZE.get(), block_size });
4009
4010 // EXTRACT EXTERNAL PROFILE
4011 std::size_t interface_nrows = m_matrix_dist_info.m_interface_nrow;
4012 std::vector<int> ext_kcol(interface_nrows + 1);
4013
4014 int ext_nnz = 0;
4015 for (std::size_t i = 0; i < interface_nrows; ++i) {
4016 auto id = m_matrix_dist_info.m_interface_rows[i];
4017 ext_kcol[i] = ext_nnz;
4018 ext_nnz += kcol[id + 1] - kcol[id] - local_row_size[id];
4019 }
4020 ext_kcol[interface_nrows] = ext_nnz;
4021
4022 std::vector<int> ext_cols(ext_nnz);
4023 {
4024 int jcol = 0;
4025 for (std::size_t i = 0; i < interface_nrows; ++i) {
4026 auto id = m_matrix_dist_info.m_interface_rows[i];
4027 for (int k = kcol[id] + local_row_size[id]; k < kcol[id + 1]; ++k) {
4028 ext_cols[jcol++] = (int)(sorted_cols[k] - nrows);
4029 }
4030 }
4031 }
4032
4033 std::size_t ext_block_nrows = ProfileInternalType::nbBlocks(interface_nrows);
4034 ProfileInternalType::computeBlockRowOffset(m_ext_block_row_offset, interface_nrows, ext_kcol.data());
4035 // clang-format off
4036 alien_debug([&] {
4037 cout() << "EXT NROWS = "<<interface_nrows;
4038 cout() << "EXT NNZ = "<<ext_kcol[interface_nrows];
4039 cout() << "EXT BNROWS = "<<ext_block_nrows;
4040 cout() << "EXT BNNZ = "<<m_block_row_offset[ext_block_nrows];
4041 });
4042 // clang-format on
4043
4044 m_ext_profilePKSIZE.reset( new ProfileInternalType{ interface_nrows,
4045 ext_kcol.data(),
4046 ext_cols.data(),
4047 m_ext_block_row_offset.data(),
4048 nullptr });
4049
4050 m_matrixPKSIZE->m_ext_profile = m_ext_profilePKSIZE.get();
4051 typedef typename MatrixInternalType::IndexBufferType IndexBufferType;
4052 m_matrixPKSIZE->m_h_interface_row_ids = m_matrix_dist_info.m_interface_rows.data();
4053 m_matrixPKSIZE->m_interface_row_ids.reset(new IndexBufferType{ m_matrix_dist_info.m_interface_rows.data(),
4054 m_matrix_dist_info.m_interface_rows.size() });
4055
4056 m_matrixPKSIZE->m_send_ids.reset(new IndexBufferType{ m_matrix_dist_info.m_send_info.m_ids.data(),
4057 m_matrix_dist_info.m_send_info.m_ids.size() });
4058 m_matrixPKSIZE->m_recv_ids.reset(new IndexBufferType{ m_matrix_dist_info.m_recv_info.m_ids.data(),
4059 m_matrix_dist_info.m_recv_info.m_ids.size() });
4060 m_matrixPKSIZE->m_recv_uids.reset(new IndexBufferType{ m_matrix_dist_info.m_recv_info.m_uids.data(),
4061 m_matrix_dist_info.m_recv_info.m_uids.size() });
4062 }
4063 else
4064 {
4065 // clang-format off
4066 m_local_offset = 0;
4067 m_local_size = nrows;
4068 m_global_size = m_local_size;
4069 m_ghost_size = 0;
4070 // clang-format on
4071
4072 std::size_t block_nrows = ProfileInternalType::nbBlocks(nrows);
4073
4074 ProfileInternalType::computeBlockRowOffset(m_block_row_offset, nrows, kcol);
4075
4076 // clang-format off
4077 alien_debug([&] {
4078 cout() << "NROWS = "<<nrows;
4079 cout() << "NNZ = "<<kcol[nrows];
4080 cout() << "BNROWS = "<<block_nrows;
4081 cout() << "BNNZ = "<<m_block_row_offset[block_nrows];
4082 });
4083 // clang-format on
4084
4085 m_profilePKSIZE.reset( new ProfileInternalType{ nrows,
4086 kcol,
4087 cols,
4088 m_block_row_offset.data(),
4089 nullptr });
4090
4091 m_matrixPKSIZE.reset( new MatrixInternalType{ m_profilePKSIZE.get(), block_size });
4092 m_matrixPKSIZE->m_first_upper_ghost_index = m_matrix_dist_info.m_first_upper_ghost_index - m_local_size;
4093 }
4094
4095 return true;
4096}
4097
4098template <typename ValueT>
4099SYCLBEllPackMatrix<ValueT>*
4100SYCLBEllPackMatrix<ValueT>::cloneTo(const MultiMatrixImpl* multi) const
4101{
4102 auto block_size = blockSize() ;
4103 SYCLBEllPackMatrix<ValueT>* matrix = new SYCLBEllPackMatrix<ValueT>(multi);
4104 matrix->setBlockSize(block_size);
4105 matrix->initMatrix(m_parallel_mng,
4106 m_local_offset,
4107 m_global_size,
4108 m_profilePKSIZE->getNRows(),
4109 m_profilePKSIZE->kcol(),
4110 m_profilePKSIZE->cols(),
4111 m_matrix_dist_info,
4112 block_size);
4113 matrix->setMatrixValues(getAddressData(), true);
4114 return matrix;
4115}
4116
4117template <typename ValueT>
4118bool SYCLBEllPackMatrix<ValueT>::setMatrixValues(Arccore::Real const* values, bool only_host)
4119{
4120 return m_matrixPKSIZE->setMatrixValues(values, only_host);
4121}
4122
4123template <typename ValueT>
4124void SYCLBEllPackMatrix<ValueT>::notifyChanges()
4125{
4126 if (m_matrixPKSIZE.get())
4127 m_matrixPKSIZE->notifyChanges();
4128}
4129
4130template <typename ValueT>
4131void SYCLBEllPackMatrix<ValueT>::endUpdate()
4132{
4133 if (m_matrixPKSIZE.get() && m_matrixPKSIZE->needUpdate()) {
4134 m_matrixPKSIZE->endUpdate();
4135 this->updateTimestamp();
4136 }
4137}
4138
4139template <typename ValueT>
4140ValueT const* SYCLBEllPackMatrix<ValueT>::getAddressData() const
4141{
4142 return m_matrixPKSIZE->getHCsrData();
4143}
4144
4145template <typename ValueT>
4146ValueT const* SYCLBEllPackMatrix<ValueT>::data() const
4147{
4148 return m_matrixPKSIZE->getHCsrData();
4149}
4150
4151template <typename ValueT>
4152ValueT* SYCLBEllPackMatrix<ValueT>::getAddressData()
4153{
4154 return m_matrixPKSIZE->getHCsrData();
4155}
4156
4157template <typename ValueT>
4158ValueT* SYCLBEllPackMatrix<ValueT>::data()
4159{
4160 return m_matrixPKSIZE->getHCsrData();
4161}
4162
4163template <typename ValueT>
4164void SYCLBEllPackMatrix<ValueT>::mult(SYCLVector<ValueT> const& x, SYCLVector<ValueT>& y) const
4165{
4166 return m_matrixPKSIZE->mult(x.internal()->values(), y.internal()->values());
4167}
4168
4169template <typename ValueT>
4170void SYCLBEllPackMatrix<ValueT>::endDistMult(SYCLVector<ValueT> const& x, SYCLVector<ValueT>& y) const
4171{
4172 m_matrixPKSIZE->addExtMult(x.internal()->ghostValues(getGhostSize()*blockSize()), y.internal()->values());
4173}
4174
4175template <typename ValueT>
4176void SYCLBEllPackMatrix<ValueT>::addLMult(ValueT alpha, SYCLVector<ValueT> const& x, SYCLVector<ValueT>& y) const
4177{
4178 m_profilePKSIZE->dcol();
4179 return m_matrixPKSIZE->addLMult(alpha, x.internal()->values(), y.internal()->values());
4180}
4181
4182template <typename ValueT>
4183void SYCLBEllPackMatrix<ValueT>::endDistLMult(ValueT alpha, SYCLVector<ValueT> const& x, SYCLVector<ValueT>& y) const
4184{
4185 m_matrixPKSIZE->addExtLMult(alpha, x.internal()->ghostValues(getGhostSize()*blockSize()), y.internal()->values());
4186}
4187
4188template <typename ValueT>
4189void SYCLBEllPackMatrix<ValueT>::addUMult(ValueT alpha, SYCLVector<ValueT> const& x, SYCLVector<ValueT>& y) const
4190{
4191 m_profilePKSIZE->dcol();
4192 return m_matrixPKSIZE->addUMult(alpha, x.internal()->values(), y.internal()->values());
4193}
4194
4195template <typename ValueT>
4196void SYCLBEllPackMatrix<ValueT>::endDistUMult(ValueT alpha, SYCLVector<ValueT> const& x, SYCLVector<ValueT>& y) const
4197{
4198 m_matrixPKSIZE->addExtUMult(alpha, x.internal()->ghostValues(getGhostSize()*blockSize()), y.internal()->values());
4199}
4200
4201template <typename ValueT>
4202void SYCLBEllPackMatrix<ValueT>::multDiag(SYCLVector<ValueType> const& x, SYCLVector<ValueType>& y) const
4203{
4204 return m_matrixPKSIZE->multDiag(x.internal()->values(), y.internal()->values());
4205}
4206
4207template <typename ValueT>
4208void SYCLBEllPackMatrix<ValueT>::multDiag(SYCLVector<ValueType>& y) const
4209{
4210 return m_matrixPKSIZE->multDiag(y.internal()->values());
4211}
4212
4213template <typename ValueT>
4214void SYCLBEllPackMatrix<ValueT>::computeDiag(SYCLVector<ValueType>& y) const
4215{
4216 auto block_size = blockSize() ;
4217 if((y.blockSize()==1)||(y.blockSize()==block_size))
4218 return m_matrixPKSIZE->computeDiag(y.internal()->values());
4219 else if(y.blockSize()==block_size*block_size)
4220 return m_matrixPKSIZE->computeBlockDiag(y.internal()->values());
4221 else
4222 throw Arccore::FatalErrorException(A_FUNCINFO, "Matrix and Vector Block Size are not compatibility");
4223}
4224
4225template <typename ValueT>
4226void SYCLBEllPackMatrix<ValueT>::multInvDiag(SYCLVector<ValueType>& y) const
4227{
4228 return m_matrixPKSIZE->multInvDiag(y.internal()->values());
4229}
4230
4231template <typename ValueT>
4232void SYCLBEllPackMatrix<ValueT>::computeInvDiag(SYCLVector<ValueType>& y) const
4233{
4234 auto block_size = blockSize() ;
4235 if((y.blockSize()==1)||(y.blockSize()==block_size))
4236 return m_matrixPKSIZE->computeInvDiag(y.internal()->values());
4237 else if(y.blockSize()==block_size*block_size)
4238 return m_matrixPKSIZE->computeInvBlockDiag(y.internal()->values());
4239 else
4240 throw Arccore::FatalErrorException(A_FUNCINFO, "Matrix and Vector Block Size are not compatibility");
4241}
4242
4243template <typename ValueT>
4244void SYCLBEllPackMatrix<ValueT>::scal(SYCLVector<ValueType> const& diag)
4245{
4246 return m_matrixPKSIZE->scal(diag.internal()->values());
4247}
4248
4249
4250template <typename ValueT>
4251void SYCLBEllPackMatrix<ValueT>::copy(SYCLBEllPackMatrix<ValueT> const& matrix)
4252{
4253 auto block_size = blockSize() ;
4254 initMatrix(matrix.m_parallel_mng,
4255 matrix.m_local_offset,
4256 matrix.m_global_size,
4257 matrix.m_profilePKSIZE->getNRows(),
4258 matrix.m_profilePKSIZE->kcol(),
4259 matrix.m_profilePKSIZE->cols(),
4260 matrix.m_matrix_dist_info,
4261 block_size);
4262 if(block_size==matrix.blockSize())
4263 {
4264 if(m_is_parallel)
4265 m_matrixPKSIZE->setMatrixValues(matrix.m_matrixPKSIZE->m_values,
4266 (*matrix.m_matrixPKSIZE->m_ext_values));
4267 else
4268 m_matrixPKSIZE->setMatrixValues(matrix.m_matrixPKSIZE->m_values);
4269 }
4270 else
4271 {
4272 assert(m_profilePKSIZE.get());
4273 assert(m_matrixPKSIZE.get());
4274 auto nb_blocks = m_profilePKSIZE->getBlockNnz();
4275 if(m_is_parallel)
4276 m_matrixPKSIZE->copy(nb_blocks,
4277 block_size,
4278 matrix.m_matrixPKSIZE->m_values,
4279 (*matrix.m_matrixPKSIZE->m_ext_values),
4280 matrix.blockSize()) ;
4281 else
4282 m_matrixPKSIZE->copy(nb_blocks,
4283 block_size,
4284 matrix.m_matrixPKSIZE->m_values,
4285 matrix.blockSize()) ;
4286 }
4287}
4288
4289template <typename ValueT>
4290void SYCLBEllPackMatrix<ValueT>::
4291allocateDevicePointers(std::size_t nrows,
4292 std::size_t nnz,
4293 int** rows,
4294 int** ncols,
4295 int** cols,
4296 ValueT** values) const
4297{
4298 auto env = SYCLEnv::instance() ;
4299 auto max_num_treads = env->maxNumThreads() ;
4300 auto& queue = env->internal()->queue() ;
4301 auto rows_ptr = malloc_device<IndexType>(nrows, queue);
4302 auto ncols_ptr = malloc_device<IndexType>(nrows, queue);
4303 auto cols_ptr = malloc_device<IndexType>(nnz, queue);
4304 auto values_ptr = malloc_device<ValueT>(nnz, queue);
4305
4306 *rows = rows_ptr ;
4307 *ncols = ncols_ptr ;
4308 *cols = cols_ptr ;
4309 *values = values_ptr ;
4310}
4311
4312template <typename ValueT>
4313void SYCLBEllPackMatrix<ValueT>::
4314freeDevicePointers(int* rows, int* ncols, int* cols, ValueT* values) const
4315{
4316 auto env = SYCLEnv::instance() ;
4317 auto& queue = env->internal()->queue() ;
4318 sycl::free(rows,queue) ;
4319 sycl::free(ncols,queue) ;
4320 sycl::free(cols,queue) ;
4321 sycl::free(values,queue) ;
4322}
4323
4324
4325template <typename ValueT>
4326void SYCLBEllPackMatrix<ValueT>::
4327copyDevicePointers(std::size_t nrows,
4328 std::size_t nnz,
4329 int* rows,
4330 int* ncols,
4331 int* cols,
4332 ValueT* values) const
4333{
4334 m_matrixPKSIZE->copyDevicePointers(m_local_offset,nrows,nnz,rows,ncols,cols,values) ;
4335}
4336
4337
4338template <typename ValueT>
4339SYCLBEllPackMatrix<ValueT>::HCSRView
4340SYCLBEllPackMatrix<ValueT>::hcsrView(BackEnd::Memory::eType memory, int nrows, int nnz) const
4341{
4342 return HCSRView(this,memory, nrows, nnz) ;
4343}
4344/*---------------------------------------------------------------------------*/
4345
4346template class ALIEN_EXPORT SYCLBEllPackMatrix<double>;
4347template class ALIEN_EXPORT BEllPackStructInfo<PKSIZE, Integer>;
4348template class ALIEN_EXPORT SYCLInternal::MatrixInternal<double,PKSIZE>;
4349
4350
4351//template bool Alien::SYCLInternal::MatrixInternal<double,1014>::setMatrixValues(Alien::SYCLInternal::MatrixInternal<double,1014>::ValueBufferType& buffer) ;
4352// to force instanciation may be there is a better way to do it
4353void force_int_instance()
4354{
4355 SYCLBEllPackMatrix<double> matrix ;
4356 SYCLBEllPackMatrix<double>::MatrixInternalType::ValueBufferType buffer{sycl::range<1>(10)} ;
4357 matrix.internal()->setMatrixValues(buffer) ;
4358 std::size_t pack_size =SYCLBEllPackMatrix<double>::MatrixInternalType::ellpack_size ;
4359}
4360/*---------------------------------------------------------------------------*/
4361
4362} // namespace Alien
4363
4364/*---------------------------------------------------------------------------*/
Alien::IMatrixImpl::IMatrixImpl
IMatrixImpl(const MultiMatrixImpl *multi_impl, BackEndId backend="")
Constructor.
Definition IMatrixImpl.cc:38
Alien::MultiMatrixImpl
Multi matrices representation container.
Definition MultiMatrixImpl.h:69
Alien
-- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature --
Definition BackEnd.h:17