Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #ifndef __PLUMED_core_ParallelTaskManager_h
23 : #define __PLUMED_core_ParallelTaskManager_h
24 :
25 : #include "ActionWithVector.h"
26 : #include "tools/Communicator.h"
27 : #include "tools/OpenMP.h"
28 : #include "tools/View.h"
29 : #include "tools/View2D.h"
30 :
31 : #include "tools/ColvarOutput.h"
32 : #include "tools/OpenACC.h"
33 :
34 : namespace PLMD {
35 :
36 : template<typename precision>
37 : struct ParActionsInput {
38 : /// Do we need to calculate the derivatives
39 : bool noderiv{false};
40 : /// Periodic boundary conditions
41 : const Pbc* pbc;
42 : /// The number of components the underlying action is computing
43 : std::size_t ncomponents{0};
44 : /// The number of scalars we are calculating for each task
45 : unsigned nscalars{0};
46 : /// The number of force scalars for each task
47 : unsigned nforcescalars{0};
48 : /// Number of derivatives for each scalar being calculated
49 : std::size_t nderivatives_per_scalar{0};
50 : /// The start of the thread unsafe forces
51 : unsigned threadunsafe_forces_start{0};
52 : /// This holds all the input data that is required to calculate all values for all tasks
53 : unsigned dataSize{0};
54 : precision *inputdata{nullptr};
55 : /// Bookeeping stuff for arguments
56 : std::size_t nargs{0};
57 : std::size_t ranks_size{0};
58 : const std::size_t* ranks{nullptr};
59 : std::size_t shapestarts_size{0};
60 : const std::size_t* shapestarts{nullptr};
61 : std::size_t shapedata_size{0};
62 : const std::size_t* shapedata{nullptr};
63 : std::size_t ncols_size{0};
64 : const std::size_t* ncols{nullptr};
65 : std::size_t bookstarts_size{0};
66 : const std::size_t* bookstarts{nullptr};
67 : std::size_t booksizes_size{0};
68 : const std::size_t* booksizes{nullptr};
69 : std::size_t bookeeping_size{0};
70 : const std::size_t* bookeeping{nullptr};
71 : std::size_t argstarts_size{0};
72 : const std::size_t* argstarts{nullptr};
73 : static ParActionsInput create( const Pbc& box ) {
74 4153 : auto toret=ParActionsInput();
75 4086 : toret.pbc=&box;
76 : return toret;
77 : }
78 :
79 : //helper function to bring data to the device in a controlled way
80 : void toACCDevice()const {
81 : #pragma acc enter data copyin(this[0:1], noderiv, pbc[0:1],ncomponents, \
82 : nscalars, nderivatives_per_scalar, threadunsafe_forces_start, \
83 : dataSize, inputdata[0:dataSize])
84 : if (nargs>0) {
85 : #pragma acc enter data copyin( nargs, \
86 : ranks_size, ranks[0:ranks_size], \
87 : shapestarts_size, shapestarts[0:shapestarts_size], \
88 : shapedata_size, shapedata[0:shapedata_size], \
89 : ncols_size, ncols[0:ncols_size], \
90 : bookstarts_size, bookstarts[0:bookstarts_size], \
91 : booksizes_size, booksizes[0:booksizes_size], \
92 : bookeeping_size, bookeeping[0:bookeeping_size], \
93 : argstarts_size, argstarts[0:argstarts_size] \
94 : )
95 : }
96 : pbc->toACCDevice();
97 : }
98 : //helper function to remove data from the device in a controlled way
99 : void removeFromACCDevice() const {
100 : pbc->removeFromACCDevice();
101 : // assuming dataSize is not changed
102 : if (nargs>0) {
103 : #pragma acc exit data delete( \
104 : shapestarts[0:shapestarts_size], shapestarts_size, \
105 : shapedata[0:shapedata_size], shapedata_size, \
106 : ncols[0:ncols_size], ncols_size, \
107 : bookstarts[0:bookstarts_size], bookstarts_size, \
108 : booksizes[0:booksizes_size], booksizes_size, \
109 : bookeeping[0:bookeeping_size], bookeeping_size, \
110 : argstarts[0:argstarts_size], argstarts_size, \
111 : ranks[0:ranks_size], ranks_size, \
112 : nargs )
113 : }
114 :
115 : #pragma acc exit data delete( \
116 : inputdata[0:dataSize], dataSize, \
117 : threadunsafe_forces_start, nderivatives_per_scalar, \
118 : nscalars, ncomponents, pbc[0:1],noderiv,this[0:1])
119 : }
120 : /// Setup the arguments
121 : void setupArguments( const ArgumentsBookkeeping& ab );
122 : unsigned sizeOfFakeVals() const {
123 : return nscalars;
124 : }
125 : };
126 :
127 : template<typename precision>
128 183432 : inline void ParActionsInput<precision>::setupArguments( const ArgumentsBookkeeping& ab ) {
129 183432 : nargs = ab.nargs;
130 183432 : ranks = ab.ranks.data();
131 183432 : ranks_size = ab.ranks.size();
132 183432 : shapestarts = ab.shapestarts.data();
133 183432 : shapestarts_size = ab.shapestarts.size();
134 183432 : shapedata = ab.shapedata.data();
135 183432 : shapedata_size = ab.shapedata.size();
136 183432 : ncols = ab.ncols.data();
137 183432 : ncols_size = ab.ncols.size();
138 183432 : bookstarts = ab.bookstarts.data();
139 183432 : bookstarts_size = ab.bookstarts.size();
140 183432 : booksizes = ab.booksizes.data();
141 183432 : booksizes_size = ab.booksizes.size();
142 183432 : bookeeping = ab.bookeeping.data();
143 183432 : bookeeping_size = ab.bookeeping.size();
144 183432 : argstarts = ab.argstarts.data();
145 183432 : argstarts_size = ab.argstarts.size();
146 183432 : }
147 :
148 : struct ArgumentBookeepingHolder {
149 : std::size_t rank;
150 : std::size_t ncols;
151 : std::size_t start;
152 : View<const std::size_t> shape;
153 : View<const std::size_t> bookeeping;
154 :
155 : template<typename precision>
156 813152 : static ArgumentBookeepingHolder create ( std::size_t argno, const ParActionsInput<precision>& inp ) {
157 : return ArgumentBookeepingHolder{
158 813152 : inp.ranks[argno], // rank
159 813152 : inp.ncols[argno], // ncols
160 813152 : inp.argstarts[argno], // start
161 813152 : View<const std::size_t>(inp.shapedata + inp.shapestarts[argno], inp.ranks[argno] ), // shape
162 813152 : View<const std::size_t>(inp.bookeeping + inp.bookstarts[argno], inp.booksizes[argno] ) // bookeeping
163 813152 : };
164 : }
165 : };
166 :
167 : template<typename precision>
168 : struct ParActionsOutput {
169 : View<precision> values;
170 : View<precision> derivatives;
171 : View<precision> buffer;
172 :
173 : static ParActionsOutput create( std::size_t ncomp,
174 : precision* v,
175 : std::size_t ndev,
176 : precision* d,
177 : std::size_t nb,
178 : precision* b ) {
179 : return ParActionsOutput{
180 : View{v,ncomp}, //values
181 : View{d,ndev}, // derivatives
182 : View{b,nb} // buffer
183 : };
184 : }
185 : };
186 :
187 : struct ForceIndexHolder {
188 : View<std::size_t> threadsafe_derivatives_end;
189 : View<std::size_t> tot_indices;
190 : View2D<std::size_t> indices;
191 :
192 : /// Constructs a ForceIndexHolder object.
193 : /// \param nc Definition here (number of components?)
194 : /// \param nd Definition here (number of derivatives?)
195 : /// \param ind Pointer to an array storing index data. It should have a size of at least 2*nc + nc*nd.
196 : static ForceIndexHolder create(const std::size_t nc,
197 : const std::size_t nd,
198 : std::size_t* ind ) {
199 : return ForceIndexHolder{
200 : View{ind,nc}, // threadsafe_derivatives_end
201 : View{ind+nc,nc}, // tot_indices
202 : View2D{ind+2*nc,nc,nd} // indices
203 : };
204 : }
205 :
206 : template<typename precision>
207 : static ForceIndexHolder create(const ParActionsInput<precision>& inp,
208 : std::size_t* ind ) {
209 : return create(inp.ncomponents,
210 : inp.nderivatives_per_scalar,ind);
211 : }
212 :
213 : /// \brief Returns the number of indexes needed by a ForceIndexHolder
214 : /// for each scalar.
215 : ///
216 : /// \param nc Definition here (number of components?)
217 : /// \param nd Definition here (number of derivatives?)
218 : ///
219 : /// \return the number of indexes needed by a ForceIndexHolder
220 : /// for each scalar.
221 : static size_t indexesPerScalar(const std::size_t nc, const std::size_t nd) {
222 : return nc // threadsafe_derivatives_end
223 : + nc // tot_indices
224 : + nc*nd; // indices
225 : }
226 : template<typename precision>
227 : static size_t indexesPerScalar(const ParActionsInput<precision>& inp) {
228 : return indexesPerScalar(inp.ncomponents,
229 : inp.nderivatives_per_scalar);
230 : }
231 : };
232 :
233 : template <typename precision>
234 : class ForcesInput {
235 : public:
236 : View<precision> force;
237 : View2D<precision> deriv;
238 : static ForcesInput create( std::size_t nv, precision* f, std::size_t nd, precision* d ) {
239 : return ForcesInput{
240 : View{f,nv}, //force
241 : View2D{d,nv,nd} //deriv
242 : };
243 : }
244 : };
245 :
246 : //There is no need to pass this as reference:
247 : template <typename precision>
248 : struct ForcesOutput {
249 : //I would suggest to invert the name or to be clearer
250 : // like something that recalls that "thread_safe" will be need to be reducted (hence it is NOT thread safe)
251 : View<precision> thread_safe;
252 : //these are the forces that we promise will not provoke races
253 : View<precision> thread_unsafe;
254 : //const T* is a ptr to const T
255 : //T* const is a conts ptr to a modifiable T
256 : static ForcesOutput create(std::vector<precision>& reduced, std::vector<precision>& notReduced) {
257 : return ForcesOutput{
258 : View{reduced.data(),reduced.size()}, // thread_safe
259 : View{notReduced.data(),notReduced.size()} // thread_unsafe
260 : };
261 : }
262 : static ForcesOutput create(precision* reduced, size_t rs, precision* notReduce, size_t nrsz) {
263 : return ForcesOutput{
264 : View{reduced,rs}, // thread_safe
265 : View{notReduce,nrsz} // thread_unsafe
266 : };
267 : }
268 : };
269 :
270 : // namespace PTMUtils {
271 : // template<class, class = void>
272 : // constexpr bool has_gatherForces_custom = false;
273 : //
274 : // //this verifies that T has a method gatherForces_custom that can be called with this signature
275 : // template<class T>
276 : // constexpr bool has_gatherForces_custom <
277 : // T,
278 : // std::void_t<
279 : // decltype(T::gatherForces_custom(
280 : // std::declval<unsigned >(),
281 : // std::declval<size_t >(),
282 : // std::declval<size_t >(),
283 : // std::declval<const typename T::input_type & >(),
284 : // std::declval<const ParActionsInput& >(),
285 : // std::declval<View<unsigned> >(),
286 : // std::declval<double *>(),
287 : // std::declval<double *>(),
288 : // std::declval<View<double> >()
289 : // ))
290 : // >
291 : // > = true;
292 : //
293 : // template<class, class = void>
294 : // constexpr bool has_gatherForces_GPU = false;
295 : //
296 : // //this verifies that T has a method gatherForces_custom that can be called with this signature
297 : // template<class T>
298 : // constexpr bool has_gatherForces_GPU <
299 : // T,
300 : // std::void_t<
301 : // decltype(T::gatherForcesGPU(
302 : // std::declval<unsigned >(),
303 : // std::declval<const typename T::input_type & >(),
304 : // std::declval<const ParActionsInput& >(),
305 : // std::declval<const ForceInput& >(),
306 : // std::declval<ForceOutput >()
307 : // ))
308 : // >
309 : // > = true;
310 : //
311 : // /// If the template class has virialSize, otherwise is 9
312 : // template<class, class=void>
313 : // constexpr size_t virialSize = 9;
314 : //
315 : // template<class T>
316 : // constexpr size_t virialSize<T, std::void_t<decltype(T::virialSize),
317 : // //this ensures that T::virialSize is a static member
318 : // std::enable_if_t<!std::is_member_pointer_v<decltype(&T::virialSize)>>
319 : // >>
320 : // = T::virialSize;
321 : // } //namespace PTMUtils
322 :
323 : template <typename CV, typename=void>
324 : struct cvprecision {
325 : typedef double type;
326 : };
327 :
328 : template <typename CV>
329 : struct cvprecision<CV,std::void_t<typename CV::precision>> {
330 : typedef typename CV::precision type;
331 : };
332 :
333 : template <typename CV>
334 : using cvprecision_t = typename cvprecision<CV>::type;
335 :
336 : template <class T>
337 : class ParallelTaskManager {
338 : public:
339 : using input_type= typename T::input_type;
340 : using precision = cvprecision_t<T>;
341 : typedef ParActionsInput<precision> ParallelActionsInput;
342 : typedef ParActionsOutput<precision> ParallelActionsOutput;
343 : typedef ForcesInput<precision> ForceInput;
344 : typedef ForcesOutput<precision> ForceOutput;
345 : protected:
346 : /// The underlying action for which we are managing parallel tasks
347 : ActionWithVector* action;
348 : /// The MPI communicator
349 : Communicator& comm;
350 : /// True if not using MPI for parllisation
351 : bool serial;
352 : /// Are we using acc for parallisation
353 : bool useacc;
354 : /// Check if the value of the derivative is known to be zero if the value
355 : /// is zero. Used to optimise force calculation
356 : bool derivativesZeroWhenValueZero;
357 : /// Number of derivatives calculated for each task
358 : std::size_t nderivatives_per_task;
359 : /// The number of forces on each thread
360 : /// The number of forces on each thread
361 : std::size_t nthreaded_forces;
362 : /// This holds the values before we pass them to the value
363 : std::vector<precision> value_stash;
364 : /// A temporary set of vectors for holding forces over threads
365 : std::vector<std::vector<precision> > omp_forces;
366 : /// This structs is used to pass data between the parallel interface and the function caller
367 : ParallelActionsInput myinput;
368 : ArgumentsBookkeeping argumentsMap;
369 : //this holds the data for myinput that will be passed though myinput
370 : std::vector<precision> input_buffer;
371 : /// This holds temporary data that we use in performTask
372 : std::size_t workspace_size;
373 : /// This holds data for that the underlying action needs to do the calculation
374 : input_type actiondata;
375 : //// This is used internally to get the number of elements in the value stash
376 : std::size_t getValueStashSize() const ;
377 : /// This is used internally to gather the forces on the threads
378 : void gatherThreads( ForceOutput forces );
379 : public:
380 : static void registerKeywords( Keywords& keys );
381 : ParallelTaskManager(ActionWithVector* av);
382 : /// Setup the parallel task manager the three arguments are
383 : /// nder = number of derivatives per scalar
384 : /// nforce_ts = number of forces that are modified by multiple tasks
385 : void setupParallelTaskManager( std::size_t nder, std::size_t nforce_ts );
386 : /// Copy the data from the underlying colvar into this parallel action
387 : void setActionInput( const input_type& adata );
388 : /// Creating the size of the workspace
389 : void setWorkspaceSize( std::size_t size );
390 : //// Set the number of force scalars per row
391 : void setNForceScalars( const unsigned& n );
392 : /// Get the action input so we can use it
393 : input_type& getActionInput();
394 : const input_type& getActionInput() const ;
395 : /// Is the calculation running in serial
396 : bool runInSerial() const {
397 9836 : return serial;
398 : }
399 : /// This runs all the tasks
400 : void runAllTasks();
401 : /// Apply the forces on the parallel object
402 : void applyForces( std::vector<double>& forcesForApply, bool reset_forces=true );
403 : /// This is used to gather forces that are thread safe
404 : static void gatherThreadSafeForces( const ParallelActionsInput& input,
405 : const ForceIndexHolder& force_indices,
406 : const ForceInput& fdata,
407 : View<precision> forces );
408 : /// This is used to gather forces that are not thread safe
409 : static void gatherThreadUnsafeForces( const ParallelActionsInput& input,
410 : const ForceIndexHolder& force_indices,
411 : const ForceInput& fdata,
412 : View<precision> forces );
413 : };
414 :
415 : struct defaultPTM {
416 : template <typename ACC>
417 : using PTM=ParallelTaskManager<ACC>;
418 : };
419 :
420 : template <class T>
421 10037 : void ParallelTaskManager<T>::registerKeywords( Keywords& keys ) {
422 10037 : keys.addFlag("SERIAL",false,"do the calculation in serial. Do not parallelize");
423 10037 : keys.addFlag("USEGPU",false,"run this calculation on the GPU");
424 20074 : keys.addLinkInDocForFlag("USEGPU","plugin_openaccPTM.md");
425 20074 : keys.addLinkInDocForFlag("SERIAL", "actions.md");
426 10037 : }
427 :
428 : template <class T>
429 3920 : ParallelTaskManager<T>::ParallelTaskManager(ActionWithVector* av):
430 3920 : action(av),
431 3920 : comm(av->comm),
432 3920 : useacc(false),
433 3920 : derivativesZeroWhenValueZero(false),
434 3920 : nderivatives_per_task(0),
435 3920 : nthreaded_forces(0),
436 : myinput(ParallelActionsInput::create(av->getPbc())),
437 3920 : workspace_size(0) {
438 3920 : action->parseFlag("USEGPU",useacc);
439 : #ifdef __PLUMED_USE_OPENACC
440 : if( useacc ) {
441 : action->log.printf(" using GPU to calculate this action\n");
442 : }
443 : #else
444 3920 : if( useacc ) {
445 0 : action->error("cannot use USEGPU flag as PLUMED has not been compiled with openacc");
446 : }
447 : #endif
448 3920 : action->parseFlag("SERIAL",serial);
449 3920 : if( serial ) {
450 0 : action->log.printf(" not using MPI to parallelise this action\n");
451 : }
452 3920 : }
453 :
454 : template <class T>
455 219329 : std::size_t ParallelTaskManager<T>::getValueStashSize() const {
456 : std::size_t valuesize=0;
457 481196 : for(unsigned i=0; i<action->getNumberOfComponents(); ++i) {
458 261867 : valuesize += action->getSizeOfValueInStash( action->getConstPntrToComponent(i) );
459 : }
460 219329 : return valuesize;
461 : }
462 :
463 :
464 : template <class T>
465 36130 : void ParallelTaskManager<T>::setupParallelTaskManager( std::size_t nder,
466 : std::size_t nforce_ts ) {
467 36130 : plumed_massert( action->getNumberOfComponents()>0, "there should be some components wen you setup the index list" );
468 36130 : myinput.ncomponents = action->getNumberOfComponents();
469 36130 : unsigned ntasks=0;
470 36130 : action->getNumberOfTasks( ntasks );
471 36130 : myinput.nscalars = 0;
472 36130 : myinput.nforcescalars = 0;
473 79585 : for(unsigned i=0; i<myinput.ncomponents; ++i) {
474 43455 : if( (action->copyOutput(i))->hasDerivatives() ) {
475 2974 : myinput.nscalars += 1 + action->getNumberOfDerivatives();
476 2974 : myinput.nforcescalars += 1;
477 40481 : } else if( (action->copyOutput(i))->getRank()==1 ) {
478 19296 : myinput.nscalars += 1;
479 19296 : myinput.nforcescalars += 1;
480 21185 : } else if( (action->copyOutput(i))->getRank()==2 ) {
481 21185 : if( ntasks==(action->copyOutput(i))->getShape()[0] ) {
482 14011 : myinput.nscalars += (action->copyOutput(i))->getNumberOfColumns();
483 14011 : myinput.nforcescalars += (action->copyOutput(i))->getNumberOfColumns();
484 : } else {
485 7174 : myinput.nscalars += 1;
486 7174 : myinput.nforcescalars += 1;
487 : }
488 : }
489 : }
490 36130 : myinput.nderivatives_per_scalar = nder;
491 36130 : nderivatives_per_task = nder*myinput.nforcescalars;
492 36130 : value_stash.resize( getValueStashSize() );
493 36130 : myinput.threadunsafe_forces_start = action->getNumberOfForceDerivatives() - nforce_ts;
494 36130 : unsigned t=OpenMP::getNumThreads();
495 36130 : if( useacc ) {
496 : t = 1;
497 : }
498 36130 : omp_forces.resize(t);
499 108390 : for(unsigned i=0; i<t; ++i) {
500 72260 : omp_forces[i].resize(nforce_ts);
501 : }
502 : // Check if we know the derivative of the component is zero if the value is zero
503 36130 : derivativesZeroWhenValueZero = true;
504 70838 : for(unsigned i=0; i<action->getNumberOfComponents(); ++i) {
505 42475 : if( !(action->getConstPntrToComponent(i))->isDerivativeZeroWhenValueIsZero() ) {
506 7767 : derivativesZeroWhenValueZero = false;
507 7767 : break;
508 : }
509 : }
510 36130 : }
511 :
512 : template <class T>
513 122 : void ParallelTaskManager<T>::setNForceScalars( const unsigned& n ) {
514 122 : myinput.nscalars = n;
515 122 : myinput.nforcescalars = n;
516 122 : nderivatives_per_task = myinput.nderivatives_per_scalar*myinput.nforcescalars;
517 122 : }
518 :
519 : template <class T>
520 1474 : void ParallelTaskManager<T>::setActionInput( const input_type& adata ) {
521 1060 : actiondata=adata;
522 1474 : }
523 :
524 : template <class T>
525 294401 : typename ParallelTaskManager<T>::input_type& ParallelTaskManager<T>::getActionInput() {
526 294401 : return actiondata;
527 : }
528 :
529 : template <class T>
530 820787 : const typename ParallelTaskManager<T>::input_type& ParallelTaskManager<T>::getActionInput() const {
531 820787 : return actiondata;
532 : }
533 :
534 : template <class T>
535 13003 : void ParallelTaskManager<T>::setWorkspaceSize( std::size_t size ) {
536 13003 : workspace_size = size;
537 13003 : }
538 :
539 : template <class T>
540 183199 : void ParallelTaskManager<T>::runAllTasks() {
541 : // Get the list of active tasks
542 183199 : std::vector<unsigned> & partialTaskList( action->getListOfActiveTasks() );
543 183199 : unsigned nactive_tasks=partialTaskList.size();
544 : // Get all the input data so we can broadcast it to the GPU
545 183199 : myinput.noderiv = true;
546 394476 : for(unsigned i=0; i<action->getNumberOfComponents(); ++i) {
547 218412 : if( (action->getConstPntrToComponent(i))->hasDerivatives() ) {
548 7135 : myinput.noderiv=false;
549 7135 : break;
550 : }
551 : }
552 183199 : action->getInputData( input_buffer, argumentsMap );
553 183199 : myinput.dataSize = input_buffer.size();
554 183199 : myinput.inputdata = input_buffer.data();
555 : // Transfer all the bookeeping information about the arguments
556 183199 : myinput.setupArguments( argumentsMap );
557 : // Reset the values at the start of the task loop
558 183199 : std::size_t totalvals=getValueStashSize();
559 183199 : if( value_stash.size()!=totalvals ) {
560 62 : value_stash.resize(totalvals);
561 : }
562 : std::fill (value_stash.begin(),value_stash.end(), precision(0.0));
563 183199 : if( useacc ) {
564 : // there should not be any available path to get here
565 0 : plumed_merror("to use the USEGPU flag you need to LOAD a plugin like \"openaccPTM\"");
566 : } else {
567 : // Get the MPI details
568 183199 : unsigned stride=comm.Get_size();
569 183199 : unsigned rank=comm.Get_rank();
570 183199 : if(serial) {
571 : stride=1;
572 : rank=0;
573 : }
574 :
575 : // Get number of threads for OpenMP
576 183199 : unsigned nt=OpenMP::getNumThreads();
577 183199 : if( nt*stride*10>nactive_tasks ) {
578 101146 : nt=nactive_tasks/stride/10;
579 : }
580 : if( nt==0 ) {
581 : nt=1;
582 : }
583 :
584 183199 : #pragma omp parallel num_threads(nt)
585 : {
586 : std::vector<precision> buffer( workspace_size );
587 : std::vector<precision> derivatives( nderivatives_per_task );
588 : #pragma omp for nowait
589 : for(unsigned i=rank; i<nactive_tasks; i+=stride) {
590 : std::size_t task_index = partialTaskList[i];
591 : std::size_t val_pos = task_index*myinput.nscalars;
592 : auto myout = ParallelActionsOutput::create ( myinput.nscalars,
593 : value_stash.data()+val_pos,
594 : nderivatives_per_task,
595 : derivatives.data(),
596 : workspace_size,
597 : buffer.data() );
598 : // Calculate the stuff in the loop for this action
599 : T::performTask( task_index, actiondata, myinput, myout );
600 : }
601 : }
602 : // MPI Gather everything
603 183199 : if( !serial ) {
604 183199 : comm.Sum( value_stash );
605 : }
606 : }
607 : // And transfer the value to the output values
608 183199 : action->transferStashToValues( partialTaskList, value_stash );
609 183199 : }
610 :
611 : template<typename prec>
612 : struct forceData {
613 : std::vector<prec> ffa;
614 : // the forces already in forcesForApply are kept: the caller may be
615 : // accumulating over several applyForces calls (reset_forces=false)
616 : forceData( std::vector<double>& forcesForApply ):
617 : ffa(forcesForApply.begin(),forcesForApply.end()),destination(forcesForApply) {}
618 : void update() {
619 : std::copy(ffa.begin(),ffa.end(),destination.begin());
620 : }
621 : private:
622 : std::vector<double>& destination;
623 : };
624 :
625 : template<>
626 : struct forceData<double> {
627 : PLMD::View<double> ffa;
628 : forceData( std::vector<double>& forcesForApply ):
629 : ffa(forcesForApply.data(),forcesForApply.size()) {}
630 : void update() {
631 : }
632 : };
633 : template <class T>
634 142436 : void ParallelTaskManager<T>::applyForces( std::vector<double>& forcesForApply, bool reset_forces ) {
635 : // Get the list of active tasks
636 142436 : std::vector<unsigned> & partialTaskList= action->getListOfActiveTasks() ;
637 142436 : unsigned nactive_tasks=partialTaskList.size();
638 : forceData<precision> forces(forcesForApply);
639 : // Only set to false if we are doing a loop over the columns of a matrix to avoid
640 : // having to gather forces across multiple nodes
641 142436 : if( reset_forces ) {
642 : // Clear force buffer
643 : std::fill (forces.ffa.begin(),forces.ffa.end(), precision(0.0));
644 : }
645 : // Get all the input data so we can broadcast it to the GPU
646 142436 : myinput.noderiv = false;
647 : // Retrieve the forces from the values
648 142436 : action->transferForcesToStash( partialTaskList, value_stash );
649 :
650 142436 : if( useacc ) {
651 : // there should not be any available path to get here
652 0 : plumed_merror("to use the USEGPU flag you need to LOAD a plugin like \"openaccPTM\"");
653 : } else {
654 : // Get the MPI details
655 142436 : unsigned stride=comm.Get_size();
656 142436 : unsigned rank=comm.Get_rank();
657 142436 : if(serial) {
658 : stride=1;
659 : rank=0;
660 : }
661 :
662 : // Get number of threads for OpenMP
663 142436 : unsigned nt=OpenMP::getNumThreads();
664 142436 : if( nt*stride*10>nactive_tasks ) {
665 77520 : nt=nactive_tasks/stride/10;
666 : }
667 : if( nt==0 ) {
668 : nt=1;
669 : }
670 142436 : #pragma omp parallel num_threads(nt)
671 : {
672 : const unsigned t=OpenMP::getThreadNum();
673 : omp_forces[t].assign( omp_forces[t].size(), 0.0 );
674 : std::vector<precision> buffer( workspace_size );
675 : std::vector<precision> fake_vals( myinput.sizeOfFakeVals() );
676 : std::vector<precision> derivatives( nderivatives_per_task );
677 : std::vector<std::size_t> indices(ForceIndexHolder::indexesPerScalar(myinput));
678 :
679 : auto force_indices = ForceIndexHolder::create( myinput,indices.data() );
680 : #pragma omp for nowait
681 : for(unsigned i=rank; i<nactive_tasks; i+=stride) {
682 : std::size_t task_index = partialTaskList[i];
683 : auto myout = ParallelActionsOutput::create( myinput.nscalars,
684 : fake_vals.data(),
685 : derivatives.size(),
686 : derivatives.data(),
687 : workspace_size,
688 : buffer.data() );
689 : // Calculate the stuff in the loop for this action
690 : T::performTask( task_index, actiondata, myinput, myout );
691 :
692 : // If this is a matrix this returns a number that isn't one as we have to loop over the columns
693 : const std::size_t nvpt = T::getNumberOfValuesPerTask( task_index, actiondata );
694 : for(unsigned j=0; j<nvpt; ++j) {
695 : bool noforces = false;
696 : // Skip the force calculation if we have been told that we can do this when the value is zero
697 : if( derivativesZeroWhenValueZero ) {
698 : noforces = true;
699 : for(unsigned k=0; k<myinput.ncomponents; ++k) {
700 : if( fabs(fake_vals[j*myinput.ncomponents+k])>epsilon ) {
701 : noforces = false;
702 : }
703 : }
704 : }
705 : if( !noforces ) {
706 : // Get the force indices
707 : T::getForceIndices( task_index,
708 : j,
709 : forces.ffa.size(),
710 : actiondata,
711 : myinput,
712 : force_indices );
713 : // Create a force input object
714 : auto finput=ForceInput::create( myinput.nforcescalars,
715 : value_stash.data()
716 : + myinput.nforcescalars*task_index
717 : + j*myinput.ncomponents,
718 : myinput.nderivatives_per_scalar,
719 : derivatives.data()
720 : + j*myinput.ncomponents*myinput.nderivatives_per_scalar );
721 :
722 : // Gather forces that are thread safe
723 : gatherThreadSafeForces( myinput,
724 : force_indices,
725 : finput,
726 : View<double>(forces.ffa.data(),
727 : forces.ffa.size()) );
728 :
729 : // Gather forces that are not thread safe
730 : gatherThreadUnsafeForces( myinput,
731 : force_indices,
732 : finput,
733 : View<precision>(omp_forces[t].data(),
734 : omp_forces[t].size()) );
735 : }
736 : }
737 : }
738 :
739 : #pragma omp critical
740 : gatherThreads( ForceOutput::create(omp_forces[t], forcesForApply ) );
741 : }
742 : forces.update();
743 : // MPI Gather everything (this must be extended to the gpu thing, after makning it mpi-aware)
744 142436 : if( !serial ) {
745 142436 : comm.Sum( forcesForApply );
746 : }
747 : }
748 142436 : }
749 :
750 : template <class T>
751 16827294 : void ParallelTaskManager<T>::gatherThreadSafeForces( const ParallelActionsInput& input,
752 : const ForceIndexHolder& force_indices,
753 : const ForceInput& fdata,
754 : View<precision> forces ) {
755 50115034 : for(unsigned i=0; i<input.ncomponents; ++i) {
756 33287740 : double ff = fdata.force[i];
757 117423093 : for(unsigned j=0; j<force_indices.threadsafe_derivatives_end[i]; ++j) {
758 84135353 : forces[ force_indices.indices[i][j] ] += ff*fdata.deriv[i][j];
759 : }
760 : }
761 16827294 : }
762 :
763 : template <class T>
764 16827294 : void ParallelTaskManager<T>::gatherThreadUnsafeForces(const ParallelActionsInput& input,
765 : const ForceIndexHolder& force_indices,
766 : const ForceInput& fdata,
767 : View<precision> forces ) {
768 50115034 : for(unsigned i=0; i<input.ncomponents; ++i) {
769 33287740 : const double ff = fdata.force[i];
770 195692138 : for(unsigned d=force_indices.threadsafe_derivatives_end[i];
771 195692138 : d<force_indices.tot_indices[i]; ++d) {
772 162404398 : forces[ force_indices.indices[i][d] - input.threadunsafe_forces_start ]
773 162404398 : += ff*fdata.deriv[i][d];
774 : }
775 : }
776 16827294 : }
777 :
778 : template <class T>
779 207352 : void ParallelTaskManager<T>::gatherThreads( ForceOutput forces ) {
780 : //Forceoutput is basically two spans, so it is ok to pass it by value
781 : unsigned k=0;
782 8692140 : for(unsigned n=forces.thread_unsafe.size()-forces.thread_safe.size(); n<forces.thread_unsafe.size(); ++n) {
783 8484788 : forces.thread_unsafe[n] += forces.thread_safe[k];
784 8484788 : ++k;
785 : }
786 207352 : }
787 :
788 : } // namespace PLMD
789 : #endif
|
