Alien  1.3.0
Developer documentation
Loading...
Searching...
No Matches
CG.h
1// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
2//-----------------------------------------------------------------------------
3// Copyright 2000-2026 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
4// See the top-level COPYRIGHT file for details.
5// SPDX-License-Identifier: Apache-2.0
6//-----------------------------------------------------------------------------
7
8
9#pragma once
10#include <ostream>
11#include <vector>
12
13namespace Alien
14{
15
16template <typename AlgebraT>
17class CG
18{
19 public:
20 // clang-format off
21 typedef AlgebraT AlgebraType;
22 typedef typename AlgebraType::Matrix MatrixType;
23 typedef typename AlgebraType::Vector VectorType;
24 typedef typename MatrixType::ValueType ValueType;
25 typedef typename AlgebraType::FutureType FutureType;
26 // clang-format on
27
28 CG(AlgebraType& algebra, ITraceMng* trace_mng = nullptr)
29 : m_algebra(algebra)
30 , m_trace_mng(trace_mng)
31 {}
32
33 virtual ~CG()
34 {}
35
36 void setOutputLevel(int level)
37 {
38 m_output_level = level;
39 }
40
41 template <typename PrecondT, typename iterT>
42 int solve(PrecondT& precond,
43 iterT& iter,
44 MatrixType const& A,
45 VectorType const& b,
46 VectorType& x)
47 {
48 if (iter.nullRhs())
49 return 0;
50 ValueType rho(0), rho1(0), alpha(0);
51 VectorType p, z, q, r;
52
53 m_algebra.allocate(AlgebraType::resource(A), p, z, q, r);
54
55 // SEQ0
56 // r = b - A * x;
57 m_algebra.copy(b, r);
58 m_algebra.mult(A, x, p);
59 m_algebra.axpy(-1., p, r);
60
61 // SEQ1
62 /*
63 * z = solve(M,r)
64 * p = z
65 * q = A p
66 * rho1 = dot(r,z)
67 * alpha = dot(p,q)
68 * alpha = rho1/alpha
69 * x += alpha*p
70 * r += alpha*q
71 */
72 m_algebra.exec(precond, r, z);
73 rho1 = m_algebra.dot(r, z);
74 m_algebra.copy(z, p);
75 m_algebra.mult(A, p, q);
76 alpha = m_algebra.dot(p, q);
77 if (m_output_level > 2)
78 _print(0, "Seq 1", "rho", rho1, "alpha", alpha);
79 if (alpha == 0) {
80 if (iter.stop(r)) {
81 ++iter;
82 m_algebra.free(p, z, q, r);
83 return 0;
84 }
85 else
86 throw typename AlgebraType::NullValueException("alpha");
87 }
88 alpha = rho1 / alpha;
89 m_algebra.axpy(alpha, p, x);
90 m_algebra.axpy(-alpha, q, r);
91 rho = rho1;
92 ++iter;
93
94 while (!iter.stop(r)) {
95 // SEQ2
96 /*
97 * z = solve(M,r)
98 * rho1 = dot(r,z)
99 * alpha = rho1/rho
100 * p = z + alpha*p
101 * q = A*p
102 * alpha = dot(q,p)
103 * alpha = rho1/alpha
104 * x += alpha*p
105 * r -= alpha*q
106 */
107 m_algebra.exec(precond, r, z);
108 rho1 = m_algebra.dot(r, z);
109 alpha = rho1 / rho;
110 m_algebra.axpy(alpha, p, z);
111 m_algebra.copy(z, p);
112 m_algebra.mult(A, p, q);
113 alpha = m_algebra.dot(q, p);
114 if (alpha == 0) {
115 if (iter.stop(r)) {
116 ++iter;
117 m_algebra.free(p, z, q, r);
118 return 0;
119 }
120 else
121 throw typename AlgebraType::NullValueException("alpha");
122 }
123 alpha = rho1 / alpha;
124 m_algebra.axpy(alpha, p, x);
125 m_algebra.axpy(-alpha, q, r);
126 rho = rho1;
127 ++iter;
128 }
129
130 m_algebra.free(p, z, q, r);
131
132 return 0;
133 }
134
135 template <typename PrecondT, typename iterT>
136 int solve2(PrecondT& precond,
137 iterT& iter,
138 MatrixType const& A,
139 VectorType const& b,
140 VectorType& x)
141 {
142
143 if (iter.nullRhs())
144 return 0;
145 ValueType rho(0), rho1(0), alpha(0);
146 FutureType frho(rho), frho1(rho1), falpha(alpha);
147 VectorType p, z, q, r;
148
149 m_algebra.allocate(AlgebraType::resource(A), p, z, q, r);
150
151 // SEQ0
152 // r = b - A * x;
153 m_algebra.copy(b, r);
154 m_algebra.mult(A, x, p);
155 m_algebra.axpy(-1., p, r);
156
157 // SEQ1
158 /*
159 * z = solve(M,r)
160 * p = z
161 * q = A p
162 * rho1 = dot(r,p)
163 * alpha = dot(p,q)
164 * alpha = rho1/alpha
165 * x += alpha*p
166 * r += alpha*q
167 */
168 m_algebra.exec(precond, r, z);
169 m_algebra.copy(z, p);
170 m_algebra.mult(A, p, q);
171 m_algebra.dot(r, p, frho1);
172 m_algebra.dot(p, q, falpha);
173 if (falpha.get() == 0) {
174 if (iter.stop(r)) {
175 ++iter;
176 m_algebra.free(p, z, q, r);
177 return 0;
178 }
179 else
180 throw typename AlgebraType::NullValueException("alpha");
181 }
182 alpha = frho1.get() / alpha;
183 m_algebra.axpy(alpha, p, x);
184 m_algebra.axpy(-alpha, q, r);
185 rho = rho1;
186 ++iter;
187
188 while (!iter.stop(r)) {
189
190 /*
191 * z = solve(M,r)
192 * rho1 = dot(r,z)
193 * alpha = rho1/rho
194 * p = z + alpha*p
195 * q = A*p
196 * alpha = dot(q,p)
197 * alpha = rho1/alpha
198 * x += alpha*p
199 * r -= alpha*q
200 */
201 m_algebra.exec(precond, r, z);
202 m_algebra.dot(r, z, frho1);
203 alpha = frho1.get() / rho;
204 m_algebra.axpy(alpha, p, z);
205 m_algebra.copy(z, p);
206 m_algebra.mult(A, p, q);
207 m_algebra.dot(p, q, falpha);
208 if (falpha.get() == 0) {
209 if (iter.stop(r)) {
210 ++iter;
211 m_algebra.free(p, z, q, r);
212 return 0;
213 }
214 else
215 throw typename AlgebraType::NullValueException("alpha");
216 }
217 alpha = rho1 / alpha;
218 m_algebra.axpy(alpha, p, x);
219 m_algebra.axpy(-alpha, q, r);
220 rho = rho1;
221 ++iter;
222 }
223
224 m_algebra.free(p, z, q, r);
225
226 return 0;
227 }
228
229 private:
230 void
231 _print(int iter, std::string const& msg, std::string const& label0,
232 ValueType value0)
233 {
234 if (m_trace_mng) {
235 m_trace_mng->info() << msg;
236 m_trace_mng->info() << "Iterate: " << iter << " " << label0 << " "
237 << value0;
238 }
239 }
240
241 void
242 _print(int iter, std::string const& msg, std::string const& label0,
243 ValueType value0, std::string const& label1, ValueType value1)
244 {
245 if (m_trace_mng) {
246 _print(iter, msg, label0, value0);
247 m_trace_mng->info() << "Iterate: " << iter << " " << label1 << " "
248 << value1;
249 }
250 }
251
252 void
253 _print(int iter, std::string const& msg, std::string const& label0,
254 ValueType value0, std::string const& label1, ValueType value1,
255 std::string const& label2, ValueType value2)
256 {
257 if (m_trace_mng) {
258 _print(iter, msg, label0, value0, label1, value1);
259 m_trace_mng->info() << "Iterate: " << iter << " " << label2 << " "
260 << value2;
261 }
262 }
263
264 void
265 _print(int iter, std::string const& msg, std::string const& label0,
266 ValueType value0, std::string const& label1, ValueType value1,
267 std::string const& label2, ValueType value2,
268 std::string const& label3, ValueType value3)
269 {
270 if (m_trace_mng) {
271 _print(iter, msg, label0, value0, label1, value1, label2, value2);
272 m_trace_mng->info() << "Iterate: " << iter << " " << label3 << " "
273 << value3;
274 }
275 }
276
277 AlgebraType& m_algebra;
278 ITraceMng* m_trace_mng = nullptr;
279 int m_output_level = 0;
280};
281} // namespace Alien
Alien
-- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature --
Definition BackEnd.h:17