Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/LinearSolver/ConjugateGradients.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2025 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief Charm++ chare array for distributed conjugate gradients.
10 : : \details Charm++ chare array for asynchronous distributed
11 : : conjugate gradients linear solver for the symmetric linear system
12 : : A * x = b, where A is in compressed sparse row storage, containing
13 : : only the upper triangular part of A.
14 : : \see Y. Saad, Iterative Methods for Sparse Linear Systems: 2nd Edition,
15 : : ISBN 9780898718003, 2003.
16 : :
17 : : Algorithm 'Preconditioned Conjugate Gradient':
18 : : ```
19 : : Compute r0 := b - A x0, z0 := M^{-1} r0, p0 := z0
20 : : For j=0,1,..., until convergence, do
21 : : alpha_j := (r_j,z_j) / (A p_j,p_j)
22 : : x_{j+1} := x_j + alpha_j p_j
23 : : r_{j+1} := r_j - alpha_j A p_j
24 : : z_{j+1} := M^{-1} r_{j+1}
25 : : beta_j := (r_{j+1},z_{j+1}) / (r_j,z_j)
26 : : p_{j+1} := z_{j+1} + beta_j p_j
27 : : end
28 : : ```
29 : : */
30 : : // *****************************************************************************
31 : :
32 : : #include <numeric>
33 : :
34 : : #include "Exception.hpp"
35 : : #include "ConjugateGradients.hpp"
36 : : #include "Vector.hpp"
37 : : #include "ContainerUtil.hpp"
38 : : #include "Reorder.hpp"
39 : : #include "Print.hpp"
40 : :
41 : : using tk::ConjugateGradients;
42 : :
43 : 810 : ConjugateGradients::ConjugateGradients(
44 : : const CSR& A,
45 : : const std::vector< tk::real >& x,
46 : : const std::vector< tk::real >& b,
47 : : const std::vector< std::size_t >& gid,
48 : : const std::unordered_map< std::size_t, std::size_t >& lid,
49 : : const std::unordered_map< int,
50 : 810 : std::unordered_set< std::size_t > >& nodecommmap ) :
51 [ + - ]: 810 : m_A( A ),
52 [ + - ]: 810 : m_An( A ),
53 [ + - ]: 810 : m_x( x ),
54 [ + - ]: 810 : m_b( b ),
55 [ + - ][ - - ]: 810 : m_pc( "none" ),
56 [ + - ]: 810 : m_gid( gid ),
57 : : m_lid( lid ),
58 : : m_nodeCommMap( nodecommmap ),
59 [ + - ][ + - ]: 810 : m_r( m_A.rsize(), 0.0 ),
60 [ + - ]: 810 : m_z( m_A.rsize(), 0.0 ),
61 [ + - ][ - - ]: 810 : m_d( m_A.rsize(), 0.0 ),
62 : 810 : m_nr( 0 ),
63 [ + - ]: 810 : m_na( 0 ),
64 : 810 : m_nb( 0 ),
65 [ + - ]: 810 : m_p( m_A.rsize(), 0.0 ),
66 [ + - ][ + + ]: 810 : m_q( m_A.rsize(), 0.0 ),
[ - - ]
67 : 810 : m_nq( 0 ),
68 [ + + ]: 810 : m_nd( 0 ),
69 : : m_initres(),
70 : : m_solved(),
71 : 810 : m_normb( 0.0 ),
72 : 810 : m_it( 0 ),
73 : 810 : m_maxit( 0 ),
74 : 810 : m_finished( false ),
75 : 810 : m_verbose( 0 ),
76 : 810 : m_rho( 0.0 ),
77 : 810 : m_rho0( 0.0 ),
78 : 810 : m_alpha( 0.0 ),
79 : 810 : m_converged( false ),
80 : 810 : m_nx( 0 ),
81 [ + - ][ + + ]: 1620 : m_apply( false )
82 : : // *****************************************************************************
83 : : // Constructor
84 : : //! \param[in] A Left hand side matrix of the linear system to solve in Ax=b
85 : : //! \param[in] x Solution (initial guess) of the linear system to solve in Ax=b
86 : : //! \param[in] b Right hand side of the linear system to solve in Ax=b
87 : : //! \param[in] gid Global node ids
88 : : //! \param[in] lid Local node ids associated to global ones
89 : : //! \param[in] nodecommmap Global mesh node IDs shared with other chares
90 : : //! associated to their chare IDs
91 : : // *****************************************************************************
92 : : {
93 : : // Fill in gid and lid for serial solve
94 [ + + ][ + - ]: 810 : if (gid.empty() || lid.empty() || nodecommmap.empty()) {
[ + + ]
95 [ + - ]: 27 : m_gid.resize( m_A.rsize()/m_A.Ncomp() );
96 : : std::iota( begin(m_gid), end(m_gid), 0 );
97 [ + - ][ + + ]: 20763 : for (auto g : m_gid) m_lid[g] = g;
98 : : }
99 : :
100 : : Assert( m_A.rsize() == m_gid.size()*A.Ncomp(), "Size mismatch" );
101 : : Assert( m_x.size() == m_gid.size()*A.Ncomp(), "Size mismatch" );
102 : : Assert( m_b.size() == m_gid.size()*A.Ncomp(), "Size mismatch" );
103 : 810 : }
104 : :
105 : : void
106 : 6288 : ConjugateGradients::setup( CkCallback c )
107 : : // *****************************************************************************
108 : : // Setup solver
109 : : //! \param[in] c Call to continue with
110 : : //! \details This function initiates computing the residual (r=b-A*x), its dot
111 : : //! product, and the rhs norm.
112 : : // *****************************************************************************
113 : : {
114 : 6288 : m_initres = c;
115 : 6288 : m_converged = false;
116 : 6288 : m_finished = false;
117 : :
118 : : // initiate computing A * x (for the initial residual)
119 [ + - ]: 6288 : thisProxy[ thisIndex ].wait4res();
120 : 6288 : residual();
121 : 6288 : pc();
122 : :
123 : : // initiate computing norm of right hand side
124 [ + - ]: 6288 : dot( m_b, m_b,
125 : 6288 : CkCallback( CkReductionTarget(ConjugateGradients,normb), thisProxy ) );
126 : 6288 : }
127 : :
128 : : void
129 : 719715 : ConjugateGradients::dot( const std::vector< tk::real >& a,
130 : : const std::vector< tk::real >& b,
131 : : CkCallback c )
132 : : // *****************************************************************************
133 : : // Initiate computation of dot product of two vectors
134 : : //! \param[in] a 1st vector of dot product
135 : : //! \param[in] b 2nd vector of dot product
136 : : //! \param[in] c Callback to target with the final result
137 : : // *****************************************************************************
138 : : {
139 : : Assert( a.size() == b.size(), "Size mismatch" );
140 : :
141 : 719715 : tk::real D = 0.0;
142 : : auto ncomp = m_A.Ncomp();
143 [ + + ]: 164109892 : for (std::size_t i=0; i<a.size()/ncomp; ++i) {
144 : 163390177 : auto incomp = i*ncomp;
145 [ + + ]: 163390177 : if (not slave(m_nodeCommMap,m_gid[i],thisIndex))
146 [ + + ]: 306286276 : for (std::size_t d=0; d<ncomp; ++d)
147 : 153552142 : D += a[incomp+d] * b[incomp+d];
148 : : }
149 : :
150 : 719715 : contribute( sizeof(tk::real), &D, CkReduction::sum_double, c );
151 : 719715 : }
152 : :
153 : : void
154 : 6288 : ConjugateGradients::normb( tk::real n )
155 : : // *****************************************************************************
156 : : // Compute the norm of the right hand side
157 : : //! \param[in] n Norm of right hand side (aggregated across all chares)
158 : : // *****************************************************************************
159 : : {
160 : 6288 : m_normb = std::sqrt(n);
161 : 6288 : normb_complete();
162 : 6288 : }
163 : :
164 : : void
165 : 6288 : ConjugateGradients::residual()
166 : : // *****************************************************************************
167 : : // Initiate A * x for computing the initial residual, r = b - A * x
168 : : // *****************************************************************************
169 : : {
170 : : // Compute own contribution to r = A * x
171 : 6288 : m_A.mult( m_x, m_r );
172 : :
173 : : // Send partial product on chare-boundary nodes to fellow chares
174 [ + + ]: 6288 : if (m_nodeCommMap.empty()) {
175 : 301 : comres_complete();
176 : : } else {
177 : : auto ncomp = m_A.Ncomp();
178 [ + + ]: 62377 : for (const auto& [c,n] : m_nodeCommMap) {
179 : 56390 : std::vector< std::vector< tk::real > > rc( n.size() );
180 : : std::size_t j = 0;
181 [ + + ]: 317852 : for (auto g : n) {
182 [ + - ]: 261462 : std::vector< tk::real > nr( ncomp );
183 : 261462 : auto i = tk::cref_find( m_lid, g );
184 [ + + ]: 523920 : for (std::size_t d=0; d<ncomp; ++d) nr[d] = m_r[ i*ncomp+d ];
185 : 261462 : rc[j++] = std::move(nr);
186 : : }
187 [ + - ][ + - ]: 112780 : thisProxy[c].comres( std::vector<std::size_t>(begin(n),end(n)), rc );
[ + - ][ - + ]
[ - - ]
188 : 56390 : }
189 : : }
190 : :
191 : 6288 : ownres_complete();
192 : 6288 : }
193 : :
194 : : void
195 : 56390 : ConjugateGradients::comres( const std::vector< std::size_t >& gid,
196 : : const std::vector< std::vector< tk::real > >& rc )
197 : : // *****************************************************************************
198 : : // Receive contributions to A * x on chare-boundaries
199 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
200 : : //! \param[in] rc Partial contributions at chare-boundary nodes
201 : : // *****************************************************************************
202 : : {
203 : : Assert( rc.size() == gid.size(), "Size mismatch" );
204 : :
205 [ + + ]: 317852 : for (std::size_t i=0; i<gid.size(); ++i) m_rc[ gid[i] ] += rc[i];
206 : :
207 [ + + ]: 56390 : if (++m_nr == m_nodeCommMap.size()) {
208 : 5987 : m_nr = 0;
209 : 5987 : comres_complete();
210 : : }
211 : 56390 : }
212 : :
213 : : void
214 [ + + ]: 6288 : ConjugateGradients::pc()
215 : : // *****************************************************************************
216 : : // Setup preconditioner
217 : : // *****************************************************************************
218 : : {
219 : : auto ncomp = m_A.Ncomp();
220 : :
221 [ + + ]: 6288 : if (m_pc == "none") {
222 : :
223 [ + + ]: 113781 : for (std::size_t i=0; i<m_q.size()/ncomp; ++i) {
224 : 109659 : auto c = tk::count( m_nodeCommMap, m_gid[i] );
225 [ + + ]: 219392 : for (std::size_t d=0; d<ncomp; ++d) {
226 : 109733 : m_q[i*ncomp+d] = 1.0 / c;
227 : : }
228 : : }
229 : :
230 [ + - ]: 2166 : } else if (m_pc == "jacobi") {
231 : :
232 : : // Extract Jacobi preconditioner from matrix
233 [ + + ]: 537872 : for (std::size_t i=0; i<m_q.size()/ncomp; ++i) {
234 [ + + ]: 1120372 : for (std::size_t d=0; d<ncomp; ++d) {
235 : 584666 : m_q[i*ncomp+d] = m_A(i,i,d);
236 : : }
237 : : }
238 : :
239 : : }
240 : :
241 : : // Send partial preconditioner on chare-boundary nodes to fellow chares
242 [ + + ]: 6288 : if (m_nodeCommMap.empty()) {
243 : 301 : comd_complete();
244 : : } else {
245 [ + + ]: 62377 : for (const auto& [c,n] : m_nodeCommMap) {
246 : 56390 : std::vector< std::vector< tk::real > > qc( n.size() );
247 : : std::size_t j = 0;
248 [ + + ]: 317852 : for (auto g : n) {
249 [ + - ]: 261462 : std::vector< tk::real > nq( ncomp );
250 : 261462 : auto i = tk::cref_find( m_lid, g );
251 [ + + ]: 523920 : for (std::size_t d=0; d<ncomp; ++d) nq[d] = m_q[ i*ncomp+d ];
252 : 261462 : qc[j++] = std::move(nq);
253 : : }
254 [ + - ][ + - ]: 112780 : thisProxy[c].comd( std::vector<std::size_t>(begin(n),end(n)), qc );
[ + - ][ - + ]
[ - - ]
255 : 56390 : }
256 : : }
257 : :
258 : 6288 : ownd_complete();
259 : 6288 : }
260 : :
261 : : void
262 : 56390 : ConjugateGradients::comd( const std::vector< std::size_t >& gid,
263 : : const std::vector< std::vector< tk::real > >& qc )
264 : : // *****************************************************************************
265 : : // Receive contributions to preconditioner on chare-boundaries
266 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
267 : : //! \param[in] qc Partial contributions at chare-boundary nodes
268 : : // *****************************************************************************
269 : : {
270 : : Assert( qc.size() == gid.size(), "Size mismatch" );
271 : :
272 [ + + ]: 317852 : for (std::size_t i=0; i<gid.size(); ++i) m_qc[ gid[i] ] += qc[i];
273 : :
274 [ + + ]: 56390 : if (++m_nd == m_nodeCommMap.size()) {
275 : 5987 : m_nd = 0;
276 : 5987 : comd_complete();
277 : : }
278 : 56390 : }
279 : :
280 : : void
281 : 6288 : ConjugateGradients::initres()
282 : : // *****************************************************************************
283 : : // Finish computing the initial residual, r = b - A * x
284 : : // *****************************************************************************
285 : : {
286 : : auto ncomp = m_A.Ncomp();
287 : :
288 : : // Combine own and communicated contributions to r = A * x
289 [ + + ]: 139277 : for (const auto& [gid,r] : m_rc) {
290 : 132989 : auto i = tk::cref_find( m_lid, gid );
291 [ + + ]: 266974 : for (std::size_t c=0; c<ncomp; ++c) m_r[i*ncomp+c] += r[c];
292 : : }
293 : 6288 : tk::destroy( m_rc );
294 : :
295 : : // Finish computing the initial residual, r = b - A * x
296 [ + + ]: 700687 : for (auto& r : m_r) r *= -1.0;
297 : 6288 : m_r += m_b;
298 : :
299 : : // Initialize p
300 : 6288 : m_p = m_r;
301 : :
302 : : // Combine own and communicated contributions to q = A * p
303 [ + + ]: 144518 : for (const auto& [gid,q] : m_qc) {
304 : 138230 : auto i = tk::cref_find( m_lid, gid );
305 [ + + ]: 277456 : for (std::size_t c=0; c<ncomp; ++c)
306 : 139226 : m_q[i*ncomp+c] += q[c];
307 : : }
308 : 6288 : tk::destroy( m_qc );
309 : :
310 : : // Set preconditioner
311 : 6288 : m_d = m_q;
312 : :
313 : : // initialize preconditioned solve: compute z = M^{-1} * r
314 [ + + ]: 700687 : for (std::size_t i=0; i<m_z.size(); ++i) m_z[i] = m_r[i] / m_d[i];
315 : :
316 : : // initiate computing the dot product of the initial residual, rho = (r,z)
317 [ + - ]: 6288 : dot( m_r, m_z,
318 : 6288 : CkCallback( CkReductionTarget(ConjugateGradients,rho), thisProxy ) );
319 : 6288 : }
320 : :
321 : : void
322 : 6288 : ConjugateGradients::rho( tk::real r )
323 : : // *****************************************************************************
324 : : // Compute rho = (r,z)
325 : : //! \param[in] r Dot product, rho = (r,z) (aggregated across all chares)
326 : : // *****************************************************************************
327 : : {
328 : : // store dot product of residual
329 : 6288 : m_rho = r;
330 : :
331 : : // send back rhs norm to caller
332 : 6288 : m_initres.send( CkDataMsg::buildNew( sizeof(tk::real), &m_normb ) );
333 : 6288 : }
334 : :
335 : : void
336 : 6282 : ConjugateGradients::init(
337 : : const std::vector< tk::real >& x,
338 : : const std::vector< tk::real >& b,
339 : : const std::vector< tk::real >& neubc,
340 : : const std::unordered_map< std::size_t,
341 : : std::vector< std::pair< int, tk::real > > >& dirbc,
342 : : bool apply,
343 : : const std::string& pc,
344 : : CkCallback cb )
345 : : // *****************************************************************************
346 : : // Initialize linear solve: set initial guess and boundary conditions
347 : : //! \param[in] x Initial guess
348 : : //! \param[in] b Right hand side vector
349 : : //! \param[in] neubc Right hand side vector with Neumann BCs partially applied
350 : : //! \param[in] dirbc Local node ids and associated Dirichlet BCs
351 : : //! \param[in] apply True to apply boundary conditions
352 : : //! \param[in] pc Preconditioner to use
353 : : //! \param[in] cb Call to continue with when initialized and ready for a solve
354 : : //! \details This function allows setting the initial guess and boundary
355 : : //! conditions, followed by computing the initial residual and the rhs norm.
356 : : //! It also performs communication of BC data.
357 : : // *****************************************************************************
358 : : {
359 : : // Configure preconditioner
360 : 6282 : m_pc = pc;
361 : :
362 : : // Optionally set initial guess
363 [ - + ]: 6282 : if (!x.empty()) m_x = x; //else std::fill( begin(m_x), end(m_x), 0.0 );
364 : :
365 : : // Optionally set rhs, assumed complete in parallel
366 [ + - ]: 6282 : if (!b.empty()) m_b = b;
367 : :
368 : : // Save if applied BCs
369 : 6282 : m_apply = apply;
370 [ - + ]: 6282 : if (not m_apply) {
371 : :
372 : : // Recompute initial residual (r=b-A*x), its dot product, and the rhs norm
373 [ - - ]: 0 : setup( cb );
374 : :
375 : : }
376 : : else {
377 : :
378 : : // Save matrix state
379 : 6282 : m_An = m_A;
380 : :
381 : : // Set Neumann BCs, partial in parallel, communication below
382 [ + + ]: 6282 : if (!neubc.empty()) m_q = neubc; else std::fill(begin(m_q), end(m_q), 0.0);
383 : :
384 : : // Store incoming Dirichlet BCs, partial in parallel, communication below
385 : 6282 : tk::destroy(m_dirbc);
386 [ + + ]: 40766 : for (auto&& [i,bcval] : dirbc) m_dirbc[i] = std::move(bcval);
387 : :
388 : : // Get ready to communicate boundary conditions. This is necessary because
389 : : // there can be nodes a chare contributes to but does not apply BCs on.
390 : : // This happens if a node is in the node communication map but not on the
391 : : // list of incoming BCs on this chare. To have all chares share the same
392 : : // view on all BC nodes, we send the global node ids together with the
393 : : // Dirichlet BCs at which BCs are set to those fellow chares that also
394 : : // contribute to those BC nodes. Only after this communication step will we
395 : : // apply the BCs, which then will correctly setup the BC rows of the matrix
396 : : // and on the rhs that (1) may contain a nonzero value, (2) may have partial
397 : : // contributions due to Neumann BCs, and (3) will be modified to apply
398 : : // Dirichlet BCs.
399 [ + - ]: 6282 : thisProxy[ thisIndex ].wait4bc();
400 [ + - ][ + + ]: 12564 : thisProxy[ thisIndex ].wait4r();
401 : :
402 : : // Send boundary conditions to those who contribute to those rows
403 [ + + ]: 6282 : if (m_nodeCommMap.empty()) {
404 : 299 : combc_complete();
405 : : } else {
406 : : auto ncomp = m_A.Ncomp();
407 [ + + ]: 62369 : for (const auto& [c,n] : m_nodeCommMap) {
408 : : decltype(m_dirbc) dbc;
409 [ + - ]: 56386 : std::vector< std::vector< tk::real > > qc( n.size() );
410 : : std::size_t k = 0;
411 [ + + ]: 317812 : for (auto g : n) {
412 : 261426 : auto i = tk::cref_find( m_lid, g );
413 : : auto j = m_dirbc.find(i);
414 [ + + ][ + - ]: 261426 : if (j != end(m_dirbc)) dbc[g] = j->second;
[ + - ]
415 [ + - ]: 261426 : std::vector< tk::real > nq( ncomp );
416 [ + + ]: 523812 : for (std::size_t d=0; d<ncomp; ++d) nq[d] = m_q[ i*ncomp+d ];
417 : 261426 : qc[k++] = std::move(nq);
418 : : }
419 [ + - ]: 56386 : std::vector< std::size_t > gid( begin(n), end(n) );
420 [ + - ][ + - ]: 112772 : thisProxy[c].combc( dbc, gid, qc );
[ + - ][ - - ]
421 : 56386 : }
422 : : }
423 : :
424 [ + - ]: 12564 : ownbc_complete( cb );
425 : :
426 : : }
427 : 6282 : }
428 : :
429 : : void
430 : 56386 : ConjugateGradients::combc(
431 : : const std::map< std::size_t, std::vector< std::pair< int, tk::real > > >& dbc,
432 : : const std::vector< std::size_t >& gid,
433 : : const std::vector< std::vector< tk::real > >& qc )
434 : : // *****************************************************************************
435 : : // Receive contributions to boundary conditions and rhs on chare-boundaries
436 : : //! \param[in] dbc Contributions to Dirichlet boundary conditions
437 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
438 : : //! \param[in] qc Contributions to Neumann boundary conditions
439 : : // *****************************************************************************
440 : : {
441 [ + + ]: 57638 : for (const auto& [g,dirbc] : dbc) m_dirbcc[g] = dirbc;
442 [ + + ]: 317812 : for (std::size_t i=0; i<gid.size(); ++i) m_qc[ gid[i] ] += qc[i];
443 : :
444 [ + + ]: 56386 : if (++m_nb == m_nodeCommMap.size()) {
445 : 5983 : m_nb = 0;
446 : 5983 : combc_complete();
447 : : }
448 : 56386 : }
449 : :
450 : : void
451 : 6282 : ConjugateGradients::apply( CkCallback cb )
452 : : // *****************************************************************************
453 : : // Apply boundary conditions
454 : : //! \param[in] cb Call to continue with after applying the BCs is complete
455 : : // *****************************************************************************
456 : : {
457 : : auto ncomp = m_A.Ncomp();
458 : :
459 : : // Merge own and received contributions to Dirichlet boundary conditions
460 [ + + ]: 7355 : for (const auto& [g,dirbc] : m_dirbcc) {
461 : 1073 : m_dirbc[ tk::cref_find(m_lid,g) ] = dirbc;
462 : : }
463 : 6282 : tk::destroy( m_dirbcc );
464 : :
465 : : // Merge own and received contributions to Neumann boundary conditions
466 [ + + ]: 144476 : for (const auto& [gid,q] : m_qc) {
467 : 138194 : auto i = tk::cref_find( m_lid, gid );
468 [ + + ]: 277348 : for (std::size_t c=0; c<ncomp; ++c) m_q[i*ncomp+c] += q[c];
469 : : }
470 : 6282 : tk::destroy( m_qc );
471 : :
472 : : // Apply Neumann BCs on rhs
473 : 6282 : m_b += m_q;
474 : :
475 : : // Apply Dirichlet BCs on matrix and rhs (with decreasing local id)
476 : : std::fill( begin(m_r), end(m_r), 0.0 );
477 [ + + ]: 40766 : for (auto bi = m_dirbc.rbegin(); bi != m_dirbc.rend(); ++bi) {
478 : 34484 : auto i = bi->first;
479 : : const auto& dirbc = bi->second;
480 [ + + ]: 117928 : for (std::size_t j=0; j<ncomp; ++j) {
481 [ + + ]: 83444 : if (dirbc[j].first) {
482 : 67124 : m_A.dirichlet( i, dirbc[j].second, m_r, m_gid, m_nodeCommMap, j );
483 : : }
484 : : }
485 : : }
486 : :
487 : : // Send rhs with Dirichlet BCs partially applied to fellow chares
488 [ + + ]: 6282 : if (m_nodeCommMap.empty()) {
489 : 299 : comr_complete();
490 : : } else {
491 [ + + ]: 62369 : for (const auto& [c,n] : m_nodeCommMap) {
492 : 56386 : std::vector< std::vector< tk::real > > rc( n.size() );
493 : : std::size_t j = 0;
494 [ + + ]: 317812 : for (auto g : n) {
495 [ + - ]: 261426 : std::vector< tk::real > nr( ncomp );
496 : 261426 : auto i = tk::cref_find( m_lid, g );
497 [ + + ]: 523812 : for (std::size_t d=0; d<ncomp; ++d) nr[d] = m_r[ i*ncomp+d ];
498 : 261426 : rc[j++] = std::move(nr);
499 : : }
500 [ + - ][ + - ]: 112772 : thisProxy[c].comr( std::vector<std::size_t>(begin(n),end(n)), rc );
[ + - ][ - + ]
[ - - ]
501 : 56386 : }
502 : : }
503 : :
504 [ + - ]: 6282 : ownr_complete( cb );
505 : 6282 : }
506 : :
507 : : void
508 : 56386 : ConjugateGradients::comr( const std::vector< std::size_t >& gid,
509 : : const std::vector< std::vector< tk::real > >& rc )
510 : : // *****************************************************************************
511 : : // Receive contributions to rhs with Dirichlet BCs applied on chare-boundaries
512 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
513 : : //! \param[in] rc Partial contributions at chare-boundary nodes
514 : : // *****************************************************************************
515 : : {
516 : : Assert( rc.size() == gid.size(), "Size mismatch" );
517 : :
518 [ + + ]: 317812 : for (std::size_t i=0; i<gid.size(); ++i) m_rc[ gid[i] ] += rc[i];
519 : :
520 [ + + ]: 56386 : if (++m_na == m_nodeCommMap.size()) {
521 : 5983 : m_na = 0;
522 : 5983 : comr_complete();
523 : : }
524 : 56386 : }
525 : :
526 : : void
527 : 6282 : ConjugateGradients::r( CkCallback cb )
528 : : // *****************************************************************************
529 : : // Finish computing rhs with BCs applied
530 : : //! \param[in] cb Call to continue with after applying the BCs is complete
531 : : // *****************************************************************************
532 : : {
533 : : auto ncomp = m_A.Ncomp();
534 : :
535 : : // Combine own and communicated contributions to Dirichlet BCs applied to rhs
536 [ + + ]: 144476 : for (const auto& [gid,r] : m_rc) {
537 : 138194 : auto i = tk::cref_find( m_lid, gid );
538 [ + + ]: 277348 : for (std::size_t c=0; c<ncomp; ++c) m_r[i*ncomp+c] += r[c];
539 : : }
540 : 6282 : tk::destroy( m_rc );
541 : :
542 : : // Subtract matrix columns * BC values from rhs at Dirichlet BC nodes
543 : 6282 : m_b -= m_r;
544 : :
545 : : // Apply Dirichlet BCs on rhs
546 [ + + ]: 40766 : for (const auto& [i,dirbc] : m_dirbc) {
547 [ + + ]: 117928 : for (std::size_t j=0; j<ncomp; ++j) {
548 [ + + ]: 83444 : if (dirbc[j].first) {
549 : 67124 : m_b[i*ncomp+j] = dirbc[j].second;
550 : : }
551 : : }
552 : : }
553 : :
554 : : // Recompute initial residual (r=b-A*x), its dot product, and the rhs norm
555 [ + - ]: 6282 : setup( cb );
556 : 6282 : }
557 : :
558 : : void
559 : 6288 : ConjugateGradients::solve( std::size_t maxit,
560 : : tk::real tol,
561 : : int pe,
562 : : uint64_t verbose,
563 : : CkCallback c )
564 : : // *****************************************************************************
565 : : // Solve linear system
566 : : //! \param[in] maxit Max iteration count
567 : : //! \param[in] tol Stop tolerance
568 : : //! \param[in] pe Processing element
569 : : //! \param[in] verbose Verbose output
570 : : //! \param[in] c Call to continue with after solve is complete
571 : : // *****************************************************************************
572 : : {
573 : 6288 : m_maxit = maxit;
574 : 6288 : m_tol = tol;
575 : 6288 : m_solved = c;
576 : 6288 : m_it = 0;
577 [ + + ]: 6288 : m_verbose = pe == 0 ? verbose : 0;
578 : :
579 [ + + ]: 6288 : if (m_verbose > 1) tk::Print() << "Xyst> Conjugate gradients start\n";
580 : :
581 [ - + ]: 6288 : if (m_converged) x(); else next();
582 : 6288 : }
583 : :
584 : : void
585 : 235713 : ConjugateGradients::next()
586 : : // *****************************************************************************
587 : : // Start next linear solver iteration
588 : : // *****************************************************************************
589 : : {
590 [ + + ]: 235713 : if (m_it == 0) m_alpha = 0.0; else m_alpha = m_rho/m_rho0;
591 : 235713 : m_rho0 = m_rho;
592 : :
593 : : // compute p = z + alpha * p
594 [ + + ]: 54511658 : for (std::size_t i=0; i<m_p.size(); ++i) m_p[i] = m_z[i] + m_alpha * m_p[i];
595 : :
596 : : // initiate computing q = A * p
597 [ + - ]: 235713 : thisProxy[ thisIndex ].wait4q();
598 : 235713 : qAp();
599 : 235713 : }
600 : :
601 : :
602 : : void
603 : 235713 : ConjugateGradients::qAp()
604 : : // *****************************************************************************
605 : : // Initiate computing q = A * p
606 : : // *****************************************************************************
607 : : {
608 : : // Compute own contribution to q = A * p
609 : 235713 : m_A.mult( m_p, m_q );
610 : :
611 : : // Send partial product on chare-boundary nodes to fellow chares
612 [ + + ]: 235713 : if (m_nodeCommMap.empty()) {
613 : 24113 : comq_complete();
614 : : } else {
615 : : auto ncomp = m_A.Ncomp();
616 [ + + ]: 1061348 : for (const auto& [c,n] : m_nodeCommMap) {
617 : 849748 : std::vector< std::vector< tk::real > > qc( n.size() );
618 : : std::size_t j = 0;
619 [ + + ]: 8869844 : for (auto g : n) {
620 [ + - ]: 8020096 : std::vector< tk::real > nq( ncomp );
621 : 8020096 : auto i = tk::cref_find( m_lid, g );
622 [ + + ]: 16045160 : for (std::size_t d=0; d<ncomp; ++d) nq[d] = m_q[ i*ncomp+d ];
623 : 8020096 : qc[j++] = std::move(nq);
624 : : }
625 [ + - ][ + - ]: 1699496 : thisProxy[c].comq( std::vector<std::size_t>(begin(n),end(n)), qc );
[ + - ][ - + ]
[ - - ]
626 : 849748 : }
627 : : }
628 : :
629 : 235713 : ownq_complete();
630 : 235713 : }
631 : :
632 : : void
633 : 849748 : ConjugateGradients::comq( const std::vector< std::size_t >& gid,
634 : : const std::vector< std::vector< tk::real > >& qc )
635 : : // *****************************************************************************
636 : : // Receive contributions to q = A * p on chare-boundaries
637 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
638 : : //! \param[in] qc Partial contributions at chare-boundary nodes
639 : : // *****************************************************************************
640 : : {
641 : : Assert( qc.size() == gid.size(), "Size mismatch" );
642 : :
643 [ + + ]: 8869844 : for (std::size_t i=0; i<gid.size(); ++i) m_qc[ gid[i] ] += qc[i];
644 : :
645 [ + + ]: 849748 : if (++m_nq == m_nodeCommMap.size()) {
646 : 211600 : m_nq = 0;
647 : 211600 : comq_complete();
648 : : }
649 : 849748 : }
650 : :
651 : : void
652 : 235713 : ConjugateGradients::q()
653 : : // *****************************************************************************
654 : : // Finish computing q = A * p
655 : : // *****************************************************************************
656 : : {
657 : : // Combine own and communicated contributions to q = A * p
658 : : auto ncomp = m_A.Ncomp();
659 [ + + ]: 6808865 : for (const auto& [gid,q] : m_qc) {
660 : 6573152 : auto i = tk::cref_find( m_lid, gid );
661 [ + + ]: 13151272 : for (std::size_t c=0; c<ncomp; ++c)
662 : 6578120 : m_q[i*ncomp+c] += q[c];
663 : : }
664 : 235713 : tk::destroy( m_qc );
665 : :
666 : : // initiate computing (p,q)
667 [ + - ]: 235713 : dot( m_p, m_q,
668 : 235713 : CkCallback( CkReductionTarget(ConjugateGradients,pq), thisProxy ) );
669 : 235713 : }
670 : :
671 : : void
672 [ + + ]: 235713 : ConjugateGradients::pq( tk::real d )
673 : : // *****************************************************************************
674 : : // Compute the dot product (p,q)
675 : : //! \param[in] d Dot product of (p,q) (aggregated across all chares)
676 : : // *****************************************************************************
677 : : {
678 : : // If (p,q)=0, then p and q are orthogonal and the system either has a trivial
679 : : // solution, x=x0, or the BCs are incomplete or wrong, in either case the
680 : : // solve cannot continue.
681 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
682 [ + + ]: 235713 : if (std::abs(d) < eps) {
683 [ - + ]: 4189 : if (m_verbose > 1) {
684 : : tk::Print() << "Xyst> Conjugate gradients it " << m_it << ", (p,q) = 0\n";
685 : : }
686 : 4189 : m_finished = true;
687 : 4189 : m_alpha = 0.0;
688 : : } else {
689 : 231524 : m_alpha = m_rho / d;
690 : : }
691 : :
692 : : // compute r = r - alpha * q
693 [ + + ]: 54511658 : for (std::size_t i=0; i<m_r.size(); ++i) m_r[i] -= m_alpha * m_q[i];
694 : :
695 : : // apply preconditioner: compute z = M^{-1} * r
696 [ + + ]: 54511658 : for (std::size_t i=0; i<m_z.size(); ++i) m_z[i] = m_r[i] / m_d[i];
697 : :
698 : : // initiate computing (r,z)
699 [ + - ]: 235713 : dot( m_r, m_z,
700 : 235713 : CkCallback( CkReductionTarget(ConjugateGradients,rz), thisProxy ) );
701 : : // initiate computing norm of residual: (r,r)
702 [ + - ]: 235713 : dot( m_r, m_r,
703 : 235713 : CkCallback( CkReductionTarget(ConjugateGradients,normres), thisProxy ) );
704 : 235713 : }
705 : :
706 : : void
707 : 235713 : ConjugateGradients::normres( tk::real r )
708 : : // *****************************************************************************
709 : : // Compute norm of residual: (r,r)
710 : : //! \param[in] r Dot product, (r,r) (aggregated across all chares)
711 : : // *****************************************************************************
712 : : {
713 : 235713 : m_normr = r;
714 : 235713 : normr_complete();
715 : 235713 : }
716 : :
717 : : void
718 : 235713 : ConjugateGradients::rz( tk::real rz )
719 : : // *****************************************************************************
720 : : // Compute (r,z)
721 : : //! \param[in] rz Dot product, (r,z) (aggregated across all chares)
722 : : // *****************************************************************************
723 : : {
724 : 235713 : m_rho = rz;
725 : :
726 : : // Advance solution: x = x + alpha * p
727 [ + + ]: 54511658 : for (std::size_t i=0; i<m_x.size(); ++i) m_x[i] += m_alpha * m_p[i];
728 : :
729 : : // Communicate solution
730 [ + - ][ + + ]: 471426 : thisProxy[ thisIndex ].wait4x();
731 : :
732 : : // Send solution on chare-boundary nodes to fellow chares
733 [ + + ]: 235713 : if (m_nodeCommMap.empty()) {
734 : 24113 : comx_complete();
735 : : } else {
736 : : auto ncomp = m_A.Ncomp();
737 [ + + ]: 1061348 : for (const auto& [c,n] : m_nodeCommMap) {
738 : 849748 : std::vector< std::vector< tk::real > > xc( n.size() );
739 : : std::size_t j = 0;
740 [ + + ]: 8869844 : for (auto g : n) {
741 [ + - ]: 8020096 : std::vector< tk::real > nx( ncomp );
742 : 8020096 : auto i = tk::cref_find( m_lid, g );
743 [ + + ]: 16045160 : for (std::size_t d=0; d<ncomp; ++d) nx[d] = m_x[ i*ncomp+d ];
744 : 8020096 : xc[j++] = std::move(nx);
745 : : }
746 [ + - ][ + - ]: 1699496 : thisProxy[c].comx( std::vector<std::size_t>(begin(n),end(n)), xc );
[ + - ][ - + ]
[ - - ]
747 : 849748 : }
748 : : }
749 : :
750 : 235713 : ownx_complete();
751 : 235713 : }
752 : :
753 : : void
754 : 849748 : ConjugateGradients::comx( const std::vector< std::size_t >& gid,
755 : : const std::vector< std::vector< tk::real > >& xc )
756 : : // *****************************************************************************
757 : : // Receive contributions to final solution on chare-boundaries
758 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
759 : : //! \param[in] xc Partial contributions at chare-boundary nodes
760 : : // *****************************************************************************
761 : : {
762 : : Assert( xc.size() == gid.size(), "Size mismatch" );
763 : :
764 [ + + ]: 8869844 : for (std::size_t i=0; i<gid.size(); ++i) m_xc[ gid[i] ] += xc[i];
765 : :
766 [ + + ]: 849748 : if (++m_nx == m_nodeCommMap.size()) {
767 : 211600 : m_nx = 0;
768 : 211600 : comx_complete();
769 : : }
770 : 849748 : }
771 : :
772 : : void
773 : 235713 : ConjugateGradients::x()
774 : : // *****************************************************************************
775 : : // Assemble solution on chare boundaries and decide what's next
776 : : // *****************************************************************************
777 : : {
778 : : // Assemble solution on chare boundaries by averaging
779 : : auto ncomp = m_A.Ncomp();
780 [ + + ]: 6808865 : for (const auto& [g,x] : m_xc) {
781 : 6573152 : auto i = tk::cref_find(m_lid,g);
782 [ + + ]: 13151272 : for (std::size_t d=0; d<ncomp; ++d) m_x[i*ncomp+d] += x[d];
783 : 6573152 : auto c = tk::count(m_nodeCommMap,g);
784 [ + + ]: 13151272 : for (std::size_t d=0; d<ncomp; ++d) m_x[i*ncomp+d] /= c;
785 : : }
786 : 235713 : tk::destroy( m_xc );
787 : :
788 : 235713 : ++m_it;
789 [ + + ]: 235713 : auto normb = m_normb > 1.0e-14 ? m_normb : 1.0;
790 : 235713 : auto normr = std::sqrt( m_normr );
791 : :
792 [ + + ][ + + ]: 235713 : if ( m_finished || normr < m_tol*normb || m_it >= m_maxit ) {
[ - + ]
793 : :
794 : 6288 : m_converged = normr > m_tol*normb ? false : true;
795 : :
796 [ + + ]: 6288 : if (m_verbose) {
797 : 51 : std::stringstream c;
798 [ + - ]: 51 : if (m_converged) {
799 [ + - ][ + - ]: 51 : c << " < " << m_tol << ", converged";
800 : : } else {
801 [ - - ]: 0 : c << " > " << m_tol << ", not converged, ||r||/||b|| = "
802 [ - - ][ - - ]: 0 : << normr << '/' << normb;
803 : : }
804 : 51 : tk::Print() << "Xyst> Conjugate gradients it " << m_it
805 [ + - ][ + - ]: 102 : << ", norm = " << normr/normb << c.str() << '\n';
806 : 51 : }
807 : :
808 : : // Restore matrix to state before init()
809 [ + + ]: 6288 : if (m_apply) m_A = m_An;
810 : :
811 : 6288 : m_solved.send( CkDataMsg::buildNew( sizeof(tk::real), &normr ) );
812 : :
813 : : } else {
814 : :
815 [ + + ]: 229425 : if (m_verbose > 1) {
816 : : tk::Print() << "Xyst> Conjugate gradients it "
817 : 114 : << m_it << ", norm = " << normr/normb << '\n';
818 : : }
819 : :
820 : 229425 : next();
821 : :
822 : : }
823 : 235713 : }
824 : :
825 : : #include "NoWarning/conjugategradients.def.h"
|
