Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/LinearSolver/BiCG.hpp
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 : :
12 : : There are a potentially large number of BiCG Charm++ chares.
13 : : Each ConjugateGradient chare gets a chunk of the full load, due to partiting
14 : : the mesh, on which the solve is performed.
15 : :
16 : : The implementation uses the Charm++ runtime system and is fully
17 : : asynchronous, overlapping computation and communication. The algorithm
18 : : utilizes the structured dagger (SDAG) Charm++ functionality.
19 : : */
20 : : // *****************************************************************************
21 : : #ifndef BiCG_h
22 : : #define BiCG_h
23 : :
24 : : #include "Types.hpp"
25 : : #include "CSR.hpp"
26 : :
27 : : #include "NoWarning/bicg.decl.h"
28 : :
29 : : namespace tk {
30 : :
31 : : //! \brief BiCG Charm++ chare array used to perform a distributed
32 : : //! linear solve with the conjugate gradients algorithm
33 : : class BiCG : public CBase_BiCG {
34 : :
35 : : public:
36 : : #if defined(__clang__)
37 : : #pragma clang diagnostic push
38 : : #pragma clang diagnostic ignored "-Wunused-parameter"
39 : : #elif defined(STRICT_GNUC)
40 : : #pragma GCC diagnostic push
41 : : #pragma GCC diagnostic ignored "-Wunused-parameter"
42 : : #endif
43 : : // Include Charm++ SDAG code. See http://charm.cs.illinois.edu/manuals/html/
44 : : // charm++/manual.html, Sec. "Structured Control Flow: Structured Dagger".
45 : : BiCG_SDAG_CODE
46 : : #if defined(__clang__)
47 : : #pragma clang diagnostic pop
48 : : #elif defined(STRICT_GNUC)
49 : : #pragma GCC diagnostic pop
50 : : #endif
51 : :
52 : : //! Constructor
53 : : explicit BiCG(
54 : : const CSR& A,
55 : : const std::vector< tk::real >& x,
56 : : const std::vector< tk::real >& b,
57 : : const std::vector< std::size_t >& gid,
58 : : const std::unordered_map< std::size_t, std::size_t >& lid,
59 : : const NodeCommMap& nodecommmap );
60 : :
61 : : #if defined(__clang__)
62 : : #pragma clang diagnostic push
63 : : #pragma clang diagnostic ignored "-Wundefined-func-template"
64 : : #endif
65 : : //! Constructor taking a tuple of {A,x,b} by rvalue reference
66 : : explicit BiCG(
67 : : std::tuple< tk::CSR,
68 : : std::vector< tk::real >,
69 : : std::vector< tk::real > >&& system,
70 : : const std::vector< std::size_t >& gid,
71 : : const std::unordered_map< std::size_t, std::size_t >& lid,
72 : 0 : const NodeCommMap& nodecommmap ) :
73 : : BiCG( std::move(std::get<0>(system)),
74 : : std::move(std::get<1>(system)),
75 : : std::move(std::get<2>(system)),
76 [ - - ][ - - ]: 0 : gid, lid, nodecommmap ) {}
77 : :
78 : : //! Migrate constructor
79 : 0 : explicit BiCG( CkMigrateMessage* ) {}
80 : : #if defined(__clang__)
81 : : #pragma clang diagnostic pop
82 : : #endif
83 : :
84 : : //! Solve linear system
85 : : void solve( std::size_t maxit, tk::real tol, CkCallback c );
86 : :
87 : : //! Initialize linear solve: set initial guess and boundary conditions
88 : : void init( const std::vector< tk::real >& x,
89 : : const std::vector< tk::real >& b,
90 : : const std::unordered_map< std::size_t,
91 : : std::vector< std::pair< bool, tk::real > > >& bc,
92 : : std::size_t ignorebc,
93 : : CkCallback cb );
94 : :
95 : : //! Setup solver
96 : : void setup( CkCallback c );
97 : :
98 : : //! Compute the norm of the right hand side
99 : : void normb( tk::real n );
100 : :
101 : : //! Compute rho = (r,r)
102 : : void rho( tk::real r );
103 : :
104 : : //! Receive contributions to r = b - A * x on chare-boundaries
105 : : void comres( const std::vector< std::size_t >& gid,
106 : : const std::vector< std::vector< tk::real > >& rc );
107 : :
108 : : //! Receive contributions to boundary conditions on chare-boundaries
109 : : void combc( const std::unordered_map< std::size_t,
110 : : std::vector< std::pair< bool, tk::real > > >& bc );
111 : :
112 : : //! Receive contributions to q = A * p on chare-boundaries
113 : : void comq( const std::vector< std::size_t >& gid,
114 : : const std::vector< std::vector< tk::real > >& qc );
115 : :
116 : : void comx( const std::vector< std::size_t >& gid,
117 : : const std::vector< std::vector< tk::real > >& xc );
118 : : void comt( const std::vector< std::size_t >& gid,
119 : : const std::vector< std::vector< tk::real > >& tc );
120 : : void comx2( const std::vector< std::size_t >& gid,
121 : : const std::vector< std::vector< tk::real > >& x2c );
122 : :
123 : : //! Compute the dot product (p,q)
124 : : void pq( tk::real d );
125 : : //! Compute the dot product (t,s)
126 : : void ts( tk::real d );
127 : : //! Compute the dot product (t,t)
128 : : void tt( tk::real d );
129 : :
130 : : //! Compute the norm of the residual: (r,r)
131 : : void normres( tk::real r );
132 : : //! Compute the norm of the residual: (r0,r)
133 : : void comprho( tk::real r );
134 : : //! Compute the norm of the residual 2nd step: (r,r)
135 : : void normresomega( tk::real r );
136 : :
137 : : //! Access solution
138 : : std::vector< tk::real > solution() const { return m_x; }
139 : :
140 : : //! Return convergence flag
141 : : bool converged() const { return m_converged; }
142 : :
143 : : /** @name Pack/unpack (Charm++ serialization) routines */
144 : : ///@{
145 : : //! \brief Pack/Unpack serialize member function
146 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
147 : 0 : void pup( PUP::er &p ) override {
148 : 0 : p | m_A;
149 : 0 : p | m_x;
150 : 0 : p | m_b;
151 : 0 : p | m_gid;
152 : 0 : p | m_lid;
153 : 0 : p | m_nodeCommMap;
154 : 0 : p | m_r;
155 : 0 : p | m_r0;
156 : 0 : p | m_rc;
157 : 0 : p | m_nr;
158 : 0 : p | m_bc;
159 : 0 : p | m_bcc;
160 : 0 : p | m_bcmask;
161 : 0 : p | m_nb;
162 : 0 : p | m_p;
163 : 0 : p | m_q;
164 : 0 : p | m_t;
165 : 0 : p | m_qc;
166 : 0 : p | m_tc;
167 : 0 : p | m_nt;
168 : 0 : p | m_nq;
169 : 0 : p | m_initres;
170 : 0 : p | m_solved;
171 : 0 : p | m_normb;
172 : 0 : p | m_it;
173 : 0 : p | m_maxit;
174 : 0 : p | m_tol;
175 : 0 : p | m_rho;
176 : 0 : p | m_rho0;
177 : 0 : p | m_alpha;
178 : 0 : p | m_omega;
179 : 0 : p | m_converged;
180 : 0 : p | m_xc;
181 : 0 : p | m_x2c;
182 : 0 : p | m_nx;
183 : 0 : p | m_nx2;
184 : 0 : }
185 : : //! \brief Pack/Unpack serialize operator|
186 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
187 : : //! \param[in,out] c BiCG object reference
188 : : friend void operator|( PUP::er& p, BiCG& c ) { c.pup(p); }
189 : : ///@}
190 : :
191 : : private:
192 : : //! Sparse matrix
193 : : CSR m_A;
194 : : //! Solution/unknown
195 : : std::vector< tk::real > m_x;
196 : : //! Right hand side
197 : : std::vector< tk::real > m_b;
198 : : //! Global node IDs
199 : : std::vector< std::size_t > m_gid;
200 : : //! Local node IDs associated to global ones
201 : : std::unordered_map< std::size_t, std::size_t > m_lid;
202 : : //! Global mesh node IDs shared with other chares associated to chare IDs
203 : : NodeCommMap m_nodeCommMap;
204 : : //! Auxiliary vector for CG solve
205 : : std::vector< tk::real > m_r;
206 : : //! Auxiliary vector for CG solve
207 : : std::vector< tk::real > m_r0;
208 : : //! Receive buffer for communication of r = b - A * x
209 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_rc;
210 : : //! Counter for assembling m_r
211 : : std::size_t m_nr;
212 : : //! Dirichlet boundary conditions
213 : : std::unordered_map< std::size_t,
214 : : std::vector< std::pair< bool, tk::real > > > m_bc;
215 : : //! Dirichlet boundary conditions communication buffer
216 : : std::unordered_map< std::size_t,
217 : : std::vector< std::pair< bool, tk::real > > > m_bcc;
218 : : //! Dirichlet boundary condition mask
219 : : std::vector< tk::real > m_bcmask;
220 : : //! Counter for assembling boundary conditions
221 : : std::size_t m_nb;
222 : : //! Auxiliary vector for CG solve
223 : : std::vector< tk::real > m_p;
224 : : //! Auxiliary vector for CG solve
225 : : std::vector< tk::real > m_q;
226 : : //! Auxiliary vector for the CG solve
227 : : std::vector< tk::real > m_t;
228 : : //! Receive buffer for communication of q = A * p
229 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_qc;
230 : : //! Receive buffer for communication of t = A * s
231 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_tc;
232 : : //! Counter for assembling m_t
233 : : std::size_t m_nt;
234 : : //! Counter for assembling m_q
235 : : std::size_t m_nq;
236 : : //! Charm++ callback to continue with when the setup is complete
237 : : CkCallback m_initres;
238 : : //! Charm++ callback to continue with when the solve is complete
239 : : CkCallback m_solved;
240 : : //! L2 norm of the right hand side
241 : : tk::real m_normb;
242 : : //! Iteration count
243 : : std::size_t m_it;
244 : : //! Max iteration count
245 : : std::size_t m_maxit;
246 : : //! Stop tolerance
247 : : tk::real m_tol;
248 : : //! Helper scalar for CG algorithm
249 : : tk::real m_rho;
250 : : //! Helper scalar for CG algorithm
251 : : tk::real m_rho0;
252 : : //! Helper scalar for CG algorithm
253 : : tk::real m_alpha;
254 : : //! Helper scalar for CG algorithm
255 : : tk::real m_omega;
256 : : //! Convergence flag: true if linear smoother converged to tolerance
257 : : bool m_converged;
258 : : //! Receive buffer for solution
259 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_xc;
260 : : //! Receive buffer for solution 2
261 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_x2c;
262 : : //! Counter for assembling the solution on chare boundaries
263 : : std::size_t m_nx;
264 : : //! Counter for assembling the solution 2 on chare boundaries
265 : : std::size_t m_nx2;
266 : :
267 : : //! Initiate computationa of dot product of two vectors
268 : : void dot( const std::vector< tk::real >& a,
269 : : const std::vector< tk::real >& b,
270 : : CkCallback c );
271 : :
272 : : //! Initiate A * x for computing the residual, r = b - A * x
273 : : void residual();
274 : : //! Finish computing the initial residual, r = b - A * x
275 : : void initres();
276 : :
277 : : //! Apply boundary conditions
278 : : void apply( CkCallback cb );
279 : :
280 : : //! Initiate computing q = A * p
281 : : void qAp();
282 : : //! Finish computing q = A * p
283 : : void q();
284 : :
285 : : //! Start next linear solver iteration
286 : : void next();
287 : :
288 : : //! Assemble solution on chare boundaries
289 : : void x();
290 : :
291 : : // Initiate computing t = A * s
292 : : void tAs();
293 : : //! Finish computing t = A * s
294 : : void t();
295 : :
296 : : //Assemble solution on chare boundaries
297 : : void x2();
298 : : };
299 : :
300 : : } // tk::
301 : :
302 : : #endif // BiCG_h
|
