Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Base/Data.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 Generic data storage with different memory layouts
9 : : \details Generic data storage with different memory layouts. See also the
10 : : rationale discussed in the [design](layout.html) document.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef Data_h
14 : : #define Data_h
15 : :
16 : : #include <array>
17 : : #include <string>
18 : : #include <cstdint>
19 : : #include <vector>
20 : : #include <set>
21 : : #include <algorithm>
22 : :
23 : : #include "Types.hpp"
24 : : #include "Exception.hpp"
25 : :
26 : : #include "NoWarning/pup_stl.hpp"
27 : :
28 : : namespace tk {
29 : :
30 : : //! Tags for selecting data layout policies
31 : : const uint8_t UnkEqComp = 0;
32 : : const uint8_t EqCompUnk = 1;
33 : :
34 : : //! Zero-runtime-cost data-layout wrappers with type-based compile-time dispatch
35 : : template< uint8_t Layout >
36 [ + - ][ + - ]: 248427 : class Data {
[ + - ][ + - ]
[ + + ][ + + ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ - + ][ - + ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ - + ][ - + ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
37 : :
38 : : private:
39 : : //! \brief Inherit type of number of components from keyword 'ncomp'
40 : : using ncomp_t = tk::ncomp_t;
41 : :
42 : : public:
43 : : //! Default constructor (required for Charm++ migration)
44 [ + - ][ - - ]: 28878 : explicit Data() : m_vec(), m_nunk(), m_nprop() {}
45 : :
46 : : //! Constructor
47 : : //! \param[in] nu Number of unknowns to allocate memory for
48 : : //! \param[in] np Total number of properties, i.e., scalar variables or
49 : : //! components, per unknown
50 [ + - ][ + - ]: 7054 : explicit Data( ncomp_t nu, ncomp_t np ) :
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
51 : : m_vec( nu*np ),
52 : : m_nunk( nu ),
53 [ + - ][ + + ]: 2110581 : m_nprop( np ) {}
[ + + ][ + + ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
54 : :
55 : : //! Const data access dispatch
56 : : //! \details Public interface to const-ref data access to a single real
57 : : //! value. Use it as Data(p,c), where p is the unknown index, and c is
58 : : //! the component index specifying the scalar equation within a system of
59 : : //! equations. Requirement: component <
60 : : //! nprop, unknown < nunk, enforced with an assert in DEBUG mode, see also
61 : : //! the constructor.
62 : : //! \param[in] unknown Unknown index
63 : : //! \param[in] component Component index, i.e., position of a scalar within
64 : : //! a system
65 : : //! \return Const reference to data of type tk::real
66 : : const tk::real&
67 : : operator()( ncomp_t unknown, ncomp_t component ) const
68 : : { return access( unknown, component, int2type< Layout >() ); }
69 : :
70 : : //! Non-const data access dispatch
71 : : //! \details Public interface to non-const-ref data access to a single real
72 : : //! value. Use it as Data(p,c), where p is the unknown index, and c is
73 : : //! the component index specifying the scalar equation within a system of
74 : : //! equations. Requirement: component <
75 : : //! nprop, unknown < nunk, enforced with an assert in DEBUG mode, see also
76 : : //! the constructor.
77 : : //! \param[in] unknown Unknown index
78 : : //! \param[in] component Component index, i.e., position of a scalar within
79 : : //! a system
80 : : //! \return Non-const reference to data of type tk::real
81 : : //! \see "Avoid Duplication in const and Non-const Member Function," and
82 : : //! "Use const whenever possible," Scott Meyers, Effective C++, 3d ed.
83 : : tk::real&
84 : : operator()( ncomp_t unknown, ncomp_t component ) {
85 : : return const_cast< tk::real& >(
86 : : static_cast< const Data& >( *this ).
87 : : operator()( unknown, component ) );
88 : : }
89 : :
90 : : //! Const ptr to physical variable access dispatch
91 : : //! \details Public interface to the first half of a physical variable
92 : : //! access. cptr() and var() are two member functions intended to be used
93 : : //! together in case when component would be expensive to
94 : : //! compute for data access via the function call operator, i.e., cptr()
95 : : //! can be used to pre-compute part of the address, which returns a
96 : : //! pointer and var() can be used to finish the data access using the
97 : : //! pointer returned by cptr(). In other words, cptr() returns part of the
98 : : //! address known based on component and intended to be used in
99 : : //! a setup phase. Then var() takes this partial address and finishes the
100 : : //! address calculation given the unknown id. Thus the following two data
101 : : //! accesses are equivalent (modulo constness):
102 : : //! * real& value = operator()( unk, comp ); and
103 : : //! * const real* p = cptr( comp ); and
104 : : //! const real& value = var( p, unk ); or real& value = var( p, unk );
105 : : //! Requirement: component < nprop, enforced with an assert in
106 : : //! DEBUG mode, see also the constructor.
107 : : //! \param[in] component Component index, i.e., position of a scalar within
108 : : //! a system
109 : : //! \return Pointer to data of type tk::real for use with var()
110 : : //! \see Example client code in Statistics::setupOrdinary() and
111 : : //! Statistics::accumulateOrd() in Statistics/Statistics.C.
112 : : const tk::real*
113 : : cptr( ncomp_t component ) const
114 : : { return cptr( component, int2type< Layout >() ); }
115 : :
116 : : //! Const-ref data-access dispatch
117 : : //! \details Public interface to the second half of a physical variable
118 : : //! access. cptr() and var() are two member functions intended to be used
119 : : //! together in case when component would be expensive to
120 : : //! compute for data access via the function call operator, i.e., cptr()
121 : : //! can be used to pre-compute part of the address, which returns a
122 : : //! pointer and var() can be used to finish the data access using the
123 : : //! pointer returned by cptr(). In other words, cptr() returns part of the
124 : : //! address known based on component and intended to be used in
125 : : //! a setup phase. Then var() takes this partial address and finishes the
126 : : //! address calculation given the unknown id. Thus the following two data
127 : : //! accesses are equivalent (modulo constness):
128 : : //! * real& value = operator()( unk, comp ); and
129 : : //! * const real* p = cptr( comp ); and
130 : : //! const real& value = var( p, unk ); or real& value = var( p, unk );
131 : : //! Requirement: unknown < nunk, enforced with an assert in DEBUG mode,
132 : : //! see also the constructor.
133 : : //! \param[in] pt Pointer to data of type tk::real as returned from cptr()
134 : : //! \param[in] unknown Unknown index
135 : : //! \return Const reference to data of type tk::real
136 : : //! \see Example client code in Statistics::setupOrdinary() and
137 : : //! Statistics::accumulateOrd() in Statistics/Statistics.C.
138 : : const tk::real&
139 : : var( const tk::real* pt, ncomp_t unknown ) const
140 : 1365923362 : { return var( pt, unknown, int2type< Layout >() ); }
141 : :
142 : : //! Non-const-ref data-access dispatch
143 : : //! \details Public interface to the second half of a physical variable
144 : : //! access. cptr() and var() are two member functions intended to be used
145 : : //! together in case when component would be expensive to
146 : : //! compute for data access via the function call operator, i.e., cptr()
147 : : //! can be used to pre-compute part of the address, which returns a
148 : : //! pointer and var() can be used to finish the data access using the
149 : : //! pointer returned by cptr(). In other words, cptr() returns part of the
150 : : //! address known based on component and intended to be used in
151 : : //! a setup phase. Then var() takes this partial address and finishes the
152 : : //! address calculation given the unknown id. Thus the following two data
153 : : //! accesses are equivalent (modulo constness):
154 : : //! * real& value = operator()( unk, comp ); and
155 : : //! * const real* p = cptr( comp ); and
156 : : //! const real& value = var( p, unk ); or real& value = var( p, unk );
157 : : //! Requirement: unknown < nunk, enforced with an assert in DEBUG mode,
158 : : //! see also the constructor.
159 : : //! \param[in] pt Pointer to data of type tk::real as returned from cptr()
160 : : //! \param[in] unknown Unknown index
161 : : //! \return Non-const reference to data of type tk::real
162 : : //! \see Example client code in Statistics::setupOrdinary() and
163 : : //! Statistics::accumulateOrd() in Statistics/Statistics.C.
164 : : //! \see "Avoid Duplication in const and Non-const Member Function," and
165 : : //! "Use const whenever possible," Scott Meyers, Effective C++, 3d ed.
166 : : tk::real&
167 : : var( const tk::real* pt, ncomp_t unknown ) {
168 : : return const_cast< tk::real& >(
169 : : static_cast< const Data& >( *this ).var( pt, unknown ) );
170 : : }
171 : :
172 : : //! Access to number of unknowns
173 : : //! \return Number of unknowns
174 : : ncomp_t nunk() const noexcept { return m_nunk; }
175 : :
176 : : //! Access to number of properties
177 : : //! \details This is the total number of scalar components per unknown
178 : : //! \return Number of propertes/unknown
179 : : ncomp_t nprop() const noexcept { return m_nprop; }
180 : :
181 : : //! Extract flat vector of all unknowns
182 : : //! \return Flat vector of reals
183 : : std::vector< tk::real >
184 : 679 : flat() const {
185 : 679 : std::vector< tk::real > w( m_nunk * m_nprop );
186 [ + + ]: 2716 : for (std::size_t j=0; j<m_nprop; ++j)
187 [ + + ]: 3754755 : for (std::size_t i=0; i<m_nunk; ++i)
188 : 3752718 : w[i*m_nprop+j] = operator()( i, j );
189 : 679 : return w;
190 : : }
191 : :
192 : : //! Assign from flat vector
193 : : //! \param[in] rhs Flat vector to assign from
194 : : //! \note This is the opposite of flat().
195 : 1 : void operator=( const std::vector< tk::real >& rhs ) {
196 : : Assert( rhs.size() == m_nunk * m_nprop, "Size mismatch" );
197 [ + + ][ + + ]: 8 : for (std::size_t j=0; j<m_nprop; ++j)
198 [ + + ][ + + ]: 18 : for (std::size_t i=0; i<m_nunk; ++i)
199 : 12 : operator()( i, j ) = rhs[i*m_nprop+j];
200 : 1 : }
201 : :
202 : : //! Assign from array(3) of vectors
203 : : //! \param[in] rhs Array(3) of vectors to assign from
204 : : void operator=( const std::array< std::vector< tk::real >, 3 >& rhs ) {
205 : : Assert( m_nprop == 3, "Size mismatch" );
206 : : Assert( rhs[0].size() == m_nunk, "Size mismatch" );
207 : : Assert( rhs[1].size() == m_nunk, "Size mismatch" );
208 : : Assert( rhs[2].size() == m_nunk, "Size mismatch" );
209 [ + + ][ + + ]: 8 : for (std::size_t j=0; j<3; ++j)
210 [ + + ][ + + ]: 18 : for (std::size_t i=0; i<m_nunk; ++i)
211 : 12 : operator()( i, j ) = rhs[j][i];
212 : : }
213 : :
214 : : //! Extract vector of unknowns given component
215 : : //! \details Requirement: component < nprop, enforced with an
216 : : //! assert in DEBUG mode, see also the constructor.
217 : : //! \param[in] component Component index, i.e., position of a scalar within
218 : : //! a system
219 : : //! \return A vector of unknowns given by component (length:
220 : : //! nunk(), i.e., the first constructor argument)
221 : : std::vector< tk::real >
222 : 99622 : extract_comp( ncomp_t component ) const {
223 : 99622 : std::vector< tk::real > w( m_nunk );
224 [ + + ]: 25194293 : for (ncomp_t i=0; i<m_nunk; ++i)
225 : 25094671 : w[i] = operator()( i, component );
226 : 99622 : return w;
227 : : }
228 : :
229 : : //! Extract (a copy of) all components for an unknown
230 : : //! \details Requirement: unknown < nunk, enforced with an assert in DEBUG
231 : : //! mode, see also the constructor.
232 : : //! \param[in] unknown Index of unknown
233 : : //! \return A vector of components for a single unknown (length: nprop,
234 : : //! i.e., the second constructor argument)
235 : : std::vector< tk::real >
236 : 67254236 : extract( ncomp_t unknown ) const {
237 : 67254236 : std::vector< tk::real > w( m_nprop );
238 [ + + ]: 667801057 : for (ncomp_t i=0; i<m_nprop; ++i) w[i] = operator()( unknown, i );
239 : 67254236 : return w;
240 : : }
241 : :
242 : : //! Extract all components for unknown
243 : : //! \details Requirement: unknown < nunk, enforced with an assert in DEBUG
244 : : //! mode, see also the constructor.
245 : : //! \param[in] unknown Index of unknown
246 : : //! \return A vector of components for a single unknown (length: nprop,
247 : : //! i.e., the second constructor argument)
248 : : //! \note This is simply an alias for extract( unknown )
249 : : std::vector< tk::real >
250 [ + - ][ + - ]: 66540870 : operator[]( ncomp_t unknown ) const { return extract( unknown ); }
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
251 : :
252 : : //! Extract (a copy of) four values of unknowns
253 : : //! \details Requirement: component < nprop, [A,B,C,D] < nunk,
254 : : //! enforced with an assert in DEBUG mode, see also the constructor.
255 : : //! \param[in] component Component index, i.e., position of a scalar within
256 : : //! a system
257 : : //! \param[in] A Index of 1st unknown
258 : : //! \param[in] B Index of 2nd unknown
259 : : //! \param[in] C Index of 3rd unknown
260 : : //! \param[in] D Index of 4th unknown
261 : : //! \return Array of the four values of component
262 : : std::array< tk::real, 4 >
263 : : extract( ncomp_t component,
264 : : ncomp_t A, ncomp_t B, ncomp_t C, ncomp_t D ) const
265 : : {
266 : : auto p = cptr( component );
267 [ + - ]: 28802926 : return {{ var(p,A), var(p,B), var(p,C), var(p,D) }};
268 : : }
269 : :
270 : : //! Extract (a copy of) four values of unknowns
271 : : //! \details Requirement: component < nprop, for all N[i] < nunk,
272 : : //! enforced with an assert in DEBUG mode, see also the constructor.
273 : : //! \param[in] component Component index, i.e., position of a scalar within
274 : : //! a system
275 : : //! \param[in] N Indices of the 4 unknowns
276 : : //! \return Array of the four values of component
277 : : std::array< tk::real, 4 >
278 : 28802929 : extract( ncomp_t component,
279 : : const std::array< ncomp_t, 4 >& N ) const
280 : : {
281 : 28802929 : return extract( component, N[0], N[1], N[2], N[3] );
282 : : }
283 : :
284 : : //! Extract (a copy of) three values of unknowns
285 : : //! \details Requirement: component < nprop, [A,B,C] < nunk,
286 : : //! enforced with an assert in DEBUG mode, see also the constructor.
287 : : //! \param[in] component Component index, i.e., position of a scalar within
288 : : //! a system
289 : : //! \param[in] A Index of 1st unknown
290 : : //! \param[in] B Index of 2nd unknown
291 : : //! \param[in] C Index of 3rd unknown
292 : : //! \return Array of the four values of component
293 : : std::array< tk::real, 3 >
294 : : extract( ncomp_t component,
295 : : ncomp_t A, ncomp_t B, ncomp_t C ) const
296 : : {
297 : : auto p = cptr( component );
298 [ + - ][ + - ]: 2 : return {{ var(p,A), var(p,B), var(p,C) }};
299 : : }
300 : :
301 : : //! Extract (a copy of) three values of unknowns
302 : : //! \details Requirement: component < nprop, for all N[i] < nunk,
303 : : //! enforced with an assert in DEBUG mode, see also the constructor.
304 : : //! \param[in] component Component index, i.e., position of a scalar within
305 : : //! a system
306 : : //! \param[in] N Indices of the 3 unknowns
307 : : //! \return Array of the three values of component
308 : : std::array< tk::real, 3 >
309 : : extract( ncomp_t component,
310 : : const std::array< ncomp_t, 3 >& N ) const
311 : : {
312 : : return extract( component, N[0], N[1], N[2] );
313 : : }
314 : :
315 : : //! Const-ref accessor to underlying raw data as a std::vector
316 : : //! \return Constant reference to underlying raw data
317 : : const std::vector< tk::real >& vec() const { return m_vec; }
318 : :
319 : : //! Non-const-ref accessor to underlying raw data as a std::vector
320 : : //! \return Non-constant reference to underlying raw data
321 : : std::vector< tk::real >& vec() { return m_vec; }
322 : :
323 : : //! Compound operator-=
324 : : //! \param[in] rhs Data object to subtract
325 : : //! \return Reference to ourselves after subtraction
326 : : Data< Layout >& operator-= ( const Data< Layout >& rhs ) {
327 : : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
328 : : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
329 : 2 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
330 : : m_vec.cbegin(), m_vec.begin(),
331 : 24 : []( tk::real s, tk::real d ){ return d-s; } );
332 : : return *this;
333 : : }
334 : : //! Operator -
335 : : //! \param[in] rhs Data object to subtract
336 : : //! \return Copy of Data object after rhs has been subtracted
337 : : //! \details Implemented in terms of compound operator-=
338 : 2 : Data< Layout > operator- ( const Data< Layout >& rhs )
339 : 2 : const { return Data< Layout >( *this ) -= rhs; }
340 : :
341 : : //! Compound operator+=
342 : : //! \param[in] rhs Data object to add
343 : : //! \return Reference to ourselves after addition
344 : : Data< Layout >& operator+= ( const Data< Layout >& rhs ) {
345 : : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
346 : : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
347 : 42551 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
348 : : m_vec.cbegin(), m_vec.begin(),
349 : 18809511 : []( tk::real s, tk::real d ){ return d+s; } );
350 : : return *this;
351 : : }
352 : : //! Operator +
353 : : //! \param[in] rhs Data object to add
354 : : //! \return Copy of Data object after rhs has been multiplied with
355 : : //! \details Implemented in terms of compound operator+=
356 : 42551 : Data< Layout > operator+ ( const Data< Layout >& rhs )
357 : 42551 : const { return Data< Layout >( *this ) += rhs; }
358 : :
359 : : //! Compound operator*= multiplying by another Data object item by item
360 : : //! \param[in] rhs Data object to multiply with
361 : : //! \return Reference to ourselves after multiplication
362 : : Data< Layout >& operator*= ( const Data< Layout >& rhs ) {
363 : : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
364 : : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
365 : 2 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
366 : : m_vec.cbegin(), m_vec.begin(),
367 : 24 : []( tk::real s, tk::real d ){ return d*s; } );
368 : : return *this;
369 : : }
370 : : //! Operator * multiplying by another Data object item by item
371 : : //! \param[in] rhs Data object to multiply with
372 : : //! \return Copy of Data object after rhs has been multiplied with
373 : : //! \details Implemented in terms of compound operator*=
374 : 2 : Data< Layout > operator* ( const Data< Layout >& rhs )
375 : 2 : const { return Data< Layout >( *this ) *= rhs; }
376 : :
377 : : //! Compound operator*= multiplying all items by a scalar
378 : : //! \param[in] rhs Scalar to multiply with
379 : : //! \return Reference to ourselves after multiplication
380 : : Data< Layout >& operator*= ( tk::real rhs ) {
381 : : // cppcheck-suppress useStlAlgorithm
382 [ + + ][ + + ]: 18852078 : for (auto& v : m_vec) v *= rhs;
[ + + ][ + + ]
[ + + ][ + + ]
383 : : return *this;
384 : : }
385 : : //! Operator * multiplying all items by a scalar
386 : : //! \param[in] rhs Scalar to multiply with
387 : : //! \return Copy of Data object after rhs has been multiplied with
388 : : //! \details Implemented in terms of compound operator*=
389 : 2 : Data< Layout > operator* ( tk::real rhs )
390 : 2 : const { return Data< Layout >( *this ) *= rhs; }
391 : :
392 : : //! Compound operator/=
393 : : //! \param[in] rhs Data object to divide by
394 : : //! \return Reference to ourselves after division
395 : : Data< Layout >& operator/= ( const Data< Layout >& rhs ) {
396 : : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
397 : : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
398 : 2 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
399 : : m_vec.cbegin(), m_vec.begin(),
400 : 24 : []( tk::real s, tk::real d ){ return d/s; } );
401 : : return *this;
402 : : }
403 : : //! Operator /
404 : : //! \param[in] rhs Data object to divide by
405 : : //! \return Copy of Data object after rhs has been divided by
406 : : //! \details Implemented in terms of compound operator/=
407 : 2 : Data< Layout > operator/ ( const Data< Layout >& rhs )
408 : 2 : const { return Data< Layout >( *this ) /= rhs; }
409 : :
410 : : //! Compound operator/= dividing all items by a scalar
411 : : //! \param[in] rhs Scalar to divide with
412 : : //! \return Reference to ourselves after division
413 : : Data< Layout >& operator/= ( tk::real rhs ) {
414 : : // cppcheck-suppress useStlAlgorithm
415 [ + + ][ + + ]: 28 : for (auto& v : m_vec) v /= rhs;
[ + + ][ + + ]
416 : : return *this;
417 : : }
418 : : //! Operator / dividing all items by a scalar
419 : : //! \param[in] rhs Scalar to divide with
420 : : //! \return Copy of Data object after rhs has been divided by
421 : : //! \details Implemented in terms of compound operator/=
422 : 2 : Data< Layout > operator/ ( tk::real rhs )
423 : 2 : const { return Data< Layout >( *this ) /= rhs; }
424 : :
425 : : //! Add new unknown at the end of the container
426 : : //! \param[in] prop Vector of properties to initialize the new unknown with
427 : : void push_back( const std::vector< tk::real >& prop )
428 [ + - ]: 158930 : { return push_back( prop, int2type< Layout >() ); }
429 : :
430 : : //! Resize data store to contain 'count' elements
431 : : //! \param[in] count Resize store to contain count * nprop elements
432 : : //! \param[in] value Value to initialize new data with (default: 0.0)
433 : : //! \note This works for both shrinking and enlarging, as this simply
434 : : //! translates to std::vector::resize(). Note that count changes, nprop
435 : : //! does not, see the private overload resize().
436 : : void resize( std::size_t count, tk::real value = 0.0 )
437 : 19073 : { resize( count, value, int2type< Layout >() ); }
438 : :
439 : : //! Remove a number of unknowns
440 : : //! \param[in] unknown Set of indices of unknowns to remove
441 : 3 : void rm( const std::set< ncomp_t >& unknown ) {
442 : : auto remove = [ &unknown ]( std::size_t i ) -> bool {
443 [ + + ]: 11 : if (unknown.find(i) != end(unknown)) return true;
444 : : return false;
445 : : };
446 : : std::size_t last = 0;
447 [ + + ]: 7 : for(std::size_t i=0; i<m_nunk; ++i, ++last) {
448 : 5 : while( remove(i) ) ++i;
449 [ + + ]: 6 : if (i >= m_nunk) break;
450 [ + + ]: 11 : for (ncomp_t p = 0; p<m_nprop; ++p)
451 : 7 : m_vec[ last*m_nprop+p ] = m_vec[ i*m_nprop+p ];
452 : : }
453 : 3 : m_vec.resize( last*m_nprop );
454 : 3 : m_nunk -= unknown.size();
455 : 3 : }
456 : :
457 : : //! Fill vector of unknowns with the same value
458 : : //! \details Requirement: component < nprop, enforced with an
459 : : //! assert in DEBUG mode, see also the constructor.
460 : : //! \param[in] component Component index, i.e., position of a scalar within
461 : : //! a system
462 : : //! \param[in] value Value to fill vector of unknowns with
463 : : inline void fill( ncomp_t component, tk::real value ) {
464 : : auto p = cptr( component );
465 [ + + ][ + + ]: 24849253 : for (ncomp_t i=0; i<m_nunk; ++i) var(p,i) = value;
[ + + ][ + + ]
[ + + ][ + + ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
466 : : }
467 : :
468 : : //! Fill full data storage with value
469 : : //! \param[in] value Value to fill data with
470 : : void fill( tk::real value )
471 : : { std::fill( begin(m_vec), end(m_vec), value ); }
472 : :
473 : : //! Check if vector of unknowns is empty
474 : : bool empty() const noexcept { return m_vec.empty(); }
475 : :
476 : : //! Layout name dispatch
477 : : //! \return The name of the data layout used
478 : : static std::string layout() { return layout( int2type< Layout >() ); }
479 : :
480 : : /** @name Pack/Unpack: Serialize Data object for Charm++ */
481 : : ///@{
482 : : //! \brief Pack/Unpack serialize member function
483 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
484 : 300793 : void pup( PUP::er &p ) {
485 : 300793 : p | m_vec;
486 : 300793 : p | m_nunk;
487 : 300793 : p | m_nprop;
488 : 300793 : }
489 : : //! \brief Pack/Unpack serialize operator|
490 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
491 : : //! \param[in,out] d DataLyaout object reference
492 : 104834 : friend void operator|( PUP::er& p, Data& d ) { d.pup(p); }
493 : : //@}
494 : :
495 : : private:
496 : : //! Transform a compile-time uint8_t into a type, used for dispatch
497 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
498 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
499 : : template< uint8_t m > struct int2type { enum { value = m }; };
500 : :
501 : : //! Overloads for the various const data accesses
502 : : //! \details Requirement: component < nprop, unknown < nunk,
503 : : //! enforced with an assert in DEBUG mode, see also the constructor.
504 : : //! \param[in] unknown Unknown index
505 : : //! \param[in] component Component index, i.e., position of a scalar within
506 : : //! a system
507 : : //! \return Const reference to data of type tk::real
508 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
509 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
510 : : const tk::real&
511 : : access( ncomp_t unknown, ncomp_t component, int2type< UnkEqComp > ) const
512 : : {
513 : : Assert( component < m_nprop, "Out-of-bounds access: "
514 : : "component < number of properties" );
515 : : Assert( unknown < m_nunk, "Out-of-bounds access: unknown < number of "
516 : : "unknowns" );
517 [ + + ][ + + ]:13116822236 : return m_vec[ unknown*m_nprop + component ];
[ + + ][ + + ]
[ + + ][ + + ]
[ + + ][ + + ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ + + ]
[ + - ][ + + ]
[ + + ][ + + ]
[ + + ][ - - ]
[ + + ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ][ + - ]
[ - - ][ - - ]
[ - - ][ + - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ + - ][ + + ]
[ + - ][ + + ]
[ + + ][ + - ]
[ + + ][ + - ]
[ + + ][ + + ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ + - ][ + + ]
[ + - ][ + + ]
[ + + ][ + - ]
[ + + ][ + - ]
[ + + ][ + + ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
518 : : }
519 : : const tk::real&
520 : : access( ncomp_t unknown, ncomp_t component, int2type< EqCompUnk > ) const
521 : : {
522 : : Assert( component < m_nprop, "Out-of-bounds access: "
523 : : "component < number of properties" );
524 : : Assert( unknown < m_nunk, "Out-of-bounds access: unknown < number of "
525 : : "unknowns" );
526 [ + - ][ - + ]: 396 : return m_vec[ (component)*m_nunk + unknown ];
[ - + ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
527 : : }
528 : :
529 : : // Overloads for the various const ptr to physical variable accesses
530 : : //! \details Requirement: component < nprop, unknown < nunk,
531 : : //! enforced with an assert in DEBUG mode, see also the constructor.
532 : : //! \param[in] component Component index, i.e., position of a scalar within
533 : : //! a system
534 : : //! \return Pointer to data of type tk::real for use with var()
535 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
536 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
537 : : const tk::real*
538 : : cptr( ncomp_t component, int2type< UnkEqComp > ) const {
539 : : Assert( component < m_nprop, "Out-of-bounds access: "
540 : : "component < number of properties" );
541 [ - - ][ - - ]: 4268194 : return m_vec.data() + component;
542 : : }
543 : : const tk::real*
544 : : cptr( ncomp_t component, int2type< EqCompUnk > ) const {
545 : : Assert( component < m_nprop, "Out-of-bounds access: "
546 : : "component < number of properties" );
547 [ + - ][ + - ]: 10 : return m_vec.data() + (component)*m_nunk;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
548 : : }
549 : :
550 : : // Overloads for the various const physical variable accesses
551 : : //! Requirement: unknown < nunk, enforced with an assert in DEBUG mode,
552 : : //! see also the constructor.
553 : : //! \param[in] pt Pointer to data of type tk::real as returned from cptr()
554 : : //! \param[in] unknown Unknown index
555 : : //! \return Const reference to data of type tk::real
556 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
557 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
558 : : inline const tk::real&
559 : : var( const tk::real* const pt, ncomp_t unknown, int2type< UnkEqComp > )
560 : : const {
561 : : Assert( unknown < m_nunk, "Out-of-bounds access: unknown < number of "
562 : : "unknowns" );
563 [ + - ][ + - ]: 1340840302 : return *(pt + unknown*m_nprop);
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
564 : : }
565 : : inline const tk::real&
566 : : var( const tk::real* const pt, ncomp_t unknown, int2type< EqCompUnk > )
567 : : const {
568 : : Assert( unknown < m_nunk, "Out-of-bounds access: unknown < number of "
569 : : "unknowns" );
570 [ + - ][ + - ]: 15 : return *(pt + unknown);
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
571 : : }
572 : :
573 : : //! Add new unknown
574 : : //! \param[in] prop Vector of properties to initialize the new unknown with
575 : : //! \note Only the UnkEqComp overload is provided as this operation would be
576 : : //! too inefficient with the EqCompUnk data layout.
577 : : void push_back( const std::vector< tk::real >& prop, int2type< UnkEqComp > )
578 : : {
579 : : Assert( prop.size() == m_nprop, "Incorrect number of properties" );
580 : : m_vec.resize( (m_nunk+1) * m_nprop );
581 : : ncomp_t u = m_nunk;
582 : : ++m_nunk;
583 : : for (ncomp_t i=0; i<m_nprop; ++i) operator()( u, i ) = prop[i];
584 : : }
585 : :
586 : : void push_back( const std::vector< tk::real >&, int2type< EqCompUnk > )
587 : : { Throw( "Not implented. It would be inefficient" ); }
588 : :
589 : : //! Resize data store to contain 'count' elements
590 : : //! \param[in] count Resize store to contain 'count' elements
591 : : //! \param[in] value Value to initialize new data with
592 : : //! \note Only the UnkEqComp overload is provided as this operation would be
593 : : //! too inefficient with the EqCompUnk data layout.
594 : : //! \note This works for both shrinking and enlarging, as this simply
595 : : //! translates to std::vector::resize().
596 : : void resize( std::size_t count, tk::real value, int2type< UnkEqComp > ) {
597 [ + - ][ + - ]: 14973 : m_vec.resize( count * m_nprop, value );
598 [ + - ][ + - ]: 12923 : m_nunk = count;
599 : : }
600 : :
601 : : void resize( std::size_t, tk::real, int2type< EqCompUnk > ) {
602 : : Throw( "Not implemented. It would be inefficient" );
603 : : }
604 : :
605 : : // Overloads for the name-queries of data lauouts
606 : : //! \return The name of the data layout used
607 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
608 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
609 : : static std::string layout( int2type< UnkEqComp > )
610 [ + - ][ + - ]: 190 : { return "unknown-major"; }
611 : : static std::string layout( int2type< EqCompUnk > )
612 [ + - ][ + - ]: 1 : { return "equation-major"; }
613 : :
614 : : std::vector< tk::real > m_vec; //!< Data pointer
615 : : ncomp_t m_nunk; //!< Number of unknowns
616 : : ncomp_t m_nprop; //!< Number of properties/unknown
617 : : };
618 : :
619 : : //! Operator * multiplying all items by a scalar from the left
620 : : //! \param[in] lhs Scalar to multiply with
621 : : //! \param[in] rhs Date object to multiply
622 : : //! \return New Data object with all items multipled with lhs
623 : : template< uint8_t Layout >
624 : 42551 : Data< Layout > operator* ( tk::real lhs, const Data< Layout >& rhs ) {
625 : 42551 : return Data< Layout >( rhs ) *= lhs;
626 : : }
627 : :
628 : : //! Operator min between two Data objects
629 : : //! \param[in] a 1st Data object
630 : : //! \param[in] b 2nd Data object
631 : : //! \return New Data object containing the minimum of all values for each
632 : : //! value in _a_ and _b_
633 : : //! \note The Data objects _a_ and _b_ must have the same number of
634 : : //! unknowns and properties.
635 : : //! \note As opposed to std::min, this function creates and returns a new object
636 : : //! instead of returning a reference to one of the operands.
637 : : template< uint8_t Layout >
638 : 2 : Data< Layout > min( const Data< Layout >& a, const Data< Layout >& b ) {
639 : : Assert( a.nunk() == b.nunk(), "Number of unknowns unequal" );
640 : : Assert( a.nprop() == b.nprop(), "Number of properties unequal" );
641 : 2 : Data< Layout > r( a.nunk(), a.nprop() );
642 : 2 : std::transform( a.vec().cbegin(), a.vec().cend(),
643 : : b.vec().cbegin(), r.vec().begin(),
644 : 12 : []( tk::real s, tk::real d ){ return std::min(s,d); } );
645 : :
646 : 2 : return r;
647 : : }
648 : :
649 : : //! Operator max between two Data objects
650 : : //! \param[in] a 1st Data object
651 : : //! \param[in] b 2nd Data object
652 : : //! \return New Data object containing the maximum of all values for each
653 : : //! value in _a_ and _b_
654 : : //! \note The Data objects _a_ and _b_ must have the same number of
655 : : //! unknowns and properties.
656 : : //! \note As opposed to std::max, this function creates and returns a new object
657 : : //! instead of returning a reference to one of the operands.
658 : : template< uint8_t Layout >
659 : 2 : Data< Layout > max( const Data< Layout >& a, const Data< Layout >& b ) {
660 : : Assert( a.nunk() == b.nunk(), "Number of unknowns unequal" );
661 : : Assert( a.nprop() == b.nprop(), "Number of properties unequal" );
662 : 2 : Data< Layout > r( a.nunk(), a.nprop() );
663 : 2 : std::transform( a.vec().cbegin(), a.vec().cend(),
664 : : b.vec().cbegin(), r.vec().begin(),
665 : 12 : []( tk::real s, tk::real d ){ return std::max(s,d); } );
666 : 2 : return r;
667 : : }
668 : :
669 : : //! Operator == between two Data objects
670 : : //! \param[in] lhs Data object to compare
671 : : //! \param[in] rhs Data object to compare
672 : : //! \return True if all entries are equal up to epsilon
673 : : template< uint8_t Layout >
674 : : bool operator== ( const Data< Layout >& lhs, const Data< Layout >& rhs ) {
675 : : Assert( rhs.nunk() == lhs.nunk(), "Incorrect number of unknowns" );
676 : : Assert( rhs.nprop() == lhs.nprop(), "Incorrect number of properties" );
677 : 4 : auto l = lhs.vec().cbegin();
678 : 4 : auto r = rhs.vec().cbegin();
679 [ + + ][ + - ]: 32 : while (l != lhs.vec().cend()) {
[ + - ][ + + ]
[ + + ][ + - ]
[ + - ][ + + ]
680 [ + - ][ - + ]: 28 : if (std::abs(*l - *r) > std::numeric_limits< tk::real >::epsilon())
[ - + ][ + - ]
[ + - ][ - + ]
[ - + ][ + - ]
681 : : return false;
682 : : ++l; ++r;
683 : : }
684 : : return true;
685 : : }
686 : :
687 : : //! Operator != between two Data objects
688 : : //! \param[in] lhs Data object to compare
689 : : //! \param[in] rhs Data object to compare
690 : : //! \return True if all entries are unequal up to epsilon
691 : : template< uint8_t Layout >
692 : : bool operator!= ( const Data< Layout >& lhs, const Data< Layout >& rhs )
693 [ + - ][ + - ]: 4 : { return !(lhs == rhs); }
[ + - ][ + - ]
694 : :
695 : : //! Compute the maximum difference between the elements of two Data objects
696 : : //! \param[in] lhs 1st Data object
697 : : //! \param[in] rhs 2nd Data object
698 : : //! \return The index, i.e., the raw position, of and the largest absolute value
699 : : //! of the difference between all corresponding elements of _lhs_ and _rhs_.
700 : : //! \details The position returned is the position in the underlying raw data
701 : : //! structure, independent of components, etc. If lhs == rhs with
702 : : //! precision std::numeric_limits< tk::real >::epsilon(), a pair of (0,0.0)
703 : : //! is returned.
704 : : //! \note The Data objects _lhs_ and _rhs_ must have the same number of
705 : : //! unknowns and properties.
706 : : template< uint8_t Layout >
707 : : std::pair< std::size_t, tk::real >
708 : 4 : maxdiff( const Data< Layout >& lhs, const Data< Layout >& rhs ) {
709 : : Assert( lhs.nunk() == rhs.nunk(), "Number of unknowns unequal" );
710 : : Assert( lhs.nprop() == rhs.nprop(), "Number of properties unequal" );
711 : 4 : auto l = lhs.vec().cbegin();
712 : 4 : auto r = rhs.vec().cbegin();
713 : 4 : std::pair< std::size_t, tk::real > m( 0, std::abs(*l - *r) );
714 : : ++l; ++r;
715 [ + + ]: 24 : while (l != lhs.vec().cend()) {
716 [ + + ]: 20 : const auto d = std::abs(*l - *r);
717 [ + + ]: 20 : if (d > m.second) m = { std::distance(lhs.vec().cbegin(),l), d };
718 : : ++l; ++r;
719 : : }
720 : 4 : return m;
721 : : }
722 : :
723 : : } // tk::
724 : :
725 : : #endif // Data_h
|
