Arcane  4.1.12.0
Developer documentation
Loading...
Searching...
No Matches
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/*
116 * Allows matrix-vector products with mixed scalar/nonscalar types.
117 * Reinterprets pointers to the vectors data into appropriate types.
118 */
119template <class Alpha, class Matrix, class Vector1, class Beta, class Vector2>
120struct spmv_impl<Alpha, Matrix, Vector1, Beta, Vector2,
121 std::enable_if_t<detail::use_builtin_matrix_ops<Matrix>::value &&
122 (math::static_rows<typename value_type<Matrix>::type>::value != math::static_rows<typename value_type<Vector1>::type>::value ||
123 math::static_rows<typename value_type<Matrix>::type>::value != math::static_rows<typename value_type<Vector2>::type>::value)>>
124{
125 static void apply(Alpha alpha, const Matrix& A, const Vector1& x, Beta beta, Vector2& y)
126 {
127 typedef typename value_type<Matrix>::type V;
128
129 auto X = backend::reinterpret_as_rhs<V>(x);
130 auto Y = backend::reinterpret_as_rhs<V>(y);
131
132 spmv(alpha, A, X, beta, Y);
133 }
134};
135
136/*---------------------------------------------------------------------------*/
137/*---------------------------------------------------------------------------*/
138
139template <class Matrix, class Vector1, class Vector2, class Vector3>
140struct residual_impl<Matrix, Vector1, Vector2, Vector3,
141 std::enable_if_t<detail::use_builtin_matrix_ops<Matrix>::value &&
142 (math::static_rows<typename value_type<Matrix>::type>::value != math::static_rows<typename value_type<Vector1>::type>::value ||
143 math::static_rows<typename value_type<Matrix>::type>::value != math::static_rows<typename value_type<Vector2>::type>::value ||
144 math::static_rows<typename value_type<Matrix>::type>::value != math::static_rows<typename value_type<Vector3>::type>::value)>>
145{
146 static void apply(Vector1 const& f,
147 Matrix const& A,
148 Vector2 const& x,
149 Vector3& r)
150 {
151 typedef typename value_type<Matrix>::type V;
152
153 auto X = backend::reinterpret_as_rhs<V>(x);
154 auto F = backend::reinterpret_as_rhs<V>(f);
155 auto R = backend::reinterpret_as_rhs<V>(r);
156
157 residual(F, A, X, R);
158 }
159};
160
161/*---------------------------------------------------------------------------*/
162/*---------------------------------------------------------------------------*/
163
164} // namespace Arcane::Alina::backend
165
166/*---------------------------------------------------------------------------*/
167/*---------------------------------------------------------------------------*/
168
169#endif
Arcane::ForLoopRunInfo
Loop execution information.
Definition ForLoopRunInfo.h:38
Arcane::Matrix
Matrix class, to be used by user.
Definition matrix/Matrix.h:36
Arcane::Vector2
Class managing a 2-dimensional vector of type T.
Definition Vector2.h:37
Arcane::Vector3
Class managing a 3-dimensional vector of type T.
Definition Vector3.h:37
Arcane::arccoreParallelFor
void arccoreParallelFor(const ComplexForLoopRanges< RankValue > &loop_ranges, const ForLoopRunInfo &run_info, const LambdaType &lambda_function, const ReducerArgs &... reducer_args)
Applies the lambda function lambda_function concurrently over the iteration interval given by loop_ra...
Definition ParallelFor.h:87
Arcane::Int32
std::int32_t Int32
Signed integer type of 32 bits.
Definition ArccoreGlobal.h:233
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