Arcane  v4.1.10.0
Documentation développeur
Chargement...
Recherche...
Aucune correspondance
MatrixOperationsImpl.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/* MatrixOperationsImpl.h (C) 2000-2026 */
9/* */
10/* Sparse matrix operations for matrices that provide row_iterator. */
11/*---------------------------------------------------------------------------*/
12#ifndef ARCCORE_ALINA_MATRIXOPERATIONSIMPL_H
13#define ARCCORE_ALINA_MATRIXOPERATIONSIMPL_H
14/*---------------------------------------------------------------------------*/
15/*---------------------------------------------------------------------------*/
16/*
17 * This file is based on the work on AMGCL library (version march 2026)
18 * which can be found at https://github.com/ddemidov/amgcl.
19 *
20 * Copyright (c) 2012-2022 Denis Demidov <dennis.demidov@gmail.com>
21 * SPDX-License-Identifier: MIT
22 */
23/*---------------------------------------------------------------------------*/
24/*---------------------------------------------------------------------------*/
25
26#include "arccore/alina/BackendInterface.h"
27
28/*---------------------------------------------------------------------------*/
29/*---------------------------------------------------------------------------*/
30
31namespace Arcane::Alina::backend
32{
33
34/*---------------------------------------------------------------------------*/
35/*---------------------------------------------------------------------------*/
36
37namespace detail
38{
39
40 template <class Matrix, class Enable = void>
41 struct use_builtin_matrix_ops : std::false_type
42 {};
43
44} // namespace detail
45
46/*---------------------------------------------------------------------------*/
47/*---------------------------------------------------------------------------*/
48
49template <class Alpha, class Matrix, class Vector1, class Beta, class Vector2>
50struct spmv_impl<Alpha, Matrix, Vector1, Beta, Vector2,
51 std::enable_if_t<
52 detail::use_builtin_matrix_ops<Matrix>::value &&
53 math::static_rows<typename value_type<Matrix>::type>::value == math::static_rows<typename value_type<Vector1>::type>::value &&
54 math::static_rows<typename value_type<Matrix>::type>::value == math::static_rows<typename value_type<Vector2>::type>::value>>
55{
56 static void apply(Alpha alpha, const Matrix& A, const Vector1& x, Beta beta, Vector2& y)
57 {
58 typedef typename value_type<Vector2>::type V;
59
60 const ptrdiff_t n = static_cast<ptrdiff_t>(nbRow(A));
61
62 arccoreParallelFor(0, n, ForLoopRunInfo{}, [&](Int32 begin, Int32 size) {
63 if (!math::is_zero(beta)) {
64 for (ptrdiff_t i = begin; i < (begin + size); ++i) {
65 V sum = math::zero<V>();
66 for (typename row_iterator<Matrix>::type a = row_begin(A, i); a; ++a)
67 sum += a.value() * x[a.col()];
68 y[i] = alpha * sum + beta * y[i];
69 }
70 }
71 else {
72 for (ptrdiff_t i = begin; i < (begin + size); ++i) {
73 V sum = math::zero<V>();
74 for (typename row_iterator<Matrix>::type a = row_begin(A, i); a; ++a)
75 sum += a.value() * x[a.col()];
76 y[i] = alpha * sum;
77 }
78 }
79 });
80 }
81};
82
83/*---------------------------------------------------------------------------*/
84/*---------------------------------------------------------------------------*/
85
86template <class Matrix, class Vector1, class Vector2, class Vector3>
87struct residual_impl<Matrix, Vector1, Vector2, Vector3,
88 std::enable_if_t<detail::use_builtin_matrix_ops<Matrix>::value &&
89 math::static_rows<typename value_type<Matrix>::type>::value == math::static_rows<typename value_type<Vector1>::type>::value &&
90 math::static_rows<typename value_type<Matrix>::type>::value == math::static_rows<typename value_type<Vector2>::type>::value &&
91 math::static_rows<typename value_type<Matrix>::type>::value == math::static_rows<typename value_type<Vector3>::type>::value>>
92{
93 static void apply(Vector1 const& rhs,
94 Matrix const& A,
95 Vector2 const& x,
96 Vector3& res)
97 {
98 typedef typename value_type<Vector3>::type V;
99 const ptrdiff_t n = static_cast<ptrdiff_t>(nbRow(A));
100
101 arccoreParallelFor(0, n, ForLoopRunInfo{}, [&](Int32 begin, Int32 size) {
102 for (ptrdiff_t i = begin; i < (begin+size); ++i) {
103 V sum = math::zero<V>();
104 for (typename row_iterator<Matrix>::type a = row_begin(A, i); a; ++a)
105 sum += a.value() * x[a.col()];
106 res[i] = rhs[i] - sum;
107 }
108 });
109 }
110};
111
112/*---------------------------------------------------------------------------*/
113/*---------------------------------------------------------------------------*/
114/*
115 * Allows to do matrix-vector products with mixed scalar/nonscalar types.
116 * Reinterprets pointers to the vectors data into appropriate types.
117 */
118template <class Alpha, class Matrix, class Vector1, class Beta, class Vector2>
119struct spmv_impl<Alpha, Matrix, Vector1, Beta, Vector2,
120 std::enable_if_t<detail::use_builtin_matrix_ops<Matrix>::value &&
121 (math::static_rows<typename value_type<Matrix>::type>::value != math::static_rows<typename value_type<Vector1>::type>::value ||
122 math::static_rows<typename value_type<Matrix>::type>::value != math::static_rows<typename value_type<Vector2>::type>::value)>>
123{
124 static void apply(Alpha alpha, const Matrix& A, const Vector1& x, Beta beta, Vector2& y)
125 {
126 typedef typename value_type<Matrix>::type V;
127
128 auto X = backend::reinterpret_as_rhs<V>(x);
129 auto Y = backend::reinterpret_as_rhs<V>(y);
130
131 spmv(alpha, A, X, beta, Y);
132 }
133};
134
135/*---------------------------------------------------------------------------*/
136/*---------------------------------------------------------------------------*/
137
138template <class Matrix, class Vector1, class Vector2, class Vector3>
139struct residual_impl<Matrix, Vector1, Vector2, Vector3,
140 std::enable_if_t<detail::use_builtin_matrix_ops<Matrix>::value &&
141 (math::static_rows<typename value_type<Matrix>::type>::value != math::static_rows<typename value_type<Vector1>::type>::value ||
142 math::static_rows<typename value_type<Matrix>::type>::value != math::static_rows<typename value_type<Vector2>::type>::value ||
143 math::static_rows<typename value_type<Matrix>::type>::value != math::static_rows<typename value_type<Vector3>::type>::value)>>
144{
145 static void apply(Vector1 const& f,
146 Matrix const& A,
147 Vector2 const& x,
148 Vector3& r)
149 {
150 typedef typename value_type<Matrix>::type V;
151
152 auto X = backend::reinterpret_as_rhs<V>(x);
153 auto F = backend::reinterpret_as_rhs<V>(f);
154 auto R = backend::reinterpret_as_rhs<V>(r);
155
156 residual(F, A, X, R);
157 }
158};
159
160/*---------------------------------------------------------------------------*/
161/*---------------------------------------------------------------------------*/
162
163} // namespace Arcane::Alina::backend
164
165/*---------------------------------------------------------------------------*/
166/*---------------------------------------------------------------------------*/
167
168#endif
Arcane::ForLoopRunInfo
Informations d'exécution d'une boucle.
Definition ForLoopRunInfo.h:37
Arcane::Matrix
Matrix class, to be used by user.
Definition matrix/Matrix.h:36
Arcane::Vector2
Classe gérant un vecteur de dimension 2 de type T.
Definition Vector2.h:36
Arcane::Vector3
Classe gérant un vecteur de dimension 3 de type T.
Definition Vector3.h:36
Arcane::arccoreParallelFor
void arccoreParallelFor(const ComplexForLoopRanges< RankValue > &loop_ranges, const ForLoopRunInfo &run_info, const LambdaType &lambda_function, const ReducerArgs &... reducer_args)
Applique en concurrence la fonction lambda lambda_function sur l'intervalle d'itération donné par loo...
Definition ParallelFor.h:85
Arcane::Int32
std::int32_t Int32
Type entier signé sur 32 bits.
Definition ArccoreGlobal.h:225
Arcane::Alina::backend::detail::use_builtin_matrix_ops
Definition MatrixOperationsImpl.h:42
Arcane::Alina::backend::residual_impl
Implementation for residual error compuatation.
Definition BackendInterface.h:263
Arcane::Alina::backend::spmv_impl
Implementation for matrix-vector product.
Definition BackendInterface.h:250