Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/LinearSolver/ConjugateGradients.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 ConjugateGradients 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 ConjugateGradients_h
22 : : #define ConjugateGradients_h
23 : :
24 : : #include "Types.hpp"
25 : : #include "CSR.hpp"
26 : :
27 : : #include "NoWarning/conjugategradients.decl.h"
28 : :
29 : : namespace tk {
30 : :
31 : : //! \brief ConjugateGradients Charm++ chare array used to perform a distributed
32 : : //! linear solve with the conjugate gradients algorithm
33 : : class ConjugateGradients : public CBase_ConjugateGradients {
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 : : ConjugateGradients_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 ConjugateGradients(
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 ConjugateGradients(
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 : 21 : const NodeCommMap& nodecommmap ) :
73 : : ConjugateGradients( std::move(std::get<0>(system)),
74 : : std::move(std::get<1>(system)),
75 : : std::move(std::get<2>(system)),
76 [ - - ][ + - ]: 21 : gid, lid, nodecommmap ) {}
77 : :
78 : : //! Migrate constructor
79 : 57 : explicit ConjugateGradients( 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 : :
119 : : //! Compute the dot product (p,q)
120 : : void pq( tk::real d );
121 : :
122 : : //! Compute the norm of the residual: (r,r)
123 : : void normres( tk::real r );
124 : :
125 : : //! Access solution
126 [ + - ][ + - ]: 801 : std::vector< tk::real > solution() const { return m_x; }
127 : :
128 : : //! Return convergence flag
129 : : bool converged() const { return m_converged; }
130 : :
131 : : /** @name Pack/unpack (Charm++ serialization) routines */
132 : : ///@{
133 : : //! \brief Pack/Unpack serialize member function
134 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
135 : 192 : void pup( PUP::er &p ) override {
136 : 192 : p | m_A;
137 : 192 : p | m_x;
138 : 192 : p | m_b;
139 : 192 : p | m_gid;
140 : 192 : p | m_lid;
141 : 192 : p | m_nodeCommMap;
142 : 192 : p | m_r;
143 : 192 : p | m_rc;
144 : 192 : p | m_nr;
145 : 192 : p | m_bc;
146 : 192 : p | m_bcc;
147 : 192 : p | m_bcmask;
148 : 192 : p | m_nb;
149 : 192 : p | m_p;
150 : 192 : p | m_q;
151 : 192 : p | m_qc;
152 : 192 : p | m_nq;
153 : 192 : p | m_initres;
154 : 192 : p | m_solved;
155 : 192 : p | m_normb;
156 : 192 : p | m_it;
157 : 192 : p | m_maxit;
158 : 192 : p | m_tol;
159 : 192 : p | m_rho;
160 : 192 : p | m_rho0;
161 : 192 : p | m_alpha;
162 : 192 : p | m_converged;
163 : 192 : p | m_xc;
164 : 192 : p | m_nx;
165 : 192 : }
166 : : //! \brief Pack/Unpack serialize operator|
167 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
168 : : //! \param[in,out] c ConjugateGradients object reference
169 : : friend void operator|( PUP::er& p, ConjugateGradients& c ) { c.pup(p); }
170 : : ///@}
171 : :
172 : : private:
173 : : //! Sparse matrix
174 : : CSR m_A;
175 : : //! Solution/unknown
176 : : std::vector< tk::real > m_x;
177 : : //! Right hand side
178 : : std::vector< tk::real > m_b;
179 : : //! Global node IDs
180 : : std::vector< std::size_t > m_gid;
181 : : //! Local node IDs associated to global ones
182 : : std::unordered_map< std::size_t, std::size_t > m_lid;
183 : : //! Global mesh node IDs shared with other chares associated to chare IDs
184 : : NodeCommMap m_nodeCommMap;
185 : : //! Auxiliary vector for CG solve
186 : : std::vector< tk::real > m_r;
187 : : //! Receive buffer for communication of r = b - A * x
188 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_rc;
189 : : //! Counter for assembling m_r
190 : : std::size_t m_nr;
191 : : //! Dirichlet boundary conditions
192 : : std::unordered_map< std::size_t,
193 : : std::vector< std::pair< bool, tk::real > > > m_bc;
194 : : //! Dirichlet boundary conditions communication buffer
195 : : std::unordered_map< std::size_t,
196 : : std::vector< std::pair< bool, tk::real > > > m_bcc;
197 : : //! Dirichlet boundary condition mask
198 : : std::vector< tk::real > m_bcmask;
199 : : //! Counter for assembling boundary conditions
200 : : std::size_t m_nb;
201 : : //! Auxiliary vector for CG solve
202 : : std::vector< tk::real > m_p;
203 : : //! Auxiliary vector for CG solve
204 : : std::vector< tk::real > m_q;
205 : : //! Receive buffer for communication of q = A * p
206 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_qc;
207 : : //! Counter for assembling m_q
208 : : std::size_t m_nq;
209 : : //! Charm++ callback to continue with when the setup is complete
210 : : CkCallback m_initres;
211 : : //! Charm++ callback to continue with when the solve is complete
212 : : CkCallback m_solved;
213 : : //! L2 norm of the right hand side
214 : : tk::real m_normb;
215 : : //! Iteration count
216 : : std::size_t m_it;
217 : : //! Max iteration count
218 : : std::size_t m_maxit;
219 : : //! Stop tolerance
220 : : tk::real m_tol;
221 : : //! Helper scalar for CG algorithm
222 : : tk::real m_rho;
223 : : //! Helper scalar for CG algorithm
224 : : tk::real m_rho0;
225 : : //! Helper scalar for CG algorithm
226 : : tk::real m_alpha;
227 : : //! Convergence flag: true if linear smoother converged to tolerance
228 : : bool m_converged;
229 : : //! Receive buffer for solution
230 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_xc;
231 : : //! Counter for assembling the solution on chare boundaries
232 : : std::size_t m_nx;
233 : :
234 : : //! Initiate computationa of dot product of two vectors
235 : : void dot( const std::vector< tk::real >& a,
236 : : const std::vector< tk::real >& b,
237 : : CkCallback c );
238 : :
239 : : //! Initiate A * x for computing the residual, r = b - A * x
240 : : void residual();
241 : : //! Finish computing the initial residual, r = b - A * x
242 : : void initres();
243 : :
244 : : //! Apply boundary conditions
245 : : void apply( CkCallback cb );
246 : :
247 : : //! Initiate computing q = A * p
248 : : void qAp();
249 : : //! Finish computing q = A * p
250 : : void q();
251 : :
252 : : //! Start next linear solver iteration
253 : : void next();
254 : :
255 : : //! Assemble solution on chare boundaries
256 : : void x();
257 : : };
258 : :
259 : : } // tk::
260 : :
261 : : #endif // ConjugateGradients_h
|
