Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/LinearSolver/BiCG.cpp
4 : : created 2025, Christopher Long
5 : : \copyright 2012-2015 J. Bakosi,
6 : : 2016-2018 Los Alamos National Security, LLC.,
7 : : 2019-2021 Triad National Security, LLC.
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.
12 : : \see Y. Saad, Iterative Methods for Sparse Linear Systems: Second Edition,
13 : : ISBN 9780898718003, 2003, Algorithm 6.18, conjugate gradients to solve the
14 : : linear system A * x = b, reproduced here:
15 : :
16 : : Compute r0:=b-A*x0, p0:=r0
17 : : Pick rhat0 = r0
18 : : Let rho0 = rhat0^T \dot r0
19 : : For j=0,1,..., until convergence, do
20 : : alpha := rho_j / (rhat0, Ap_j)
21 : : h := x_j + alpha p_j
22 : : s := r_j - alpha A p_j
23 : : Exit if s is sufficient
24 : : t = As
25 : : w = (t,s)/(t,t) //replaces beta
26 : :
27 : : x_{j+1} := h + ws
28 : : r_{j+1} := s - wt
29 : : Exit if r_{j+1} is sufficient
30 : : rho_{j+1} := (rhat0, r_{j+1})
31 : : beta := (rho_{j+1}/rho_{j})*(alpha/w)
32 : : p_{j+1} := r_{j+1} + beta(p_j - wAp_j)
33 : : end
34 : :
35 : : Usage notes and interfaces into BiCG:
36 : : (1) BiCG::converged() is an accessor used to query for convergence
37 : : (2) BiCG::init() is a charm entry method that initializes the solver, and
38 : : returns execution to client code using CkCallback input arg
39 : : (3) BiCG::solve() actually solves the linear system, and returns execution
40 : : to client code using CkCallback input arg
41 : : (4) BiCG::solution() is an accessor to the solution of linear solver
42 : : CAUTION! solution() is not a ref, so it creates a copy of the entire
43 : : solution vector
44 : : */
45 : : // *****************************************************************************
46 : :
47 : : #include <numeric>
48 : : #include <iostream>
49 : :
50 : : #include "Exception.hpp"
51 : : #include "BiCG.hpp"
52 : : #include "Vector.hpp"
53 : :
54 : : using tk::BiCG;
55 : :
56 : 0 : BiCG::BiCG(
57 : : const CSR& A,
58 : : const std::vector< tk::real >& x,
59 : : const std::vector< tk::real >& b,
60 : : const std::vector< std::size_t >& gid,
61 : : const std::unordered_map< std::size_t, std::size_t >& lid,
62 : 0 : const NodeCommMap& nodecommmap ) :
63 [ - - ]: 0 : m_A( A ),
64 [ - - ]: 0 : m_x( x ),
65 [ - - ]: 0 : m_b( b ),
66 [ - - ]: 0 : m_gid( gid ),
67 : : m_lid( lid ),
68 : : m_nodeCommMap( nodecommmap ),
69 [ - - ]: 0 : m_r( m_A.rsize(), 0.0 ),
70 [ - - ][ - - ]: 0 : m_r0( m_A.rsize(), 0.0 ),
71 [ - - ]: 0 : m_rc(),
72 : 0 : m_nr( 0 ),
73 : 0 : m_bc(),
74 : 0 : m_bcc(),
75 [ - - ]: 0 : m_bcmask( m_A.rsize(), 1.0 ),
76 : 0 : m_nb( 0 ),
77 [ - - ][ - - ]: 0 : m_p( m_A.rsize(), 0.0 ),
78 [ - - ][ - - ]: 0 : m_q( m_A.rsize(), 0.0 ),
79 [ - - ][ - - ]: 0 : m_t( m_A.rsize(), 0.0 ),
80 [ - - ]: 0 : m_qc(),
81 : 0 : m_tc(), //can I reuse this?
82 : 0 : m_nt( 0 ),
83 : 0 : m_nq( 0 ),
84 [ - - ]: 0 : m_initres(),
85 : 0 : m_solved(),
86 : 0 : m_normb( 0.0 ),
87 : 0 : m_it( 0 ),
88 : 0 : m_maxit( 0 ),
89 : 0 : m_rho( 0.0 ),
90 : 0 : m_rho0( 0.0 ),
91 : 0 : m_alpha( 0.0 ),
92 : 0 : m_omega( 0.0 ),
93 : 0 : m_converged( false ),
94 : 0 : m_xc(),
95 : 0 : m_x2c(),
96 : 0 : m_nx( 0 ),
97 [ - - ][ - - ]: 0 : m_nx2( 0 )
98 : : // *****************************************************************************
99 : : // Constructor
100 : : //! \param[in] A Left hand side matrix of the linear system to solve in Ax=b
101 : : //! \param[in] x Solution (initial guess) of the linear system to solve in Ax=b
102 : : //! \param[in] b Right hand side of the linear system to solve in Ax=b
103 : : //! \param[in] gid Global node ids
104 : : //! \param[in] lid Local node ids associated to global ones
105 : : //! \param[in] nodecommmap Global mesh node IDs shared with other chares
106 : : //! associated to their chare IDs
107 : : // *****************************************************************************
108 : : {
109 : : // Fill in gid and lid for serial solve
110 [ - - ][ - - ]: 0 : if (gid.empty() || lid.empty() || nodecommmap.empty()) {
[ - - ]
111 [ - - ]: 0 : m_gid.resize( m_A.rsize()/m_A.Ncomp() );
112 : : std::iota( begin(m_gid), end(m_gid), 0 );
113 [ - - ][ - - ]: 0 : for (auto g : m_gid) m_lid[g] = g;
114 : : }
115 : :
116 : : Assert( m_A.rsize() == m_gid.size()*A.Ncomp(), "Size mismatch" );
117 : : Assert( m_x.size() == m_gid.size()*A.Ncomp(), "Size mismatch" );
118 : : Assert( m_b.size() == m_gid.size()*A.Ncomp(), "Size mismatch" );
119 : 0 : }
120 : :
121 : : void
122 : 0 : BiCG::setup( CkCallback c )
123 : : // *****************************************************************************
124 : : // Setup solver
125 : : //! \param[in] c Call to continue with after initialization is complete
126 : : //! \details This function initiates computing the residual (r=b-A*x), its dot
127 : : //! product, and the rhs norm.
128 : : // *****************************************************************************
129 : : {
130 : 0 : m_initres = c;
131 : :
132 : : // initiate computing A * x (for the initial residual)
133 [ - - ]: 0 : thisProxy[ thisIndex ].wait4res();
134 : 0 : residual();
135 : :
136 : : // initiate computing norm of right hand side
137 [ - - ]: 0 : dot( m_b, m_b,
138 : 0 : CkCallback( CkReductionTarget(BiCG,normb), thisProxy ) );
139 : 0 : }
140 : :
141 : : void
142 : 0 : BiCG::dot( const std::vector< tk::real >& a,
143 : : const std::vector< tk::real >& b,
144 : : CkCallback c )
145 : : // *****************************************************************************
146 : : // Initiate computation of dot product of two vectors
147 : : //! \param[in] a 1st vector of dot product
148 : : //! \param[in] b 2nd vector of dot product
149 : : //! \param[in] c Callback to target with the final result
150 : : // *****************************************************************************
151 : : {
152 : : Assert( a.size() == b.size(), "Size mismatch" );
153 : :
154 : 0 : tk::real D = 0.0;
155 : : auto ncomp = m_A.Ncomp();
156 [ - - ]: 0 : for (std::size_t i=0; i<a.size()/ncomp; ++i) {
157 : 0 : auto incomp = i*ncomp;
158 [ - - ]: 0 : if (not slave(m_nodeCommMap,m_gid[i],thisIndex))
159 [ - - ]: 0 : for (std::size_t d=0; d<ncomp; ++d)
160 : 0 : D += a[incomp+d] * b[incomp+d];
161 : : }
162 : :
163 : 0 : contribute( sizeof(tk::real), &D, CkReduction::sum_double, c );
164 : 0 : }
165 : :
166 : : void
167 : 0 : BiCG::normb( tk::real n )
168 : : // *****************************************************************************
169 : : // Compute the norm of the right hand side
170 : : //! \param[in] n Norm of right hand side (aggregated across all chares)
171 : : // *****************************************************************************
172 : : {
173 : 0 : m_normb = std::sqrt(n);
174 : 0 : normb_complete();
175 : 0 : }
176 : :
177 : : void
178 : 0 : BiCG::residual()
179 : : // *****************************************************************************
180 : : // Initiate A * x for computing the initial residual, r = b - A * x
181 : : // *****************************************************************************
182 : : {
183 : : // Compute own contribution to r = A * x
184 : 0 : m_A.mult( m_x, m_r, m_bcmask );
185 : :
186 : : // Send partial product on chare-boundary nodes to fellow chares
187 [ - - ]: 0 : if (m_nodeCommMap.empty()) {
188 : 0 : comres_complete();
189 : : } else {
190 : : auto ncomp = m_A.Ncomp();
191 [ - - ]: 0 : for (const auto& [c,n] : m_nodeCommMap) {
192 : 0 : std::vector< std::vector< tk::real > > rc( n.size() );
193 : : std::size_t j = 0;
194 [ - - ]: 0 : for (auto g : n) {
195 [ - - ]: 0 : std::vector< tk::real > nr( ncomp );
196 : 0 : auto i = tk::cref_find( m_lid, g );
197 [ - - ]: 0 : for (std::size_t d=0; d<ncomp; ++d) nr[d] = m_r[ i*ncomp+d ];
198 : 0 : rc[j++] = std::move(nr);
199 : : }
200 [ - - ][ - - ]: 0 : thisProxy[c].comres( std::vector<std::size_t>(begin(n),end(n)), rc );
[ - - ][ - - ]
[ - - ]
201 : 0 : }
202 : : }
203 : :
204 : 0 : ownres_complete();
205 : 0 : }
206 : :
207 : : void
208 : 0 : BiCG::comres( const std::vector< std::size_t >& gid,
209 : : const std::vector< std::vector< tk::real > >& rc )
210 : : // *****************************************************************************
211 : : // Receive contributions to A * x on chare-boundaries
212 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
213 : : //! \param[in] rc Partial contributions at chare-boundary nodes
214 : : // *****************************************************************************
215 : : {
216 : : Assert( rc.size() == gid.size(), "Size mismatch" );
217 : :
218 : : using tk::operator+=;
219 : :
220 [ - - ]: 0 : for (std::size_t i=0; i<gid.size(); ++i)
221 : 0 : m_rc[ gid[i] ] += rc[i];
222 : :
223 [ - - ]: 0 : if (++m_nr == m_nodeCommMap.size()) {
224 : 0 : m_nr = 0;
225 : 0 : comres_complete();
226 : : }
227 : 0 : }
228 : :
229 : : void
230 : 0 : BiCG::initres()
231 : : // *****************************************************************************
232 : : // Finish computing the initial residual, r = b - A * x
233 : : // *****************************************************************************
234 : : {
235 : : // Combine own and communicated contributions to r = A * x
236 : : auto ncomp = m_A.Ncomp();
237 [ - - ]: 0 : for (const auto& [gid,r] : m_rc) {
238 : 0 : auto i = tk::cref_find( m_lid, gid );
239 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) m_r[i*ncomp+c] += r[c];
240 : : }
241 : 0 : tk::destroy( m_rc );
242 : :
243 : : // Finish computing the initial residual, r = b - A * x
244 [ - - ]: 0 : for (auto& r : m_r) r *= -1.0;
245 : 0 : m_r += m_b;
246 : :
247 : : // Initialize p
248 : 0 : m_p = m_r;
249 : : //Initialize r0
250 : 0 : m_r0 = m_r;
251 : : // initiate computing the dot product of the initial residual, rho = (r,r)
252 [ - - ]: 0 : dot( m_r0, m_r,
253 : 0 : CkCallback( CkReductionTarget(BiCG,rho), thisProxy ) );
254 : 0 : }
255 : :
256 : : void
257 : 0 : BiCG::rho( tk::real r )
258 : : // *****************************************************************************
259 : : // Compute rho = (r,r)
260 : : //! \param[in] r Dot product, rho = (r,r) (aggregated across all chares)
261 : : // *****************************************************************************
262 : : {
263 : : // store dot product of residual
264 : 0 : m_rho = r;
265 : :
266 : : // send back rhs norm to caller
267 : 0 : m_initres.send( CkDataMsg::buildNew( sizeof(tk::real), &m_normb ) );
268 : 0 : }
269 : :
270 : : void
271 [ - - ]: 0 : BiCG::init(
272 : : const std::vector< tk::real >& x,
273 : : const std::vector< tk::real >& b,
274 : : const std::unordered_map< std::size_t,
275 : : std::vector< std::pair< bool, tk::real > > >& bc,
276 : : std::size_t ignorebc,
277 : : CkCallback cb )
278 : : // *****************************************************************************
279 : : // Initialize linear solve: set initial guess and boundary conditions
280 : : //! \param[in] x Initial guess
281 : : //! \param[in] b Right hand side vector
282 : : //! \param[in] bc Local node ids and associated Dirichlet BCs
283 : : //! \param[in] ignorebc True if applyin BCs should be skipped
284 : : //! \param[in] cb Call to continue with when initialized and ready for a solve
285 : : //! \details This function allows setting the initial guess and boundary
286 : : //! conditions, followed by computing the initial residual and the rhs norm.
287 : : // *****************************************************************************
288 : : {
289 : : // Optionally set initial guess
290 [ - - ]: 0 : if (not x.empty()) m_x = x;
291 : :
292 : : // Optionally update rhs
293 [ - - ]: 0 : if (not b.empty()) m_b = b;
294 : :
295 [ - - ]: 0 : if (ignorebc) {
296 : :
297 [ - - ]: 0 : setup( cb );
298 : :
299 : : } else {
300 : :
301 : : // Store incoming BCs
302 : : m_bc = bc;
303 : :
304 : : // Get ready to communicate boundary conditions. This is necessary because
305 : : // there can be nodes a chare contributes to but does not apply BCs on. This
306 : : // happens if a node is in the node communication map but not on the list of
307 : : // incoming BCs on this chare. To have all chares share the same view on all
308 : : // BC nodes, we send the global node ids together with the Dirichlet BCs at
309 : : // which BCs are set to those fellow chares that also contribute to those BC
310 : : // nodes. Only after this communication step we apply the BCs on the matrix,
311 : : // which then will correctly setup the BC rows that exist on multiple chares
312 : : // (which now will be the same as the results of making the BCs consistent
313 : : // across all chares that contribute.
314 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4bc();
315 : :
316 : : // Send boundary conditions to those who contribute to those rows
317 [ - - ]: 0 : if (m_nodeCommMap.empty()) {
318 : 0 : combc_complete();
319 : : } else {
320 [ - - ]: 0 : for (const auto& [c,n] : m_nodeCommMap) {
321 : : std::unordered_map< std::size_t,
322 : : std::vector< std::pair< bool, tk::real > > > expbc;
323 [ - - ][ - - ]: 0 : for (auto g : n) {
324 : 0 : auto i = tk::cref_find( m_lid, g );
325 : : auto j = bc.find(i);
326 [ - - ][ - - ]: 0 : if (j != end(bc)) expbc[g] = j->second;
[ - - ]
327 : : }
328 [ - - ][ - - ]: 0 : thisProxy[c].combc( expbc );
329 : : }
330 : : }
331 : :
332 [ - - ]: 0 : ownbc_complete( cb );
333 : :
334 : : }
335 : 0 : }
336 : :
337 : : void
338 : 0 : BiCG::combc(
339 : : const std::unordered_map< std::size_t,
340 : : std::vector< std::pair< bool, tk::real > > >& bc )
341 : : // *****************************************************************************
342 : : // Receive contributions to boundary conditions on chare-boundaries
343 : : //! \param[in] bc Contributions to boundary conditions
344 : : // *****************************************************************************
345 : : {
346 [ - - ]: 0 : for (const auto& [g,dirbc] : bc) m_bcc[ tk::cref_find(m_lid,g) ] = dirbc;
347 : :
348 [ - - ]: 0 : if (++m_nb == m_nodeCommMap.size()) {
349 : 0 : m_nb = 0;
350 : 0 : combc_complete();
351 : : }
352 : 0 : }
353 : :
354 : : void
355 : 0 : BiCG::apply( CkCallback cb )
356 : : // *****************************************************************************
357 : : // Apply boundary conditions
358 : : //! \param[in] cb Call to continue with after applying the BCs is complete
359 : : // *****************************************************************************
360 : : {
361 : : // Merge own and received contributions to boundary conditions
362 [ - - ]: 0 : for (const auto& [i,dirbc] : m_bcc) m_bc[i] = dirbc;
363 : 0 : tk::destroy( m_bcc );
364 : :
365 : : auto ncomp = m_A.Ncomp();
366 : :
367 : : // Setup Dirichlet BC map as contiguous mask
368 [ - - ]: 0 : for (const auto& [i,bc] : m_bc)
369 [ - - ]: 0 : for (std::size_t j=0; j<ncomp; ++j)
370 : 0 : m_bcmask[i*ncomp+j] = 0.0;
371 : :
372 : : // Apply Dirichlet BCs on matrix and rhs
373 [ - - ]: 0 : for (const auto& [i,dirbc] : m_bc) {
374 [ - - ]: 0 : for (std::size_t j=0; j<ncomp; ++j) {
375 [ - - ]: 0 : if (dirbc[j].first) {
376 : 0 : m_A.dirichlet( i, m_gid, m_nodeCommMap, j );
377 : 0 : m_b[i*ncomp+j] = dirbc[j].second;
378 : : }
379 : : }
380 : : }
381 : :
382 : : // Recompute initial residual (r=b-A*x), its dot product, and the rhs norm
383 [ - - ]: 0 : setup( cb );
384 : 0 : }
385 : :
386 : : void
387 : 0 : BiCG::solve( std::size_t maxit, tk::real tol, CkCallback c )
388 : : // *****************************************************************************
389 : : // Solve linear system
390 : : //! \param[in] maxit Max iteration count
391 : : //! \param[in] tol Stop tolerance
392 : : //! \param[in] c Call to continue with after solve is complete
393 : : // *****************************************************************************
394 : : {
395 : 0 : m_maxit = maxit;
396 : 0 : m_tol = tol;
397 : 0 : m_solved = c;
398 : 0 : m_it = 0;
399 : :
400 : 0 : next();
401 : 0 : }
402 : :
403 : : void
404 : 0 : BiCG::next()
405 : : // *****************************************************************************
406 : : // Start next linear solver iteration
407 : : // *****************************************************************************
408 : : {
409 [ - - ]: 0 : if (m_it == 0) {
410 : 0 : m_alpha = 0.0;
411 : 0 : m_omega=0.0;
412 : :
413 : : }else{
414 : 0 : m_alpha =( m_rho/m_rho0 ) * ( m_alpha/m_omega ) ; //alpha functions as Beta??
415 : : }
416 : 0 : m_rho0 = m_rho;
417 : :
418 : : // compute p = r + alpha * p
419 [ - - ]: 0 : for (std::size_t i=0; i<m_p.size(); ++i) m_p[i] = m_r[i] + m_alpha * ( m_p[i] - m_omega * m_q[i] );
420 : :
421 : : // initiate computing q = A * p
422 [ - - ]: 0 : thisProxy[ thisIndex ].wait4q();
423 : 0 : qAp();
424 : 0 : }
425 : :
426 : : void
427 : 0 : BiCG::qAp()
428 : : // *****************************************************************************
429 : : // Initiate computing q = A * p
430 : : // *****************************************************************************
431 : : {
432 : : // Compute own contribution to q = A * p
433 : 0 : m_A.mult( m_p, m_q, m_bcmask );
434 : :
435 : : // Send partial product on chare-boundary nodes to fellow chares
436 [ - - ]: 0 : if (m_nodeCommMap.empty()) {
437 : 0 : comq_complete();
438 : : } else {
439 : : auto ncomp = m_A.Ncomp();
440 [ - - ]: 0 : for (const auto& [c,n] : m_nodeCommMap) {
441 : 0 : std::vector< std::vector< tk::real > > qc( n.size() );
442 : : std::size_t j = 0;
443 [ - - ]: 0 : for (auto g : n) {
444 [ - - ]: 0 : std::vector< tk::real > nq( ncomp );
445 : 0 : auto i = tk::cref_find( m_lid, g );
446 [ - - ]: 0 : for (std::size_t d=0; d<ncomp; ++d) nq[d] = m_q[ i*ncomp+d ];
447 : 0 : qc[j++] = std::move(nq);
448 : : }
449 [ - - ][ - - ]: 0 : thisProxy[c].comq( std::vector<std::size_t>(begin(n),end(n)), qc );
[ - - ][ - - ]
[ - - ]
450 : 0 : }
451 : : }
452 : :
453 : 0 : ownq_complete();
454 : 0 : }
455 : :
456 : : void
457 : 0 : BiCG::comq( const std::vector< std::size_t >& gid,
458 : : const std::vector< std::vector< tk::real > >& qc )
459 : : // *****************************************************************************
460 : : // Receive contributions to q = A * p on chare-boundaries
461 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
462 : : //! \param[in] qc Partial contributions at chare-boundary nodes
463 : : // *****************************************************************************
464 : : {
465 : : Assert( qc.size() == gid.size(), "Size mismatch" );
466 : :
467 : : using tk::operator+=;
468 : :
469 [ - - ]: 0 : for (std::size_t i=0; i<gid.size(); ++i)
470 : 0 : m_qc[ gid[i] ] += qc[i];
471 : :
472 [ - - ]: 0 : if (++m_nq == m_nodeCommMap.size()) {
473 : 0 : m_nq = 0;
474 : 0 : comq_complete();
475 : : }
476 : 0 : }
477 : :
478 : : void
479 : 0 : BiCG::q()
480 : : // *****************************************************************************
481 : : // Finish computing q = A * p
482 : : // *****************************************************************************
483 : : {
484 : : // Combine own and communicated contributions to q = A * p
485 : : auto ncomp = m_A.Ncomp();
486 [ - - ]: 0 : for (const auto& [gid,q] : m_qc) {
487 : 0 : auto i = tk::cref_find( m_lid, gid );
488 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
489 : 0 : m_q[i*ncomp+c] += q[c];
490 : : }
491 : 0 : tk::destroy( m_qc );
492 : :
493 : : //BiCGStab uses (rhat_0,q), initiate here
494 [ - - ]: 0 : dot( m_r0, m_q,
495 : 0 : CkCallback( CkReductionTarget(BiCG,pq), thisProxy ) );
496 : 0 : }
497 : :
498 : : void
499 [ - - ]: 0 : BiCG::pq( tk::real d )
500 : : // *****************************************************************************
501 : : // Compute the dot product (p,q)
502 : : //! \param[in] d Dot product of (p,q) (aggregated across all chares)
503 : : // *****************************************************************************
504 : : {
505 : : // If (p,q)=0, then p and q are orthogonal and the system either has a trivial
506 : : // solution, x=x0, or the BCs are incomplete or wrong, in either case the
507 : : // solve cannot continue.
508 : : // sod_pe4 DOES have a trivial solution, yet we are using it for a test case...
509 : : // we have to tolerate small eps values and rely on the numerator also being small.
510 : : // We should generally allow for the mesh to stay still though, so need to think about this some.
511 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
512 [ - - ]: 0 : if (std::abs(d) < eps) {
513 [ - - ]: 0 : if ( m_it > 0) {
514 : 0 : m_alpha = m_rho / d;
515 : : } else {
516 : 0 : m_it = m_maxit;
517 : 0 : m_alpha = 0.0;
518 : : }
519 : : } else {
520 : 0 : m_alpha = m_rho / d;
521 : : }
522 : :
523 : : // compute s = r - alpha * q
524 [ - - ]: 0 : for (std::size_t i=0; i<m_r.size(); ++i) m_r[i] -= m_alpha * m_q[i];
525 : :
526 : : // initiate computing norm of residual: (r,r)
527 [ - - ]: 0 : dot( m_r, m_r,
528 : 0 : CkCallback( CkReductionTarget(BiCG,normres), thisProxy ) );
529 : 0 : }
530 : :
531 : : void
532 : 0 : BiCG::normres( tk::real r )
533 : : // *****************************************************************************
534 : : // Compute norm of residual: (r,r)
535 : : //! \param[in] r Dot product, (r,r) (aggregated across all chares)
536 : : // *****************************************************************************
537 : : {
538 : 0 : m_rho = r;
539 : :
540 : : // Advance solution: h = x + alpha * p
541 [ - - ]: 0 : for (std::size_t i=0; i<m_x.size(); ++i) m_x[i] += m_alpha * m_p[i];
542 : :
543 : : // Communicate solution
544 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4x();
545 : :
546 : : // Send solution on chare-boundary nodes to fellow chares
547 [ - - ]: 0 : if (m_nodeCommMap.empty()) {
548 : 0 : comx_complete();
549 : : } else {
550 : : auto ncomp = m_A.Ncomp();
551 [ - - ]: 0 : for (const auto& [c,n] : m_nodeCommMap) {
552 : 0 : std::vector< std::vector< tk::real > > xc( n.size() );
553 : : std::size_t j = 0;
554 [ - - ]: 0 : for (auto g : n) {
555 [ - - ]: 0 : std::vector< tk::real > nx( ncomp );
556 : 0 : auto i = tk::cref_find( m_lid, g );
557 [ - - ]: 0 : for (std::size_t d=0; d<ncomp; ++d) nx[d] = m_x[ i*ncomp+d ];
558 : 0 : xc[j++] = std::move(nx);
559 : : }
560 [ - - ][ - - ]: 0 : thisProxy[c].comx( std::vector<std::size_t>(begin(n),end(n)), xc );
[ - - ][ - - ]
[ - - ]
561 : 0 : }
562 : : }
563 : :
564 : 0 : ownx_complete();
565 : 0 : }
566 : :
567 : : void
568 : 0 : BiCG::comx( const std::vector< std::size_t >& gid,
569 : : const std::vector< std::vector< tk::real > >& xc )
570 : : // *****************************************************************************
571 : : // Receive contributions to final solution on chare-boundaries
572 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
573 : : //! \param[in] xc Partial contributions at chare-boundary nodes
574 : : // *****************************************************************************
575 : : {
576 : : Assert( xc.size() == gid.size(), "Size mismatch" );
577 : :
578 [ - - ]: 0 : for (std::size_t i=0; i<gid.size(); ++i) m_xc[ gid[i] ] += xc[i];
579 : :
580 [ - - ]: 0 : if (++m_nx == m_nodeCommMap.size()) {
581 : 0 : m_nx = 0;
582 : 0 : comx_complete();
583 : : }
584 : 0 : }
585 : :
586 : : void
587 : 0 : BiCG::x()
588 : : // *****************************************************************************
589 : : // Assemble solution on chare boundaries
590 : : // *****************************************************************************
591 : : {
592 : : // Assemble solution on chare boundaries by averaging
593 : : auto ncomp = m_A.Ncomp();
594 [ - - ]: 0 : for (const auto& [g,x] : m_xc) {
595 : 0 : auto i = tk::cref_find(m_lid,g);
596 [ - - ]: 0 : for (std::size_t d=0; d<ncomp; ++d) m_x[i*ncomp+d] += x[d];
597 : 0 : auto c = tk::count(m_nodeCommMap,g);
598 [ - - ]: 0 : for (std::size_t d=0; d<ncomp; ++d) m_x[i*ncomp+d] /= c;
599 : : }
600 : 0 : tk::destroy( m_xc );
601 : :
602 : : //Don't iterate counter yet!
603 [ - - ]: 0 : auto normb = m_normb > 1.0e-14 ? m_normb : 1.0;
604 : 0 : auto normr = std::sqrt( m_rho );
605 [ - - ][ - - ]: 0 : if ( m_it < m_maxit && normr > m_tol*normb ) { //If we are not solved, continue, else exit
606 : : //Here is where we significantly diverge from regular CG as we dont' go to "next()" yet
607 : : // initiate computing t = A * s ;
608 : :
609 [ - - ]: 0 : thisProxy[ thisIndex ].wait4t();
610 : 0 : tAs();
611 : :
612 : : } else {
613 : :
614 [ - - ][ - - ]: 0 : m_converged = m_it == m_maxit && normr > m_tol*normb ? false : true;
615 : 0 : m_solved.send( CkDataMsg::buildNew( sizeof(tk::real), &normr ) );
616 : :
617 : : }
618 : 0 : }
619 : :
620 : : void
621 : 0 : BiCG::tAs()
622 : : // *****************************************************************************
623 : : // Initiate computing t = A * s
624 : : // *****************************************************************************
625 : : {
626 : : // Compute own contribution to t = A * s
627 : 0 : m_A.mult( m_r, m_t, m_bcmask );
628 : :
629 : : // Send partial product on chare-boundary nodes to fellow chares
630 [ - - ]: 0 : if (m_nodeCommMap.empty()) {
631 : 0 : comt_complete();
632 : : } else {
633 : : auto ncomp = m_A.Ncomp();
634 [ - - ]: 0 : for (const auto& [c,n] : m_nodeCommMap) {
635 : 0 : std::vector< std::vector< tk::real > > tc( n.size() );
636 : : std::size_t j = 0;
637 [ - - ]: 0 : for (auto g : n) {
638 [ - - ]: 0 : std::vector< tk::real > nt( ncomp );
639 : 0 : auto i = tk::cref_find( m_lid, g );
640 [ - - ]: 0 : for (std::size_t d=0; d<ncomp; ++d) nt[d] = m_t[ i*ncomp+d ];
641 : 0 : tc[j++] = std::move(nt);
642 : : }
643 [ - - ][ - - ]: 0 : thisProxy[c].comt( std::vector<std::size_t>(begin(n),end(n)), tc );
[ - - ][ - - ]
[ - - ]
644 : 0 : }
645 : : }
646 : :
647 : 0 : ownt_complete();
648 : 0 : }
649 : :
650 : : void
651 : 0 : BiCG::comt( const std::vector< std::size_t >& gid,
652 : : const std::vector< std::vector< tk::real > >& tc )
653 : : // *****************************************************************************
654 : : // Receive contributions to t = A * s on chare-boundaries
655 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
656 : : //! \param[in] tc Partial contributions at chare-boundary nodes
657 : : // *****************************************************************************
658 : : {
659 : : Assert( tc.size() == gid.size(), "Size mismatch" );
660 : :
661 : : using tk::operator+=;
662 : :
663 [ - - ]: 0 : for (std::size_t i=0; i<gid.size(); ++i)
664 : 0 : m_tc[ gid[i] ] += tc[i];
665 : :
666 [ - - ]: 0 : if (++m_nt == m_nodeCommMap.size()) {
667 : 0 : m_nt = 0;
668 : 0 : comt_complete();
669 : : }
670 : 0 : }
671 : :
672 : :
673 : : void
674 : 0 : BiCG::t()
675 : : // *****************************************************************************
676 : : // Finish computing t = A * s
677 : : // *****************************************************************************
678 : : {
679 : : // Combine own and communicated contributions to t = A * s
680 : : auto ncomp = m_A.Ncomp();
681 [ - - ]: 0 : for (const auto& [gid,t] : m_tc) {
682 : 0 : auto i = tk::cref_find( m_lid, gid );
683 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
684 : 0 : m_t[i*ncomp+c] += t[c];
685 : : }
686 : 0 : tk::destroy( m_tc );
687 : :
688 : : //omega = (t,s)/(t,t). Compute numerator
689 [ - - ]: 0 : dot( m_t, m_r,
690 : 0 : CkCallback( CkReductionTarget(BiCG,ts), thisProxy ) );
691 : 0 : }
692 : : void
693 : 0 : BiCG::ts( tk::real d )
694 : : // *****************************************************************************
695 : : // Compute the dot product (t,s)
696 : : //! \param[in] d Dot product of (t,s) (aggregated across all chares)
697 : : // *****************************************************************************
698 : : {
699 : :
700 : 0 : m_omega = d; //numerator set
701 : :
702 [ - - ]: 0 : dot( m_t, m_t, //compute denominator
703 : 0 : CkCallback( CkReductionTarget(BiCG,tt), thisProxy ) );
704 : :
705 : 0 : }
706 : : void
707 [ - - ]: 0 : BiCG::tt( tk::real d )
708 : : // *****************************************************************************
709 : : // Compute the dot product (t,t)
710 : : //! \param[in] d Dot product of (t,t) (aggregated across all chares)
711 : : // *****************************************************************************
712 : : {
713 : : // If (t,t)=0, then p and q are orthogonal and the system either has a trivial
714 : : // solution, x=x0, or the BCs are incomplete or wrong, in either case the
715 : : // solve cannot continue.
716 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
717 [ - - ][ - - ]: 0 : if (std::abs(d) < eps && std::abs(m_omega) > 1000*eps) { //Unsure how to check for smallness here since m_omega may also be very small
718 : : //m_it = m_maxit;
719 : 0 : m_omega = m_omega/d; //d is denominator
720 : : //m_omega = 0.0;
721 : : }
722 : : else {
723 : 0 : m_omega = m_omega/d; //d is denominator
724 : : }
725 : : //compute x = h + omega * s
726 [ - - ]: 0 : for (std::size_t i=0; i<m_x.size(); ++i) m_x[i] += m_omega * m_r[i];
727 : : // Now update r: compute r = s - omega * t
728 [ - - ]: 0 : for (std::size_t i=0; i<m_r.size(); ++i) m_r[i] -= m_omega * m_t[i];
729 : : // initiate computing norm of residual: (r,r) for exit criteria
730 [ - - ]: 0 : dot( m_r, m_r,
731 : 0 : CkCallback( CkReductionTarget(BiCG,normresomega), thisProxy ) );
732 : 0 : }
733 : : void
734 : 0 : BiCG::normresomega( tk::real r )
735 : : // *****************************************************************************
736 : : // Compute norm of residual: (r,r) for second half of bicgstab
737 : : //! \param[in] r Dot product, (r,r) (aggregated across all chares)
738 : : // *****************************************************************************
739 : : {
740 : 0 : m_rho = r;
741 : :
742 : : // Communicate solution
743 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4x2();
744 : :
745 : : // Send solution on chare-boundary nodes to fellow chares
746 [ - - ]: 0 : if (m_nodeCommMap.empty()) {
747 : 0 : comx2_complete(); //Does this need to be unique?
748 : : } else {
749 : : auto ncomp = m_A.Ncomp();
750 [ - - ]: 0 : for (const auto& [c,n] : m_nodeCommMap) {
751 : 0 : std::vector< std::vector< tk::real > > x2c( n.size() );
752 : : std::size_t j = 0;
753 [ - - ]: 0 : for (auto g : n) {
754 [ - - ]: 0 : std::vector< tk::real > nx2( ncomp );
755 : 0 : auto i = tk::cref_find( m_lid, g );
756 [ - - ]: 0 : for (std::size_t d=0; d<ncomp; ++d) nx2[d] = m_x[ i*ncomp+d ];
757 : 0 : x2c[j++] = std::move(nx2);
758 : : }
759 [ - - ][ - - ]: 0 : thisProxy[c].comx2( std::vector<std::size_t>(begin(n),end(n)), x2c );
[ - - ][ - - ]
[ - - ]
760 : 0 : }
761 : : }
762 : :
763 : 0 : ownx2_complete();
764 : 0 : }
765 : : void
766 : 0 : BiCG::x2()
767 : : // *****************************************************************************
768 : : // Assemble solution on chare boundaries
769 : : // *****************************************************************************
770 : : {
771 : : // Assemble solution on chare boundaries by averaging
772 : : auto ncomp = m_A.Ncomp();
773 [ - - ]: 0 : for (const auto& [g,x] : m_x2c) {
774 : 0 : auto i = tk::cref_find(m_lid,g);
775 [ - - ]: 0 : for (std::size_t d=0; d<ncomp; ++d) m_x[i*ncomp+d] += x[d];
776 : 0 : auto c = tk::count(m_nodeCommMap,g);
777 [ - - ]: 0 : for (std::size_t d=0; d<ncomp; ++d) m_x[i*ncomp+d] /= c;
778 : : }
779 : 0 : tk::destroy( m_x2c );
780 : :
781 : 0 : ++m_it;
782 [ - - ]: 0 : auto normb = m_normb > 1.0e-14 ? m_normb : 1.0;
783 : 0 : auto normr = std::sqrt( m_rho );
784 : :
785 [ - - ][ - - ]: 0 : if ( m_it < m_maxit && normr > m_tol*normb ) { //If we are not solved, continue, else exit
786 : : //get next rho value
787 [ - - ]: 0 : dot( m_r0, m_r,
788 : 0 : CkCallback( CkReductionTarget(BiCG,comprho), thisProxy ) );
789 : :
790 : : } else {
791 : :
792 [ - - ][ - - ]: 0 : m_converged = m_it == m_maxit && normr > m_tol*normb ? false : true;
793 : 0 : m_solved.send( CkDataMsg::buildNew( sizeof(tk::real), &normr ) );
794 : :
795 : : }
796 : 0 : }
797 : :
798 : : void
799 : 0 : BiCG::comprho( tk::real r )
800 : : {
801 : 0 : m_rho = r;
802 : 0 : next();
803 : 0 : }
804 : :
805 : : void
806 : 0 : BiCG::comx2( const std::vector< std::size_t >& gid,
807 : : const std::vector< std::vector< tk::real > >& x2c )
808 : : // *****************************************************************************
809 : : // Receive contributions to final solution on chare-boundaries
810 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
811 : : //! \param[in] xc Partial contributions at chare-boundary nodes
812 : : // *****************************************************************************
813 : : {
814 : : Assert( x2c.size() == gid.size(), "Size mismatch" );
815 : :
816 [ - - ]: 0 : for (std::size_t i=0; i<gid.size(); ++i) m_x2c[ gid[i] ] += x2c[i];
817 : :
818 [ - - ]: 0 : if (++m_nx2 == m_nodeCommMap.size()) {
819 : 0 : m_nx2 = 0;
820 : 0 : comx2_complete();
821 : : }
822 : 0 : }
823 : :
824 : : #include "NoWarning/bicg.def.h"
|
