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 : #include "DataPassingObject.h"
23 : #include "tools/OpenMP.h"
24 : #include "tools/Tools.h"
25 :
26 : namespace PLMD {
27 :
28 : template<typename T>
29 515212 : static void getPointer(const TypesafePtr & p, const std::vector<unsigned>& shape, const unsigned& start, const unsigned& stride, T*&pp ) {
30 515212 : if( shape.size()==1 && stride==1 ) {
31 37214 : pp=p.get<T*>( {shape[0]} );
32 477998 : } else if( shape.size()==1 ) {
33 367437 : auto p_=p.get<T*>( {shape[0],stride} );
34 367419 : pp = p_+start;
35 110561 : } else if( shape.size()==2 ) {
36 110561 : pp=p.get<T*>( {shape[0],shape[1]} );
37 : }
38 515194 : }
39 :
40 : template <class T>
41 : class DataPassingObjectTyped : public DataPassingObject {
42 : private:
43 : /// A pointer to the value
44 : TypesafePtr v;
45 : /// A pointer to the force
46 : TypesafePtr f;
47 : public:
48 : /// This convers a number from the MD code into a double
49 : double MD2double(const TypesafePtr &) const override ;
50 : /// This is used when you want to save the passed object to a double variable in PLUMED rather than the pointer
51 : /// this can be used even when you don't pass a pointer from the MD code
52 : void saveValueAsDouble( const TypesafePtr & val ) override;
53 : /// Set the pointer to the value
54 : void setValuePointer( const TypesafePtr & p, const std::vector<unsigned>& shape, const bool& isconst ) override;
55 : /// Set the pointer to the force
56 : void setForcePointer( const TypesafePtr & p, const std::vector<unsigned>& shape ) override;
57 : /// This gets the data in the pointer and passes it to the output value
58 : void share_data( std::vector<double>& values ) const override ;
59 : /// Share the data and put it in the value from sequential data
60 : void share_data( const unsigned& j, const unsigned& k, Value* value ) override;
61 : /// Share the data and put it in the value from a scattered data
62 : void share_data( const std::vector<AtomNumber>&index, const std::vector<unsigned>& i, Value* value ) override;
63 : /// Pass the force from the value to the output value
64 : void add_force( Value* vv ) override;
65 : void add_force( const std::vector<int>& index, Value* value ) override;
66 : void add_force( const std::vector<AtomNumber>& index, const std::vector<unsigned>& i, Value* value ) override;
67 : /// Rescale the force on the output value
68 : void rescale_force( const unsigned& n, const double& factor, Value* value ) override;
69 : /// This transfers everything to the output
70 : void setData( Value* value ) override;
71 : };
72 :
73 8580 : std::unique_ptr<DataPassingObject> DataPassingObject::create(unsigned n) {
74 8580 : if(n==sizeof(double)) {
75 8580 : return std::make_unique<DataPassingObjectTyped<double>>();
76 0 : } else if(n==sizeof(float)) {
77 0 : return std::make_unique<DataPassingObjectTyped<float>>();
78 : }
79 : std::string pp;
80 0 : Tools::convert(n,pp);
81 0 : plumed_merror("cannot create an MD interface with sizeof(real)=="+ pp);
82 : return NULL;
83 : }
84 :
85 : template <class T>
86 0 : double DataPassingObjectTyped<T>::MD2double(const TypesafePtr & m) const {
87 0 : double d=double(m.template get<T>());
88 0 : return d;
89 : }
90 :
91 : template <class T>
92 5097 : void DataPassingObjectTyped<T>::saveValueAsDouble( const TypesafePtr & val ) {
93 5097 : hasbackup=true;
94 5097 : bvalue=double(val.template get<T>());
95 5097 : }
96 :
97 : template <class T>
98 459865 : void DataPassingObjectTyped<T>::setValuePointer( const TypesafePtr & val, const std::vector<unsigned>& shape, const bool& isconst ) {
99 459865 : if( shape.size()==0 ) {
100 365804 : if( isconst ) {
101 365676 : val.get<const T*>();
102 : } else {
103 128 : val.get<T*>(); // just check type and discard pointer
104 : }
105 94061 : } else if( shape.size()==1 ) {
106 23470 : if( isconst )
107 23456 : val.get<const T*>({shape[0]});
108 : else
109 14 : val.get<T*>({shape[0]}); // just check type and discard pointer
110 70591 : } else if( shape.size()==2 ) {
111 70591 : if( isconst )
112 70588 : val.get<const T*>({shape[0],shape[1]});
113 : else
114 3 : val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
115 : }
116 459864 : v=val.copy();
117 459864 : }
118 :
119 : template <class T>
120 291131 : void DataPassingObjectTyped<T>::setForcePointer( const TypesafePtr & val, const std::vector<unsigned>& shape ) {
121 291131 : if( shape.size()==0 ) {
122 224229 : val.get<T*>(); // just check type and discard pointer
123 66902 : } else if( shape.size()==1 )
124 0 : val.get<T*>({shape[0]}); // just check type and discard pointer
125 66902 : else if( shape.size()==2 )
126 66902 : val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
127 291115 : f=val.copy();
128 291113 : }
129 :
130 : template <class T>
131 12447 : void DataPassingObjectTyped<T>::setData( Value* value ) {
132 12447 : if( value->getRank()==0 ) {
133 703 : *v.template get<T*>() = static_cast<T>(value->get()) / unit;
134 703 : return;
135 : }
136 : T* pp;
137 11744 : getPointer( v, value->getShape(), start, stride, pp );
138 11744 : unsigned nvals=value->getNumberOfValues();
139 429326 : for(unsigned i=0; i<nvals; ++i) {
140 417582 : pp[i] = T( value->get(i) );
141 : }
142 : }
143 :
144 : template <class T>
145 248389 : void DataPassingObjectTyped<T>::share_data( const unsigned& j, const unsigned& k, Value* value ) {
146 248389 : if( value->getRank()==0 ) {
147 7415 : if( hasbackup ) {
148 7343 : value->set( unit*bvalue );
149 : } else {
150 72 : value->set( unit*double(v.template get<T>()) );
151 : }
152 7415 : return;
153 : }
154 240974 : std::vector<unsigned> s(value->getShape());
155 240974 : if( s.size()==1 ) {
156 172378 : s[0]=k-j;
157 : }
158 : const T* pp;
159 240974 : getPointer( v, s, start, stride, pp );
160 240958 : std::vector<double> & d=value->data;
161 240958 : #pragma omp parallel for num_threads(value->getGoodNumThreads(j,k))
162 : for(unsigned i=j; i<k; ++i) {
163 : d[i]=unit*pp[i*stride];
164 : }
165 : }
166 :
167 : template <class T>
168 1350 : void DataPassingObjectTyped<T>::share_data( std::vector<double>& values ) const {
169 1350 : std::vector<unsigned> maxel(1,values.size());
170 : const T* pp;
171 1350 : getPointer( v, maxel, start, stride, pp );
172 1350 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(values))
173 : for(unsigned i=0; i<values.size(); ++i) {
174 : values[i]=unit*pp[i*stride];
175 : }
176 1350 : }
177 :
178 : template <class T>
179 70680 : void DataPassingObjectTyped<T>::share_data( const std::vector<AtomNumber>&index, const std::vector<unsigned>& i, Value* value ) {
180 : plumed_dbg_assert( value->getRank()==1 );
181 70680 : std::vector<unsigned> maxel(1,index.size());
182 : #ifndef NDEBUG
183 : // bounds are only checked in debug mode since they require this extra step that is potentially expensive
184 : maxel[0]=(i.size()>0?*std::max_element(i.begin(),i.end())+1:0);
185 : #else
186 70680 : maxel[0]=0;
187 : #endif
188 : const T* pp;
189 70680 : getPointer( v, maxel, start, stride, pp );
190 : // cannot be parallelized with omp because access to data is not ordered
191 : unsigned k=0;
192 2507350 : for(const auto & p : index) {
193 2436670 : value->data[p.index()]=unit*pp[i[k]*stride];
194 2436670 : k++;
195 : }
196 70680 : }
197 :
198 : template <class T>
199 41993 : void DataPassingObjectTyped<T>::add_force( Value* value ) {
200 41993 : if( value->getRank()==0 ) {
201 40 : *f.template get<T*>() += funit*static_cast<T>(value->getForce(0));
202 40 : return;
203 : }
204 : T* pp;
205 41953 : getPointer( f, value->getShape(), start, stride, pp );
206 41953 : unsigned nvals=value->getNumberOfValues();
207 41953 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,nvals))
208 : for(unsigned i=0; i<nvals; ++i) {
209 : pp[i*stride] += funit*T(value->getForce(i));
210 : }
211 : }
212 :
213 : template <class T>
214 85772 : void DataPassingObjectTyped<T>::add_force( const std::vector<int>& index, Value* value ) {
215 85772 : plumed_assert( value->getRank()==1 );
216 85772 : std::vector<unsigned> s(1,index.size());
217 : T* pp;
218 85772 : getPointer( f, s, start, stride, pp );
219 85770 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,index.size()))
220 : for(unsigned i=0; i<index.size(); ++i) {
221 : pp[i*stride] += funit*T(value->getForce(index[i]));
222 : }
223 85770 : }
224 :
225 : template <class T>
226 62589 : void DataPassingObjectTyped<T>::add_force( const std::vector<AtomNumber>& index, const std::vector<unsigned>& i, Value* value ) {
227 : plumed_dbg_assert( value->getRank()==1 );
228 62589 : std::vector<unsigned> maxel(1,index.size());
229 : #ifndef NDEBUG
230 : // bounds are only checked in debug mode since they require this extra step that is potentially expensive
231 : maxel[0]=(i.size()>0?*std::max_element(i.begin(),i.end())+1:0);
232 : #else
233 62589 : maxel[0]=0;
234 : #endif
235 : T* pp;
236 62589 : getPointer( f, maxel, start, stride, pp );
237 : unsigned k=0;
238 1856013 : for(const auto & p : index) {
239 1793424 : pp[stride*i[k]] += funit*T(value->getForce(p.index()));
240 1793424 : k++;
241 : }
242 62589 : }
243 :
244 : template <class T>
245 150 : void DataPassingObjectTyped<T>::rescale_force( const unsigned& n, const double& factor, Value* value ) {
246 150 : plumed_assert( value->getRank()>0 );
247 150 : std::vector<unsigned> s( value->getShape() );
248 150 : if( s.size()==1 ) {
249 150 : s[0] = n;
250 : }
251 : T* pp;
252 150 : getPointer( f, s, start, stride, pp );
253 150 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,n))
254 : for(unsigned i=0; i<n; ++i) {
255 : pp[i*stride] *= T(factor);
256 : }
257 150 : }
258 :
259 : }
|
