Alien  1.3.0
Developer documentation
Loading...
Searching...
No Matches
SimpleCSRMatrixMultT.h
1// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
2//-----------------------------------------------------------------------------
3// Copyright 2000-2026 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
4// See the top-level COPYRIGHT file for details.
5// SPDX-License-Identifier: Apache-2.0
6//-----------------------------------------------------------------------------
7
8#pragma once
9
10#include <arccore/collections/Array2.h>
11
12#include <alien/handlers/scalar/CSRModifierViewT.h>
13
14/*---------------------------------------------------------------------------*/
15
16namespace Alien::SimpleCSRInternal
17{
18
19/*---------------------------------------------------------------------------*/
20
21template <typename ValueT>
22SimpleCSRMatrixMultT<ValueT>::SimpleCSRMatrixMultT(const MatrixType& matrix)
23: m_matrix_impl(matrix)
24{}
25
26/*---------------------------------------------------------------------------*/
27
28template <typename ValueT>
29void SimpleCSRMatrixMultT<ValueT>::synchronize(VectorType& x) const
30{
31 if (m_matrix_impl.block()) {
32 if (m_matrix_impl.m_is_parallel)
33 _synchronizeBlock(x);
34 }
35 else if (m_matrix_impl.vblock()) {
36 if (m_matrix_impl.m_is_parallel)
37 _synchronizeVariableBlock(x);
38 }
39 else {
40 if (m_matrix_impl.m_is_parallel)
41 _synchronize(x);
42 }
43}
44
45
46
47
48template <typename ValueT>
49void SimpleCSRMatrixMultT<ValueT>::mult(const VectorType& x, VectorType& y) const
50{
51 if (m_matrix_impl.block()) {
52 if (m_matrix_impl.m_is_parallel)
53 _parallelMultBlock(x, y);
54 else
55 _seqMultBlock(x, y);
56 }
57 else if (m_matrix_impl.vblock()) {
58 if (m_matrix_impl.m_is_parallel)
59 _parallelMultVariableBlock(x, y);
60 else
61 _seqMultVariableBlock(x, y);
62 }
63 else {
64 if (m_matrix_impl.m_is_parallel)
65 _parallelMult(x, y);
66 else
67 _seqMult(x, y);
68 }
69}
70
71template <typename ValueT>
72void SimpleCSRMatrixMultT<ValueT>::addLMult(Real alpha, const VectorType& x, VectorType& y) const
73{
74 _seqAddLMult(alpha, x, y);
75}
76
77template <typename ValueT>
78void SimpleCSRMatrixMultT<ValueT>::addUMult(Real alpha, const VectorType& x, VectorType& y) const
79{
80 _seqAddUMult(alpha, x, y);
81}
82
83template <typename ValueT>
84void SimpleCSRMatrixMultT<ValueT>::mult(const UniqueArray<Real>& x, UniqueArray<Real>& y) const
85{
86 if (m_matrix_impl.m_is_parallel)
87 _parallelMult(x, y);
88 else
89 _seqMult(x, y);
90}
91
92
93template <typename ValueT>
94void SimpleCSRMatrixMultT<ValueT>::lowerUpperSynchronize(const VectorType& x_impl, VectorType& y_impl) const
95{
96 Real* x_ptr = const_cast<Real*>(x_impl.getDataPtr());
97 Real* y_ptr = y_impl.getDataPtr();
98 SendRecvOp<Real> op(x_ptr,
99 m_matrix_impl.m_matrix_dist_info.m_send_info,
100 m_matrix_impl.m_send_policy,
101 y_ptr,
102 m_matrix_impl.m_matrix_dist_info.m_recv_info,
103 m_matrix_impl.m_recv_policy,
104 m_matrix_impl.m_parallel_mng,
105 m_matrix_impl.m_trace);
106 op.lowerRecv();
107 op.upperSend();
108}
109
110
111
112template <typename ValueT>
113void SimpleCSRMatrixMultT<ValueT>::upperLowerSynchronize(const VectorType& x_impl, VectorType& y_impl) const
114{
115 Real* x_ptr = const_cast<Real*>(x_impl.getDataPtr());
116 Real* y_ptr = y_impl.getDataPtr();
117 SendRecvOp<Real> op(x_ptr,
118 m_matrix_impl.m_matrix_dist_info.m_send_info,
119 m_matrix_impl.m_send_policy,
120 y_ptr,
121 m_matrix_impl.m_matrix_dist_info.m_recv_info,
122 m_matrix_impl.m_recv_policy,
123 m_matrix_impl.m_parallel_mng,
124 m_matrix_impl.m_trace);
125 op.upperRecv();
126 op.lowerSend();
127}
128/*---------------------------------------------------------------------------*/
129
130template <typename ValueT>
131void SimpleCSRMatrixMultT<ValueT>::_synchronize(VectorType& x_impl) const
132{
133 Integer alloc_size = m_matrix_impl.m_local_size + m_matrix_impl.m_ghost_size;
134 x_impl.resize(alloc_size);
135 Real* x_ptr = (Real*)x_impl.getDataPtr();
136 ConstArrayView<Real> matrix = m_matrix_impl.m_matrix.getValues();
137 ConstArrayView<Integer> cols = m_matrix_impl.getDistStructInfo().m_cols;
138 ConstArrayView<Integer> row_offset =
139 m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
140 SendRecvOp<Real> op(x_ptr, m_matrix_impl.m_matrix_dist_info.m_send_info,
141 m_matrix_impl.m_send_policy, x_ptr, m_matrix_impl.m_matrix_dist_info.m_recv_info,
142 m_matrix_impl.m_recv_policy, m_matrix_impl.m_parallel_mng, m_matrix_impl.m_trace);
143 op.start();
144 op.end();
145}
146
147template <typename ValueT>
148void SimpleCSRMatrixMultT<ValueT>::_parallelMult(
149const VectorType& x_impl, VectorType& y_impl) const
150{
151 Integer alloc_size = m_matrix_impl.m_local_size + m_matrix_impl.m_ghost_size;
152 x_impl.resize(alloc_size);
153 Real* y_ptr = y_impl.getDataPtr();
154 Real* x_ptr = (Real*)x_impl.getDataPtr();
155 ConstArrayView<Real> matrix = m_matrix_impl.m_matrix.getValues();
156 // ConstArrayView<Integer> cols2 =
157 // m_matrix_impl.m_matrix.getCSRProfile().getCols();
158 ConstArrayView<Integer> cols = m_matrix_impl.getDistStructInfo().m_cols;
159 ConstArrayView<Integer> row_offset =
160 m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
161 SendRecvOp<Real> op(x_ptr, m_matrix_impl.m_matrix_dist_info.m_send_info,
162 m_matrix_impl.m_send_policy, x_ptr, m_matrix_impl.m_matrix_dist_info.m_recv_info,
163 m_matrix_impl.m_recv_policy, m_matrix_impl.m_parallel_mng, m_matrix_impl.m_trace);
164 op.start();
165 ConstArrayView<Integer> local_row_size =
166 m_matrix_impl.m_matrix_dist_info.m_local_row_size;
167 for (Integer irow = 0; irow < m_matrix_impl.m_local_size; ++irow) {
168 Integer off = row_offset[irow];
169 Integer off2 = off + local_row_size[irow];
170 Real tmpy = 0.;
171 for (Integer j = off; j < off2; ++j) {
172 tmpy += matrix[j] * x_ptr[cols[j]];
173 }
174 y_ptr[irow] = tmpy;
175 }
176 op.end();
177
178 Integer interface_nrow = m_matrix_impl.m_matrix_dist_info.m_interface_nrow;
179 ConstArrayView<Integer> row_ids = m_matrix_impl.m_matrix_dist_info.m_interface_rows;
180 for (Integer i = 0; i < interface_nrow; ++i) {
181 Integer irow = row_ids[i];
182 Integer off = row_offset[irow] + local_row_size[irow];
183 Integer off2 = row_offset[irow + 1];
184 Real tmpy = 0.;
185 for (Integer j = off; j < off2; ++j) {
186 tmpy += matrix[j] * x_ptr[cols[j]];
187 }
188 y_ptr[irow] += tmpy;
189 }
190}
191
192template <typename ValueT>
193void SimpleCSRMatrixMultT<ValueT>::_parallelMult(
194const UniqueArray<Real>& x_impl, UniqueArray<Real>& y_impl) const
195{
196 Real* y_ptr = dataPtr(y_impl);
197 Real* x_ptr = (Real*)dataPtr(x_impl);
198 ConstArrayView<Real> matrix = m_matrix_impl.m_matrix.getValues();
199 ConstArrayView<Integer> cols = m_matrix_impl.getDistStructInfo().m_cols;
200 ConstArrayView<Integer> row_offset =
201 m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
202 SendRecvOp<Real> op(x_ptr, m_matrix_impl.m_matrix_dist_info.m_send_info,
203 m_matrix_impl.m_send_policy, x_ptr, m_matrix_impl.m_matrix_dist_info.m_recv_info,
204 m_matrix_impl.m_recv_policy, m_matrix_impl.m_parallel_mng, m_matrix_impl.m_trace);
205 op.start();
206 ConstArrayView<Integer> local_row_size =
207 m_matrix_impl.m_matrix_dist_info.m_local_row_size;
208 for (Integer irow = 0; irow < m_matrix_impl.m_local_size; ++irow) {
209 Integer off = row_offset[irow];
210 Integer off2 = off + local_row_size[irow];
211 Real tmpy = 0.;
212 for (Integer j = off; j < off2; ++j) {
213 tmpy += matrix[j] * x_ptr[cols[j]];
214 }
215 y_ptr[irow] = tmpy;
216 }
217 op.end();
218
219 Integer interface_nrow = m_matrix_impl.m_matrix_dist_info.m_interface_nrow;
220 ConstArrayView<Integer> row_ids = m_matrix_impl.m_matrix_dist_info.m_interface_rows;
221 for (Integer i = 0; i < interface_nrow; ++i) {
222 Integer irow = row_ids[i];
223 Integer off = row_offset[irow] + local_row_size[irow];
224 Integer off2 = row_offset[irow + 1];
225 Real tmpy = 0.;
226 for (Integer j = off; j < off2; ++j) {
227 tmpy += matrix[j] * x_ptr[cols[j]];
228 }
229 y_ptr[irow] += tmpy;
230 }
231}
232/*---------------------------------------------------------------------------*/
233
234template <typename ValueT>
235void SimpleCSRMatrixMultT<ValueT>::_seqMult(const VectorType& x_impl, VectorType& y_impl) const
236{
237#ifdef ALIEN_USE_PERF_TIMER
238 typename MatrixType::SentryType sentry(m_matrix_impl.timer(), "CSR-SPMV");
239#endif
240 Real* y_ptr = y_impl.getDataPtr();
241 Real* x_ptr = (Real*)x_impl.getDataPtr();
242 ConstArrayView<Real> matrix = m_matrix_impl.m_matrix.getValues();
243 ConstArrayView<Integer> cols = m_matrix_impl.m_matrix.getCSRProfile().getCols();
244 ConstArrayView<Integer> row_offset =
245 m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
246 for (Integer irow = 0; irow < m_matrix_impl.m_local_size; ++irow) {
247 Real tmpy = 0.;
248 for (Integer j = row_offset[irow]; j < row_offset[irow + 1]; ++j) {
249 tmpy += matrix[j] * x_ptr[cols[j]];
250 }
251 y_ptr[irow] = tmpy;
252 }
253}
254
255template <typename ValueT>
256void SimpleCSRMatrixMultT<ValueT>::_seqAddLMult(Real alpha, const VectorType& x_impl, VectorType& y_impl) const
257{
258#ifdef ALIEN_USE_PERF_TIMER
259 typename MatrixType::SentryType sentry(m_matrix_impl.timer(), "CSR-AddLMult");
260#endif
261 Real* y_ptr = y_impl.getDataPtr();
262 Real* x_ptr = (Real*)x_impl.getDataPtr();
263 ConstArrayView<Real> matrix = m_matrix_impl.m_matrix.getValues();
264 ConstArrayView<Integer> cols = m_matrix_impl.m_matrix.getCSRProfile().getCols();
265 ConstArrayView<Integer> row_offset = m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
266 auto diag_offset = m_matrix_impl.m_matrix.getCSRProfile().getUpperDiagOffset();
267 for (Integer irow = 0; irow < m_matrix_impl.m_local_size; ++irow) {
268 Real tmpy = y_ptr[irow];
269 for (Integer j = row_offset[irow]; j < diag_offset[irow]; ++j) {
270 tmpy += alpha * matrix[j] * x_ptr[cols[j]];
271 }
272 y_ptr[irow] = tmpy;
273 }
274}
275
276template <typename ValueT>
277void SimpleCSRMatrixMultT<ValueT>::_seqAddUMult(Real alpha, const VectorType& x_impl, VectorType& y_impl) const
278{
279#ifdef ALIEN_USE_PERF_TIMER
280 typename MatrixType::SentryType sentry(m_matrix_impl.timer(), "CSR-AddUMult");
281#endif
282 Real* y_ptr = y_impl.getDataPtr();
283 Real* x_ptr = (Real*)x_impl.getDataPtr();
284 ConstArrayView<Real> matrix = m_matrix_impl.m_matrix.getValues();
285 ConstArrayView<Integer> cols = m_matrix_impl.m_matrix.getCSRProfile().getCols();
286 ConstArrayView<Integer> row_offset = m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
287 auto diag_offset = m_matrix_impl.m_matrix.getCSRProfile().getUpperDiagOffset();
288 for (Integer irow = 0; irow < m_matrix_impl.m_local_size; ++irow) {
289 Real tmpy = y_ptr[irow];
290 for (Integer j = diag_offset[irow] + 1; j < row_offset[irow + 1]; ++j) {
291 tmpy += alpha * matrix[j] * x_ptr[cols[j]];
292 }
293 y_ptr[irow] = tmpy;
294 }
295}
296
297template <typename ValueT>
298void SimpleCSRMatrixMultT<ValueT>::_seqMult(
299const UniqueArray<Real>& x_impl, UniqueArray<Real>& y_impl) const
300{
301 Real* y_ptr = dataPtr(y_impl);
302 Real* x_ptr = (Real*)dataPtr(x_impl);
303 ConstArrayView<Real> matrix = m_matrix_impl.m_matrix.getValues();
304 ConstArrayView<Integer> cols = m_matrix_impl.m_matrix.getCSRProfile().getCols();
305 ConstArrayView<Integer> row_offset =
306 m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
307 for (Integer irow = 0; irow < m_matrix_impl.m_local_size; ++irow) {
308 Real tmpy = 0.;
309 for (Integer j = row_offset[irow]; j < row_offset[irow + 1]; ++j) {
310 tmpy += matrix[j] * x_ptr[cols[j]];
311 }
312 y_ptr[irow] = tmpy;
313 }
314}
315
316/*---------------------------------------------------------------------------*/
317
318template <typename ValueT>
319void SimpleCSRMatrixMultT<ValueT>::_synchronizeBlock(VectorType& x) const
320{
321 Integer alloc_size = m_matrix_impl.m_local_size + m_matrix_impl.m_ghost_size;
322 const Integer block_size = m_matrix_impl.block()->size();
323 x.resize(alloc_size * block_size);
324 ConstArrayView<Real> x_ptr = x.fullValues();
325 Real const* matrix = m_matrix_impl.m_matrix.getDataPtr();
326 ConstArrayView<Integer> cols = m_matrix_impl.getDistStructInfo().m_cols;
327 ConstArrayView<Integer> row_offset =
328 m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
329 SimpleCSRInternal::SendRecvOp<Real> op(x.getDataPtr(),
330 m_matrix_impl.m_matrix_dist_info.m_send_info, m_matrix_impl.m_send_policy,
331 (Real*)x.getDataPtr(), m_matrix_impl.m_matrix_dist_info.m_recv_info,
332 m_matrix_impl.m_recv_policy, m_matrix_impl.m_parallel_mng, m_matrix_impl.m_trace,
333 block_size);
334 op.start();
335 op.end();
336}
337
338
339
340
341
342template <typename ValueT>
343void SimpleCSRMatrixMultT<ValueT>::_parallelMultBlock(const VectorType& x, VectorType& y) const
344{
345 Integer alloc_size = m_matrix_impl.m_local_size + m_matrix_impl.m_ghost_size;
346 const Integer block_size = m_matrix_impl.block()->size();
347 x.resize(alloc_size * block_size);
348 ArrayView<Real> _y = y.fullValues();
349 ConstArrayView<Real> x_ptr = x.fullValues();
350 Real const* matrix = m_matrix_impl.m_matrix.getDataPtr();
351 ConstArrayView<Integer> cols = m_matrix_impl.getDistStructInfo().m_cols;
352 // ConstArrayView<Integer> cols2 =
353 // m_matrix_impl.m_matrix.getCSRProfile().getCols();
354 ConstArrayView<Integer> row_offset =
355 m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
356 SimpleCSRInternal::SendRecvOp<Real> op(x.getDataPtr(),
357 m_matrix_impl.m_matrix_dist_info.m_send_info, m_matrix_impl.m_send_policy,
358 (Real*)x.getDataPtr(), m_matrix_impl.m_matrix_dist_info.m_recv_info,
359 m_matrix_impl.m_recv_policy, m_matrix_impl.m_parallel_mng, m_matrix_impl.m_trace,
360 block_size);
361 op.start();
362 ConstArrayView<Integer> local_row_size =
363 m_matrix_impl.m_matrix_dist_info.m_local_row_size;
364 for (Integer irow = 0; irow < m_matrix_impl.m_local_size; ++irow) {
365 ArrayView<Real> y_ptr = _y.subView(irow * block_size, block_size);
366 Integer off = row_offset[irow];
367 Integer off2 = off + local_row_size[irow];
368 Real const* m = matrix + off * block_size * block_size;
369 for (Integer ieq = 0; ieq < block_size; ++ieq)
370 y_ptr[ieq] = 0.;
371 for (Integer jcol = off; jcol < off2; ++jcol) {
372 ConstArrayView<Real> ptr = x_ptr.subView(cols[jcol] * block_size, block_size);
373 for (Integer ieq = 0; ieq < block_size; ++ieq)
374 for (Integer iu = 0; iu < block_size; ++iu)
375 y_ptr[ieq] += m[iu + block_size * ieq] * ptr[iu];
376 m += block_size * block_size;
377 }
378 // y_ptr += block_size;
379 }
380 op.end();
381
382 Integer interface_nrow = m_matrix_impl.m_matrix_dist_info.m_interface_nrow;
383 ConstArrayView<Integer> row_ids = m_matrix_impl.m_matrix_dist_info.m_interface_rows;
384 ArrayView<Real> y_ptr = _y;
385 for (Integer i = 0; i < interface_nrow; ++i) {
386 Integer irow = row_ids[i];
387 ArrayView<Real> yptr = y_ptr.subView(irow * block_size, block_size);
388 Integer off = row_offset[irow] + local_row_size[irow];
389 Integer off2 = row_offset[irow + 1];
390 Real const* m = matrix + off * block_size * block_size;
391 for (Integer jcol = off; jcol < off2; ++jcol) {
392 ConstArrayView<Real> ptr = x_ptr.subView(cols[jcol] * block_size, block_size);
393 for (Integer ieq = 0; ieq < block_size; ++ieq)
394 for (Integer iu = 0; iu < block_size; ++iu)
395 yptr[ieq] += m[iu + block_size * ieq] * ptr[iu];
396 m += block_size * block_size;
397 }
398 }
399}
400
401/*---------------------------------------------------------------------------*/
402
403template <typename ValueT>
404void SimpleCSRMatrixMultT<ValueT>::_seqMultBlock(const VectorType& x, VectorType& y) const
405{
406 Real* y_ptr = y.getDataPtr();
407 Real const* x_ptr = x.getDataPtr();
408 Real const* matrix = m_matrix_impl.m_matrix.getDataPtr();
409 const Integer block_size = m_matrix_impl.block()->size();
410 ConstArrayView<Integer> cols = m_matrix_impl.m_matrix.getCSRProfile().getCols();
411 ConstArrayView<Integer> row_offset = m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
412 Real* yptr = y_ptr ;
413 for (Integer irow = 0; irow < m_matrix_impl.m_local_size; ++irow)
414 {
415 Integer off = row_offset[irow];
416 Integer off2 = row_offset[irow + 1];
417 Real const* m = matrix + off * block_size * block_size;
418 for (Integer ieq = 0; ieq < block_size; ++ieq)
419 yptr[ieq] = 0.;
420 for (Integer jcol = off; jcol < off2; ++jcol)
421 {
422 Real const* ptr = x_ptr + cols[jcol] * block_size;
423 for (Integer ieq = 0; ieq < block_size; ++ieq)
424 for (Integer iu = 0; iu < block_size; ++iu)
425 yptr[ieq] += m[iu + block_size * ieq] * ptr[iu];
426 m += block_size * block_size;
427 }
428 yptr += block_size;
429 }
430
431 /*
432 for(int il=0;il<m_matrix_impl.m_local_size;++il)
433 {
434 std::cout<<"X["<<il<<"]:";
435 for(int i=0;i<block_size;++i)
436 std::cout<<x_ptr[il*block_size+i]<<",";
437 std::cout<<std::endl;
438 }
439 for(int il=0;il<m_matrix_impl.m_local_size;++il)
440 {
441 std::cout<<"Y["<<il<<"]:";
442 for(int i=0;i<block_size;++i)
443 std::cout<<y_ptr[il*block_size+i]<<",";
444 std::cout<<std::endl;
445 }
446 for(std::size_t irow=0;irow<m_matrix_impl.m_local_size;++irow)
447 {
448 for(int ieq=0;ieq<block_size;++ieq)
449 {
450 std::cout<<"LINE["<<irow<<","<<ieq<<"] : ";
451 ValueType value = 0. ;
452 Integer off = row_offset[irow];
453 Integer off2 = row_offset[irow + 1];
454 Real const* m = matrix + off * block_size * block_size;
455 for (Integer jcol = off; jcol < off2; ++jcol)
456 {
457 Real const* ptr = x_ptr + cols[jcol] * block_size;
458 for (Integer iu = 0; iu < block_size; ++iu)
459 value += m[iu + block_size * ieq] * ptr[iu];
460 m += block_size * block_size;
461 }
462 std::cout<<"\nY_CPU["<<irow<<","<<ieq<<"]:"<<value<<std::endl;
463 }
464 }*/
465}
466
467/*---------------------------------------------------------------------------*/
468template <typename ValueT>
469void SimpleCSRMatrixMultT<ValueT>::_synchronizeVariableBlock(VectorType& x_impl) const
470{
471 // alien_info([&] { cout()<<"_parallelMultVariableBlock";}) ;
472
473 ConstArrayView<Integer> block_sizes = m_matrix_impl.getDistStructInfo().m_block_sizes;
474 ConstArrayView<Integer> block_offsets = m_matrix_impl.getDistStructInfo().m_block_offsets;
475 {
476 const Integer last = block_offsets.size() - 1;
477 x_impl.resize(block_offsets[last]);
478 }
479
480 const ValueT* x_ptr = x_impl.getDataPtr();
481 const ValueT* matrix_ptr = m_matrix_impl.m_matrix.getDataPtr();
482
483 ConstArrayView<Integer> row_offset =
484 m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
485 ConstArrayView<Integer> cols = m_matrix_impl.getDistStructInfo().m_cols;
486 ConstArrayView<Integer> block_cols =
487 m_matrix_impl.m_matrix.getCSRProfile().getBlockCols();
488
489 SendRecvOp<Real> op(x_ptr, m_matrix_impl.m_matrix_dist_info.m_send_info,
490 m_matrix_impl.m_send_policy, (ValueT*)x_ptr,
491 m_matrix_impl.m_matrix_dist_info.m_recv_info, m_matrix_impl.m_recv_policy,
492 m_matrix_impl.m_parallel_mng, m_matrix_impl.m_trace, block_sizes, block_offsets);
493
494 op.start();
495 op.end();
496}
497
498
499
500
501
502template <typename ValueT>
503void SimpleCSRMatrixMultT<ValueT>::_parallelMultVariableBlock(
504const VectorType& x_impl, VectorType& y_impl) const
505{
506 // alien_info([&] { cout()<<"_parallelMultVariableBlock";}) ;
507
508 ArrayView<ValueT> y = y_impl.fullValues();
509
510 ConstArrayView<Integer> block_sizes = m_matrix_impl.getDistStructInfo().m_block_sizes;
511 ConstArrayView<Integer> block_offsets = m_matrix_impl.getDistStructInfo().m_block_offsets;
512 {
513 const Integer last = block_offsets.size() - 1;
514 x_impl.resize(block_offsets[last]);
515 }
516
517 const ValueT* x_ptr = x_impl.getDataPtr();
518 const ValueT* matrix_ptr = m_matrix_impl.m_matrix.getDataPtr();
519
520 ConstArrayView<Integer> row_offset =
521 m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
522 ConstArrayView<Integer> cols = m_matrix_impl.getDistStructInfo().m_cols;
523 ConstArrayView<Integer> block_cols =
524 m_matrix_impl.m_matrix.getCSRProfile().getBlockCols();
525
526 SendRecvOp<Real> op(x_ptr, m_matrix_impl.m_matrix_dist_info.m_send_info,
527 m_matrix_impl.m_send_policy, (ValueT*)x_ptr,
528 m_matrix_impl.m_matrix_dist_info.m_recv_info, m_matrix_impl.m_recv_policy,
529 m_matrix_impl.m_parallel_mng, m_matrix_impl.m_trace, block_sizes, block_offsets);
530
531 op.start();
532
533 UniqueArray<ValueT> tmpy;
534 tmpy.reserve(m_matrix_impl.vblock()->maxBlockSize());
535
536 ConstArrayView<Integer> local_row_size =
537 m_matrix_impl.m_matrix_dist_info.m_local_row_size;
538 for (Integer irow = 0; irow < m_matrix_impl.m_local_size; ++irow) {
539 Integer off = row_offset[irow];
540 Integer off2 = off + local_row_size[irow];
541 const Integer block_size_row = block_sizes[irow];
542 tmpy.resize(block_size_row);
543 tmpy.fill(ValueT());
544 for (Integer j = off; j < off2; ++j) {
545 const Integer col = cols[j];
546 const Integer block_size_col = block_sizes[col];
547 ConstArrayView<ValueT> x(block_size_col, x_ptr + block_offsets[col]);
548 ConstArray2View<ValueT> matrix(
549 matrix_ptr + block_cols[j], block_size_row, block_size_col);
550 for (Integer krow = 0; krow < block_size_row; ++krow) {
551 for (Integer kcol = 0; kcol < block_size_col; ++kcol) {
552 tmpy[krow] += matrix[krow][kcol] * x[kcol];
553 }
554 }
555 }
556 y.subView(block_offsets[irow], block_size_row).copy(tmpy);
557 }
558
559 op.end();
560
561 Integer interface_nrow = m_matrix_impl.m_matrix_dist_info.m_interface_nrow;
562 ConstArrayView<Integer> row_ids = m_matrix_impl.m_matrix_dist_info.m_interface_rows;
563 for (Integer i = 0; i < interface_nrow; ++i) {
564 Integer irow = row_ids[i];
565 Integer off = row_offset[irow] + local_row_size[irow];
566 Integer off2 = row_offset[irow + 1];
567 const Integer block_size_row = block_sizes[irow];
568 tmpy.resize(block_size_row);
569 tmpy.fill(ValueT());
570 for (Integer j = off; j < off2; ++j) {
571 const Integer col = cols[j];
572 const Integer block_size_col = block_sizes[col];
573 ConstArrayView<ValueT> x(block_size_col, x_ptr + block_offsets[col]);
574 ConstArray2View<ValueT> matrix(
575 matrix_ptr + block_cols[j], block_size_row, block_size_col);
576 for (Integer krow = 0; krow < block_size_row; ++krow) {
577 for (Integer kcol = 0; kcol < block_size_col; ++kcol) {
578 tmpy[krow] += matrix[krow][kcol] * x[kcol];
579 }
580 }
581 }
582 ArrayView<ValueT> y_view = y.subView(block_offsets[irow], block_size_row);
583 for (Integer k = 0; k < block_size_row; ++k)
584 y_view[k] += tmpy[k];
585 }
586}
587
588/*---------------------------------------------------------------------------*/
589
590template <typename ValueT>
591void SimpleCSRMatrixMultT<ValueT>::_seqMultVariableBlock(
592const VectorType& x_impl, VectorType& y_impl) const
593{
594 // alien_info([&] { cout()<<"_seqMultVariableBlock";}) ;
595
596 ArrayView<ValueT> y = y_impl.fullValues();
597
598 const ValueT* x_ptr = x_impl.getDataPtr();
599 const ValueT* matrix_ptr = m_matrix_impl.m_matrix.getDataPtr();
600
601 ConstArrayView<Integer> row_offset =
602 m_matrix_impl.m_matrix.getCSRProfile().getRowOffset();
603 ConstArrayView<Integer> cols = m_matrix_impl.m_matrix.getCSRProfile().getCols();
604 ConstArrayView<Integer> block_cols =
605 m_matrix_impl.m_matrix.getCSRProfile().getBlockCols();
606 const VBlock* block = m_matrix_impl.vblock();
607 const VBlockImpl& block_infos = x_impl.vblockImpl();
608
609 UniqueArray<ValueT> tmpy;
610 tmpy.reserve(block->maxBlockSize());
611 for (Integer irow = 0; irow < m_matrix_impl.m_local_size;
612 ++irow) // Attention, c'est local !!!!
613 {
614 const Integer block_size_row = block->size(irow);
615 tmpy.resize(block_size_row);
616 tmpy.fill(ValueT());
617 for (Integer j = row_offset[irow]; j < row_offset[irow + 1]; ++j) {
618 const Integer col = cols[j];
619 const Integer block_size_col = block->size(col);
620 ConstArrayView<ValueT> x(block_size_col, x_ptr + block_infos.offset(col));
621 ConstArray2View<ValueT> matrix(
622 matrix_ptr + block_cols[j], block_size_row, block_size_col);
623 for (Integer krow = 0; krow < block_size_row; ++krow) {
624 for (Integer kcol = 0; kcol < block_size_col; ++kcol) {
625 tmpy[krow] += matrix[krow][kcol] * x[kcol];
626 }
627 }
628 }
629 y.subView(block_infos.offset(irow), block_size_row).copy(tmpy);
630 }
631}
632
633
634template <typename ValueT>
635void SimpleCSRMatrixMultT<ValueT>::multDiag(VectorType& y) const
636{
637 Real* y_ptr = y.getDataPtr();
638 CSRConstViewT<MatrixType> view(m_matrix_impl);
639 // clang-format off
640 auto nrows = view.nrows() ;
641 auto kcol = view.kcol() ;
642 auto dcol = view.dcol() ;
643 auto cols = view.cols() ;
644 auto values = view.data() ;
645 // clang-format on
646 for (Integer irow = 0; irow < nrows; ++irow)
647 y_ptr[irow] = y_ptr[irow] * values[dcol[irow]];
648}
649
650template <typename ValueT>
651void SimpleCSRMatrixMultT<ValueT>::multInvDiag(VectorType& y) const
652{
653 Real* y_ptr = y.getDataPtr();
654 CSRConstViewT<MatrixType> view(m_matrix_impl);
655 // clang-format off
656 auto nrows = view.nrows() ;
657 auto kcol = view.kcol() ;
658 auto dcol = view.dcol() ;
659 auto cols = view.cols() ;
660 auto values = view.data() ;
661 // clang-format on
662 for (Integer irow = 0; irow < nrows; ++irow)
663 y_ptr[irow] = y_ptr[irow] / values[dcol[irow]];
664}
665
666template <typename ValueT>
667void SimpleCSRMatrixMultT<ValueT>::computeDiag(VectorType& y) const
668{
669 Real* y_ptr = y.getDataPtr();
670
671 CSRConstViewT<MatrixType> view(m_matrix_impl);
672 // clang-format off
673 auto nrows = view.nrows() ;
674 auto kcol = view.kcol() ;
675 auto dcol = view.dcol() ;
676 auto cols = view.cols() ;
677 auto values = view.data() ;
678 // clang-format on
679 if(m_matrix_impl.blockSize()==1)
680 {
681 for (Integer irow = 0; irow < nrows; ++irow)
682 y_ptr[irow] = values[dcol[irow]];
683 }
684 else
685 {
686 Integer block_size = m_matrix_impl.blockSize();
687 Integer block2_size = block_size*block_size ;
688 for (Integer irow = 0; irow < nrows; ++irow)
689 {
690 for(Integer ieq=0;ieq<block_size;++ieq)
691 y_ptr[irow*block_size+ieq] = values[dcol[irow]*block2_size+ieq*block_size+ieq];
692 }
693 }
694}
695
696template <typename ValueT>
697void SimpleCSRMatrixMultT<ValueT>::computeInvDiag(VectorType& y) const
698{
699 Real* y_ptr = y.getDataPtr();
700
701 CSRConstViewT<MatrixType> view(m_matrix_impl);
702 // clang-format off
703 auto nrows = view.nrows() ;
704 auto kcol = view.kcol() ;
705 auto dcol = view.dcol() ;
706 auto cols = view.cols() ;
707 auto values = view.data() ;
708 // clang-format on
709 if(m_matrix_impl.blockSize()==1)
710 {
711 for (Integer irow = 0; irow < nrows; ++irow)
712 y_ptr[irow] = 1. / values[dcol[irow]];
713 }
714 else
715 {
716 Integer block_size = m_matrix_impl.blockSize();
717 Integer block2_size = block_size*block_size ;
718 for (Integer irow = 0; irow < nrows; ++irow)
719 {
720 for(Integer ieq=0;ieq<block_size;++ieq)
721 {
722 y_ptr[irow*block_size+ieq] = 1. / values[dcol[irow]*block2_size+ieq*block_size+ieq];
723 }
724 }
725 /*
726 auto kcol = view.kcol() ;
727 for(int il=0;il<nrows;++il)
728 {
729 std::cout<<"MAT["<<il<<","<<il<<"]:";
730 for(int k=kcol[il];k<kcol[il+1];++k)
731 {
732 std::cout<<"("<<cols[k]<<",";
733 for(int ieq=0;ieq<block_size;++ieq)
734 for(int ju=0;ju<block_size;++ju)
735 std::cout<<values[k*block2_size+ieq*block_size +ju]<<",";
736 std::cout<<")"<<std::endl;
737 }
738 }
739 for(int il=0;il<nrows;++il)
740 {
741 std::cout<<"DIAG["<<il<<"]:";
742 for(int i=0;i<block_size;++i)
743 std::cout<<y_ptr[il*block_size+i]<<",";
744 std::cout<<std::endl;
745 std::cout<<"MAT["<<il<<","<<il<<"]:";
746 auto mat_offset = kcol[il]*block2_size ;
747 for(int ieq=0;ieq<block_size;++ieq)
748 for(int ju=0;ju<block_size;++ju)
749 std::cout<<values[mat_offset+ieq*block_size +ju]<<",";
750 std::cout<<std::endl;
751 }*/
752 }
753}
754
755/*---------------------------------------------------------------------------*/
756
757} // namespace Alien::SimpleCSRInternal
758
759/*---------------------------------------------------------------------------*/
Alien::SimpleCSRInternal::SimpleCSRMatrixMultT::SimpleCSRMatrixMultT
SimpleCSRMatrixMultT(const MatrixType &matrix)
Constructeur de la classe.
Definition SimpleCSRMatrixMultT.h:22
Alien::SimpleCSRInternal::SimpleCSRMatrixMultT::mult
void mult(const VectorType &x, VectorType &y) const
Matrix vector product.
Definition SimpleCSRMatrixMultT.h:49