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 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Charm++ chare array for distributed conjugate gradients.
9 : : \details Charm++ chare array for asynchronous distributed
10 : : conjugate gradients linear solver.
11 : : \see Y. Saad, Iterative Methods for Sparse Linear Systems: Second Edition,
12 : : ISBN 9780898718003, 2003, Algorithm 6.18, conjugate gradients to solve the
13 : : linear system A * x = b, reproduced here:
14 : :
15 : : Compute r0:=b-A*x0, p0:=r0
16 : : For j=0,1,..., until convergence, do
17 : : alpha_j := (r_j,r_j) / (Ap_j,p_j)
18 : : x_{j+1} := x_j + alpha_j p_j
19 : : r_{j+1} := r_j - alpha_j A p_j
20 : : beta_j := (r_{j+1},r_{j+1}) / (r_j,r_j)
21 : : p_{j+1} := r_{j+1} + beta_j p_j
22 : : end
23 : : */
24 : : // *****************************************************************************
25 : :
26 : : #include <numeric>
27 : : #include <iostream>
28 : :
29 : : #include "Exception.hpp"
30 : : #include "ConjugateGradients.hpp"
31 : : #include "Vector.hpp"
32 : :
33 : : using tk::ConjugateGradients;
34 : :
35 : 27 : ConjugateGradients::ConjugateGradients(
36 : : const CSR& A,
37 : : const std::vector< tk::real >& x,
38 : : const std::vector< tk::real >& b,
39 : : const std::vector< std::size_t >& gid,
40 : : const std::unordered_map< std::size_t, std::size_t >& lid,
41 : 27 : const NodeCommMap& nodecommmap ) :
42 : : m_A( A ),
43 : : m_x( x ),
44 : : m_b( b ),
45 : : m_gid( gid ),
46 : : m_lid( lid ),
47 : : m_nodeCommMap( nodecommmap ),
48 : 0 : m_r( m_A.rsize(), 0.0 ),
49 : : m_rc(),
50 : : m_nr( 0 ),
51 : : m_bc(),
52 : : m_bcc(),
53 : 0 : m_bcmask( m_A.rsize(), 1.0 ),
54 : : m_nb( 0 ),
55 : 0 : m_p( m_A.rsize(), 0.0 ),
56 : 0 : m_q( m_A.rsize(), 0.0 ),
57 : : m_qc(),
58 : : m_nq( 0 ),
59 : : m_initres(),
60 : : m_solved(),
61 : : m_normb( 0.0 ),
62 : : m_it( 0 ),
63 : : m_maxit( 0 ),
64 : : m_rho( 0.0 ),
65 : : m_rho0( 0.0 ),
66 : : m_alpha( 0.0 ),
67 : : m_converged( false ),
68 : : m_xc(),
69 [ + - ][ + - ]: 27 : m_nx( 0 )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
70 : : // *****************************************************************************
71 : : // Constructor
72 : : //! \param[in] A Left hand side matrix of the linear system to solve in Ax=b
73 : : //! \param[in] x Solution (initial guess) of the linear system to solve in Ax=b
74 : : //! \param[in] b Right hand side of the linear system to solve in Ax=b
75 : : //! \param[in] gid Global node ids
76 : : //! \param[in] lid Local node ids associated to global ones
77 : : //! \param[in] nodecommmap Global mesh node IDs shared with other chares
78 : : //! associated to their chare IDs
79 : : // *****************************************************************************
80 : : {
81 : : // Fill in gid and lid for serial solve
82 [ + + ][ + - ]: 27 : if (gid.empty() || lid.empty() || nodecommmap.empty()) {
[ + + ][ + + ]
83 [ + - ]: 3 : m_gid.resize( m_A.rsize()/m_A.Ncomp() );
84 : 3 : std::iota( begin(m_gid), end(m_gid), 0 );
85 [ + + ][ + - ]: 7325 : for (auto g : m_gid) m_lid[g] = g;
86 : : }
87 : :
88 [ - + ][ - - ]: 27 : Assert( m_A.rsize() == m_gid.size()*A.Ncomp(), "Size mismatch" );
[ - - ][ - - ]
89 [ - + ][ - - ]: 27 : Assert( m_x.size() == m_gid.size()*A.Ncomp(), "Size mismatch" );
[ - - ][ - - ]
90 [ - + ][ - - ]: 27 : Assert( m_b.size() == m_gid.size()*A.Ncomp(), "Size mismatch" );
[ - - ][ - - ]
91 : 27 : }
92 : :
93 : : void
94 : 807 : ConjugateGradients::setup( CkCallback c )
95 : : // *****************************************************************************
96 : : // Setup solver
97 : : //! \param[in] c Call to continue with after initialization is complete
98 : : //! \details This function initiates computing the residual (r=b-A*x), its dot
99 : : //! product, and the rhs norm.
100 : : // *****************************************************************************
101 : : {
102 : 807 : m_initres = c;
103 : :
104 : : // initiate computing A * x (for the initial residual)
105 [ + - ]: 807 : thisProxy[ thisIndex ].wait4res();
106 : 807 : residual();
107 : :
108 : : // initiate computing norm of right hand side
109 [ + - ]: 807 : dot( m_b, m_b,
110 [ + - ]: 1614 : CkCallback( CkReductionTarget(ConjugateGradients,normb), thisProxy ) );
111 : 807 : }
112 : :
113 : : void
114 : 10928 : ConjugateGradients::dot( const std::vector< tk::real >& a,
115 : : const std::vector< tk::real >& b,
116 : : CkCallback c )
117 : : // *****************************************************************************
118 : : // Initiate computation of dot product of two vectors
119 : : //! \param[in] a 1st vector of dot product
120 : : //! \param[in] b 2nd vector of dot product
121 : : //! \param[in] c Callback to target with the final result
122 : : // *****************************************************************************
123 : : {
124 [ - + ][ - - ]: 10928 : Assert( a.size() == b.size(), "Size mismatch" );
[ - - ][ - - ]
125 : :
126 : 10928 : tk::real D = 0.0;
127 : 10928 : auto ncomp = m_A.Ncomp();
128 [ + + ]: 17706730 : for (std::size_t i=0; i<a.size()/ncomp; ++i) {
129 : 17695802 : auto incomp = i*ncomp;
130 [ + - ][ + + ]: 17695802 : if (not slave(m_nodeCommMap,m_gid[i],thisIndex))
131 [ + + ]: 64476604 : for (std::size_t d=0; d<ncomp; ++d)
132 : 48316866 : D += a[incomp+d] * b[incomp+d];
133 : : }
134 : :
135 [ + - ]: 10928 : contribute( sizeof(tk::real), &D, CkReduction::sum_double, c );
136 : 10928 : }
137 : :
138 : : void
139 : 807 : ConjugateGradients::normb( tk::real n )
140 : : // *****************************************************************************
141 : : // Compute the norm of the right hand side
142 : : //! \param[in] n Norm of right hand side (aggregated across all chares)
143 : : // *****************************************************************************
144 : : {
145 : 807 : m_normb = std::sqrt(n);
146 : 807 : normb_complete();
147 : 807 : }
148 : :
149 : : void
150 : 807 : ConjugateGradients::residual()
151 : : // *****************************************************************************
152 : : // Initiate A * x for computing the initial residual, r = b - A * x
153 : : // *****************************************************************************
154 : : {
155 : : // Compute own contribution to r = A * x
156 : 807 : m_A.mult( m_x, m_r, m_bcmask );
157 : :
158 : : // Send partial product on chare-boundary nodes to fellow chares
159 [ + + ]: 807 : if (m_nodeCommMap.empty()) {
160 : 63 : comres_complete();
161 : : } else {
162 : 744 : auto ncomp = m_A.Ncomp();
163 [ + + ]: 2416 : for (const auto& [c,n] : m_nodeCommMap) {
164 [ + - ]: 1672 : std::vector< std::vector< tk::real > > rc( n.size() );
165 : 1672 : std::size_t j = 0;
166 [ + + ]: 271548 : for (auto g : n) {
167 [ + - ]: 539752 : std::vector< tk::real > nr( ncomp );
168 [ + - ]: 269876 : auto i = tk::cref_find( m_lid, g );
169 [ + + ]: 1067440 : for (std::size_t d=0; d<ncomp; ++d) nr[d] = m_r[ i*ncomp+d ];
170 : 269876 : rc[j++] = std::move(nr);
171 : : }
172 [ + - ][ + - ]: 1672 : thisProxy[c].comres( std::vector<std::size_t>(begin(n),end(n)), rc );
[ + - ]
173 : : }
174 : : }
175 : :
176 : 807 : ownres_complete();
177 : 807 : }
178 : :
179 : : void
180 : 1672 : ConjugateGradients::comres( const std::vector< std::size_t >& gid,
181 : : const std::vector< std::vector< tk::real > >& rc )
182 : : // *****************************************************************************
183 : : // Receive contributions to A * x on chare-boundaries
184 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
185 : : //! \param[in] rc Partial contributions at chare-boundary nodes
186 : : // *****************************************************************************
187 : : {
188 [ - + ][ - - ]: 1672 : Assert( rc.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
189 : :
190 : : using tk::operator+=;
191 : :
192 [ + + ]: 271548 : for (std::size_t i=0; i<gid.size(); ++i)
193 : 269876 : m_rc[ gid[i] ] += rc[i];
194 : :
195 [ + + ]: 1672 : if (++m_nr == m_nodeCommMap.size()) {
196 : 744 : m_nr = 0;
197 : 744 : comres_complete();
198 : : }
199 : 1672 : }
200 : :
201 : : void
202 : 807 : ConjugateGradients::initres()
203 : : // *****************************************************************************
204 : : // Finish computing the initial residual, r = b - A * x
205 : : // *****************************************************************************
206 : : {
207 : : // Combine own and communicated contributions to r = A * x
208 : 807 : auto ncomp = m_A.Ncomp();
209 [ + + ]: 240861 : for (const auto& [gid,r] : m_rc) {
210 [ + - ]: 240054 : auto i = tk::cref_find( m_lid, gid );
211 [ + + ]: 951004 : for (std::size_t c=0; c<ncomp; ++c) m_r[i*ncomp+c] += r[c];
212 : : }
213 : 807 : tk::destroy( m_rc );
214 : :
215 : : // Finish computing the initial residual, r = b - A * x
216 [ + + ]: 3763550 : for (auto& r : m_r) r *= -1.0;
217 : 807 : m_r += m_b;
218 : :
219 : : // Initialize p
220 : 807 : m_p = m_r;
221 : :
222 : : // initiate computing the dot product of the initial residual, rho = (r,r)
223 [ + - ]: 807 : dot( m_r, m_r,
224 [ + - ]: 1614 : CkCallback( CkReductionTarget(ConjugateGradients,rho), thisProxy ) );
225 : 807 : }
226 : :
227 : : void
228 : 807 : ConjugateGradients::rho( tk::real r )
229 : : // *****************************************************************************
230 : : // Compute rho = (r,r)
231 : : //! \param[in] r Dot product, rho = (r,r) (aggregated across all chares)
232 : : // *****************************************************************************
233 : : {
234 : : // store dot product of residual
235 : 807 : m_rho = r;
236 : :
237 : : // send back rhs norm to caller
238 : 807 : m_initres.send( CkDataMsg::buildNew( sizeof(tk::real), &m_normb ) );
239 : 807 : }
240 : :
241 : : void
242 : 801 : ConjugateGradients::init(
243 : : const std::vector< tk::real >& x,
244 : : const std::vector< tk::real >& b,
245 : : const std::unordered_map< std::size_t,
246 : : std::vector< std::pair< bool, tk::real > > >& bc,
247 : : std::size_t ignorebc,
248 : : CkCallback cb )
249 : : // *****************************************************************************
250 : : // Initialize linear solve: set initial guess and boundary conditions
251 : : //! \param[in] x Initial guess
252 : : //! \param[in] b Right hand side vector
253 : : //! \param[in] bc Local node ids and associated Dirichlet BCs
254 : : //! \param[in] ignorebc True if applyin BCs should be skipped
255 : : //! \param[in] cb Call to continue with when initialized and ready for a solve
256 : : //! \details This function allows setting the initial guess and boundary
257 : : //! conditions, followed by computing the initial residual and the rhs norm.
258 : : // *****************************************************************************
259 : : {
260 : : // Optionally set initial guess
261 [ + + ]: 801 : if (not x.empty()) m_x = x;
262 : :
263 : : // Optionally update rhs
264 [ + + ]: 801 : if (not b.empty()) m_b = b;
265 : :
266 [ + + ]: 801 : if (ignorebc) {
267 : :
268 [ + - ]: 108 : setup( cb );
269 : :
270 : : } else {
271 : :
272 : : // Store incoming BCs
273 : 693 : m_bc = bc;
274 : :
275 : : // Get ready to communicate boundary conditions. This is necessary because
276 : : // there can be nodes a chare contributes to but does not apply BCs on. This
277 : : // happens if a node is in the node communication map but not on the list of
278 : : // incoming BCs on this chare. To have all chares share the same view on all
279 : : // BC nodes, we send the global node ids together with the Dirichlet BCs at
280 : : // which BCs are set to those fellow chares that also contribute to those BC
281 : : // nodes. Only after this communication step we apply the BCs on the matrix,
282 : : // which then will correctly setup the BC rows that exist on multiple chares
283 : : // (which now will be the same as the results of making the BCs consistent
284 : : // across all chares that contribute.
285 [ + - ]: 693 : thisProxy[ thisIndex ].wait4bc();
286 : :
287 : : // Send boundary conditions to those who contribute to those rows
288 [ + + ]: 693 : if (m_nodeCommMap.empty()) {
289 : 61 : combc_complete();
290 : : } else {
291 [ + + ]: 2138 : for (const auto& [c,n] : m_nodeCommMap) {
292 : : std::unordered_map< std::size_t,
293 : 1506 : std::vector< std::pair< bool, tk::real > > > expbc;
294 [ + + ]: 263462 : for (auto g : n) {
295 [ + - ]: 261956 : auto i = tk::cref_find( m_lid, g );
296 [ + - ]: 261956 : auto j = bc.find(i);
297 [ + + ][ + - ]: 261956 : if (j != end(bc)) expbc[g] = j->second;
[ + - ]
298 : : }
299 [ + - ][ + - ]: 1506 : thisProxy[c].combc( expbc );
300 : : }
301 : : }
302 : :
303 [ + - ]: 693 : ownbc_complete( cb );
304 : :
305 : : }
306 : 801 : }
307 : :
308 : : void
309 : 1506 : ConjugateGradients::combc(
310 : : const std::unordered_map< std::size_t,
311 : : std::vector< std::pair< bool, tk::real > > >& bc )
312 : : // *****************************************************************************
313 : : // Receive contributions to boundary conditions on chare-boundaries
314 : : //! \param[in] bc Contributions to boundary conditions
315 : : // *****************************************************************************
316 : : {
317 [ + + ][ + - ]: 11400 : for (const auto& [g,dirbc] : bc) m_bcc[ tk::cref_find(m_lid,g) ] = dirbc;
[ + - ][ + - ]
318 : :
319 [ + + ]: 1506 : if (++m_nb == m_nodeCommMap.size()) {
320 : 632 : m_nb = 0;
321 : 632 : combc_complete();
322 : : }
323 : 1506 : }
324 : :
325 : : void
326 : 693 : ConjugateGradients::apply( CkCallback cb )
327 : : // *****************************************************************************
328 : : // Apply boundary conditions
329 : : //! \param[in] cb Call to continue with after applying the BCs is complete
330 : : // *****************************************************************************
331 : : {
332 : : // Merge own and received contributions to boundary conditions
333 [ + + ][ + - ]: 9564 : for (const auto& [i,dirbc] : m_bcc) m_bc[i] = dirbc;
[ + - ]
334 : 693 : tk::destroy( m_bcc );
335 : :
336 : 693 : auto ncomp = m_A.Ncomp();
337 : :
338 : : // Setup Dirichlet BC map as contiguous mask
339 [ + + ]: 37157 : for (const auto& [i,bc] : m_bc)
340 [ + + ]: 129984 : for (std::size_t j=0; j<ncomp; ++j)
341 : 93520 : m_bcmask[i*ncomp+j] = 0.0;
342 : :
343 : : // Apply Dirichlet BCs on matrix and rhs
344 [ + + ]: 37157 : for (const auto& [i,dirbc] : m_bc) {
345 [ + + ]: 129984 : for (std::size_t j=0; j<ncomp; ++j) {
346 [ + + ]: 93520 : if (dirbc[j].first) {
347 [ + - ]: 93124 : m_A.dirichlet( i, m_gid, m_nodeCommMap, j );
348 : 93124 : m_b[i*ncomp+j] = dirbc[j].second;
349 : : }
350 : : }
351 : : }
352 : :
353 : : // Recompute initial residual (r=b-A*x), its dot product, and the rhs norm
354 [ + - ]: 693 : setup( cb );
355 : 693 : }
356 : :
357 : : void
358 : 807 : ConjugateGradients::solve( std::size_t maxit, tk::real tol, CkCallback c )
359 : : // *****************************************************************************
360 : : // Solve linear system
361 : : //! \param[in] maxit Max iteration count
362 : : //! \param[in] tol Stop tolerance
363 : : //! \param[in] c Call to continue with after solve is complete
364 : : // *****************************************************************************
365 : : {
366 : 807 : m_maxit = maxit;
367 : 807 : m_tol = tol;
368 : 807 : m_solved = c;
369 : 807 : m_it = 0;
370 : :
371 : 807 : next();
372 : 807 : }
373 : :
374 : : void
375 : 4657 : ConjugateGradients::next()
376 : : // *****************************************************************************
377 : : // Start next linear solver iteration
378 : : // *****************************************************************************
379 : : {
380 [ + + ]: 4657 : if (m_it == 0) m_alpha = 0.0; else m_alpha = m_rho/m_rho0;
381 : 4657 : m_rho0 = m_rho;
382 : :
383 : : // compute p = r + alpha * p
384 [ + + ]: 22677277 : for (std::size_t i=0; i<m_p.size(); ++i) m_p[i] = m_r[i] + m_alpha * m_p[i];
385 : :
386 : : // initiate computing q = A * p
387 [ + - ]: 4657 : thisProxy[ thisIndex ].wait4q();
388 : 4657 : qAp();
389 : 4657 : }
390 : :
391 : :
392 : : void
393 : 4657 : ConjugateGradients::qAp()
394 : : // *****************************************************************************
395 : : // Initiate computing q = A * p
396 : : // *****************************************************************************
397 : : {
398 : : // Compute own contribution to q = A * p
399 : 4657 : m_A.mult( m_p, m_q, m_bcmask );
400 : :
401 : : // Send partial product on chare-boundary nodes to fellow chares
402 [ + + ]: 4657 : if (m_nodeCommMap.empty()) {
403 : 385 : comq_complete();
404 : : } else {
405 : 4272 : auto ncomp = m_A.Ncomp();
406 [ + + ]: 13356 : for (const auto& [c,n] : m_nodeCommMap) {
407 [ + - ]: 9084 : std::vector< std::vector< tk::real > > qc( n.size() );
408 : 9084 : std::size_t j = 0;
409 [ + + ]: 1404522 : for (auto g : n) {
410 [ + - ]: 2790876 : std::vector< tk::real > nq( ncomp );
411 [ + - ]: 1395438 : auto i = tk::cref_find( m_lid, g );
412 [ + + ]: 5527992 : for (std::size_t d=0; d<ncomp; ++d) nq[d] = m_q[ i*ncomp+d ];
413 : 1395438 : qc[j++] = std::move(nq);
414 : : }
415 [ + - ][ + - ]: 9084 : thisProxy[c].comq( std::vector<std::size_t>(begin(n),end(n)), qc );
[ + - ]
416 : : }
417 : : }
418 : :
419 : 4657 : ownq_complete();
420 : 4657 : }
421 : :
422 : : void
423 : 9084 : ConjugateGradients::comq( const std::vector< std::size_t >& gid,
424 : : const std::vector< std::vector< tk::real > >& qc )
425 : : // *****************************************************************************
426 : : // Receive contributions to q = A * p on chare-boundaries
427 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
428 : : //! \param[in] qc Partial contributions at chare-boundary nodes
429 : : // *****************************************************************************
430 : : {
431 [ - + ][ - - ]: 9084 : Assert( qc.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
432 : :
433 : : using tk::operator+=;
434 : :
435 [ + + ]: 1404522 : for (std::size_t i=0; i<gid.size(); ++i)
436 : 1395438 : m_qc[ gid[i] ] += qc[i];
437 : :
438 [ + + ]: 9084 : if (++m_nq == m_nodeCommMap.size()) {
439 : 4272 : m_nq = 0;
440 : 4272 : comq_complete();
441 : : }
442 : 9084 : }
443 : :
444 : : void
445 : 4657 : ConjugateGradients::q()
446 : : // *****************************************************************************
447 : : // Finish computing q = A * p
448 : : // *****************************************************************************
449 : : {
450 : : // Combine own and communicated contributions to q = A * p
451 : 4657 : auto ncomp = m_A.Ncomp();
452 [ + + ]: 1251767 : for (const auto& [gid,q] : m_qc) {
453 [ + - ]: 1247110 : auto i = tk::cref_find( m_lid, gid );
454 [ + + ]: 4947376 : for (std::size_t c=0; c<ncomp; ++c)
455 : 3700266 : m_q[i*ncomp+c] += q[c];
456 : : }
457 : 4657 : tk::destroy( m_qc );
458 : :
459 : : // initiate computing (p,q)
460 [ + - ]: 4657 : dot( m_p, m_q,
461 [ + - ]: 9314 : CkCallback( CkReductionTarget(ConjugateGradients,pq), thisProxy ) );
462 : 4657 : }
463 : :
464 : : void
465 : 4657 : ConjugateGradients::pq( tk::real d )
466 : : // *****************************************************************************
467 : : // Compute the dot product (p,q)
468 : : //! \param[in] d Dot product of (p,q) (aggregated across all chares)
469 : : // *****************************************************************************
470 : : {
471 : : // If (p,q)=0, then p and q are orthogonal and the system either has a trivial
472 : : // solution, x=x0, or the BCs are incomplete or wrong, in either case the
473 : : // solve cannot continue.
474 : 4657 : const auto eps = std::numeric_limits< tk::real >::epsilon();
475 [ + + ]: 4657 : if (std::abs(d) < eps) {
476 : 9 : m_it = m_maxit;
477 : 9 : m_alpha = 0.0;
478 : : } else {
479 : 4648 : m_alpha = m_rho / d;
480 : : }
481 : :
482 : : // compute r = r - alpha * q
483 [ + + ]: 22677277 : for (std::size_t i=0; i<m_r.size(); ++i) m_r[i] -= m_alpha * m_q[i];
484 : :
485 : : // initiate computing norm of residual: (r,r)
486 [ + - ]: 4657 : dot( m_r, m_r,
487 [ + - ]: 9314 : CkCallback( CkReductionTarget(ConjugateGradients,normres), thisProxy ) );
488 : 4657 : }
489 : :
490 : : void
491 : 4657 : ConjugateGradients::normres( tk::real r )
492 : : // *****************************************************************************
493 : : // Compute norm of residual: (r,r)
494 : : //! \param[in] r Dot product, (r,r) (aggregated across all chares)
495 : : // *****************************************************************************
496 : : {
497 : 4657 : m_rho = r;
498 : :
499 : : // Advance solution: x = x + alpha * p
500 [ + + ]: 22677277 : for (std::size_t i=0; i<m_x.size(); ++i) m_x[i] += m_alpha * m_p[i];
501 : :
502 : : // Communicate solution
503 [ + - ]: 4657 : thisProxy[ thisIndex ].wait4x();
504 : :
505 : : // Send solution on chare-boundary nodes to fellow chares
506 [ + + ]: 4657 : if (m_nodeCommMap.empty()) {
507 : 385 : comx_complete();
508 : : } else {
509 : 4272 : auto ncomp = m_A.Ncomp();
510 [ + + ]: 13356 : for (const auto& [c,n] : m_nodeCommMap) {
511 [ + - ]: 9084 : std::vector< std::vector< tk::real > > xc( n.size() );
512 : 9084 : std::size_t j = 0;
513 [ + + ]: 1404522 : for (auto g : n) {
514 [ + - ]: 2790876 : std::vector< tk::real > nx( ncomp );
515 [ + - ]: 1395438 : auto i = tk::cref_find( m_lid, g );
516 [ + + ]: 5527992 : for (std::size_t d=0; d<ncomp; ++d) nx[d] = m_x[ i*ncomp+d ];
517 : 1395438 : xc[j++] = std::move(nx);
518 : : }
519 [ + - ][ + - ]: 9084 : thisProxy[c].comx( std::vector<std::size_t>(begin(n),end(n)), xc );
[ + - ]
520 : : }
521 : : }
522 : :
523 : 4657 : ownx_complete();
524 : 4657 : }
525 : :
526 : : void
527 : 9084 : ConjugateGradients::comx( const std::vector< std::size_t >& gid,
528 : : const std::vector< std::vector< tk::real > >& xc )
529 : : // *****************************************************************************
530 : : // Receive contributions to final solution on chare-boundaries
531 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
532 : : //! \param[in] xc Partial contributions at chare-boundary nodes
533 : : // *****************************************************************************
534 : : {
535 [ - + ][ - - ]: 9084 : Assert( xc.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
536 : :
537 [ + + ]: 1404522 : for (std::size_t i=0; i<gid.size(); ++i) m_xc[ gid[i] ] += xc[i];
538 : :
539 [ + + ]: 9084 : if (++m_nx == m_nodeCommMap.size()) {
540 : 4272 : m_nx = 0;
541 : 4272 : comx_complete();
542 : : }
543 : 9084 : }
544 : :
545 : : void
546 : 4657 : ConjugateGradients::x()
547 : : // *****************************************************************************
548 : : // Assemble solution on chare boundaries
549 : : // *****************************************************************************
550 : : {
551 : : // Assemble solution on chare boundaries by averaging
552 : 4657 : auto ncomp = m_A.Ncomp();
553 [ + + ]: 1251767 : for (const auto& [g,x] : m_xc) {
554 [ + - ]: 1247110 : auto i = tk::cref_find(m_lid,g);
555 [ + + ]: 4947376 : for (std::size_t d=0; d<ncomp; ++d) m_x[i*ncomp+d] += x[d];
556 [ + - ]: 1247110 : auto c = tk::count(m_nodeCommMap,g);
557 [ + + ]: 4947376 : for (std::size_t d=0; d<ncomp; ++d) m_x[i*ncomp+d] /= c;
558 : : }
559 : 4657 : tk::destroy( m_xc );
560 : :
561 : 4657 : ++m_it;
562 [ + + ]: 4657 : auto normb = m_normb > 1.0e-14 ? m_normb : 1.0;
563 : 4657 : auto normr = std::sqrt( m_rho );
564 : :
565 [ + + ][ + + ]: 4657 : if ( m_it < m_maxit && normr > m_tol*normb ) {
566 : :
567 [ + - ]: 3850 : next();
568 : :
569 : : } else {
570 : :
571 [ + + ][ - + ]: 807 : m_converged = m_it == m_maxit && normr > m_tol*normb ? false : true;
572 [ + - ][ + - ]: 807 : m_solved.send( CkDataMsg::buildNew( sizeof(tk::real), &normr ) );
573 : :
574 : : }
575 : 4657 : }
576 : :
577 : : #include "NoWarning/conjugategradients.def.h"
|
