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