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 515290 : static void getPointer(const TypesafePtr & p, const std::vector<unsigned>& shape, const unsigned& start, const unsigned& stride, T*&pp ) {
30 515290 : if( shape.size()==1 && stride==1 ) {
31 37226 : pp=p.get<T*>( {shape[0]} );
32 478064 : } else if( shape.size()==1 ) {
33 367488 : auto p_=p.get<T*>( {shape[0],stride} );
34 367470 : pp = p_+start;
35 110576 : } else if( shape.size()==2 ) {
36 110576 : pp=p.get<T*>( {shape[0],shape[1]} );
37 : }
38 515272 : }
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 8634 : std::unique_ptr<DataPassingObject> DataPassingObject::create(unsigned n) {
74 8634 : if(n==sizeof(double)) {
75 8627 : return std::make_unique<DataPassingObjectTyped<double>>();
76 7 : } else if(n==sizeof(float)) {
77 7 : 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 5102 : void DataPassingObjectTyped<T>::saveValueAsDouble( const TypesafePtr & val ) {
93 5102 : hasbackup=true;
94 5102 : bvalue=double(val.template get<T>());
95 5102 : }
96 :
97 : template <>
98 1 : void DataPassingObjectTyped<float>::saveValueAsDouble( const TypesafePtr & val ) {
99 1 : hasbackup=true;
100 1 : bvalue=double(val.template get<float>());
101 : // The following is to avoid extra digits in case the MD code uses floats
102 : // e.g.: float f=0.002 when converted to double becomes 0.002000000094995
103 : // To avoid this, we keep only up to 6 significant digits after first one
104 1 : double magnitude=std::pow(10,std::floor(std::log10(bvalue)));
105 1 : bvalue = std::round(bvalue/magnitude*1e6)/1e6*magnitude;
106 1 : }
107 :
108 : template <class T>
109 459939 : void DataPassingObjectTyped<T>::setValuePointer( const TypesafePtr & val, const std::vector<unsigned>& shape, const bool& isconst ) {
110 459939 : if( shape.size()==0 ) {
111 365865 : if( isconst ) {
112 365737 : val.get<const T*>();
113 : } else {
114 128 : val.get<T*>(); // just check type and discard pointer
115 : }
116 94074 : } else if( shape.size()==1 ) {
117 23470 : if( isconst )
118 23456 : val.get<const T*>({shape[0]});
119 : else
120 14 : val.get<T*>({shape[0]}); // just check type and discard pointer
121 70604 : } else if( shape.size()==2 ) {
122 70604 : if( isconst )
123 70601 : val.get<const T*>({shape[0],shape[1]});
124 : else
125 3 : val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
126 : }
127 459938 : v=val.copy();
128 459938 : }
129 :
130 : template <class T>
131 291183 : void DataPassingObjectTyped<T>::setForcePointer( const TypesafePtr & val, const std::vector<unsigned>& shape ) {
132 291183 : if( shape.size()==0 ) {
133 224268 : val.get<T*>(); // just check type and discard pointer
134 66915 : } else if( shape.size()==1 )
135 0 : val.get<T*>({shape[0]}); // just check type and discard pointer
136 66915 : else if( shape.size()==2 )
137 66915 : val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
138 291167 : f=val.copy();
139 291165 : }
140 :
141 : template <class T>
142 12447 : void DataPassingObjectTyped<T>::setData( Value* value ) {
143 12447 : if( value->getRank()==0 ) {
144 703 : *v.template get<T*>() = static_cast<T>(value->get()) / unit;
145 703 : return;
146 : }
147 : T* pp;
148 11744 : getPointer( v, value->getShape(), start, stride, pp );
149 11744 : unsigned nvals=value->getNumberOfValues();
150 429326 : for(unsigned i=0; i<nvals; ++i) {
151 417582 : pp[i] = T( value->get(i) );
152 : }
153 : }
154 :
155 : template <class T>
156 248449 : void DataPassingObjectTyped<T>::share_data( const unsigned& j, const unsigned& k, Value* value ) {
157 248449 : if( value->getRank()==0 ) {
158 7433 : if( hasbackup ) {
159 7361 : value->set( unit*bvalue );
160 : } else {
161 72 : value->set( unit*double(v.template get<T>()) );
162 : }
163 7433 : return;
164 : }
165 241016 : std::vector<unsigned> s(value->getShape());
166 241016 : if( s.size()==1 ) {
167 172409 : s[0]=k-j;
168 : }
169 : const T* pp;
170 241016 : getPointer( v, s, start, stride, pp );
171 241000 : std::vector<double> & d=value->data;
172 241000 : #pragma omp parallel for num_threads(value->getGoodNumThreads(j,k))
173 : for(unsigned i=j; i<k; ++i) {
174 : d[i]=unit*pp[i*stride];
175 : }
176 : }
177 :
178 : template <class T>
179 1350 : void DataPassingObjectTyped<T>::share_data( std::vector<double>& values ) const {
180 1350 : std::vector<unsigned> maxel(1,values.size());
181 : const T* pp;
182 1350 : getPointer( v, maxel, start, stride, pp );
183 1350 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(values))
184 : for(unsigned i=0; i<values.size(); ++i) {
185 : values[i]=unit*pp[i*stride];
186 : }
187 1350 : }
188 :
189 : template <class T>
190 70700 : void DataPassingObjectTyped<T>::share_data( const std::vector<AtomNumber>&index, const std::vector<unsigned>& i, Value* value ) {
191 : plumed_dbg_assert( value->getRank()==1 );
192 70700 : std::vector<unsigned> maxel(1,index.size());
193 : #ifndef NDEBUG
194 : // bounds are only checked in debug mode since they require this extra step that is potentially expensive
195 : maxel[0]=(i.size()>0?*std::max_element(i.begin(),i.end())+1:0);
196 : #else
197 70700 : maxel[0]=0;
198 : #endif
199 : const T* pp;
200 70700 : getPointer( v, maxel, start, stride, pp );
201 : // cannot be parallelized with omp because access to data is not ordered
202 : unsigned k=0;
203 2507486 : for(const auto & p : index) {
204 2436786 : value->data[p.index()]=unit*pp[i[k]*stride];
205 2436786 : k++;
206 : }
207 70700 : }
208 :
209 : template <class T>
210 41997 : void DataPassingObjectTyped<T>::add_force( Value* value ) {
211 41997 : if( value->getRank()==0 ) {
212 40 : *f.template get<T*>() += funit*static_cast<T>(value->getForce(0));
213 40 : return;
214 : }
215 : T* pp;
216 41957 : getPointer( f, value->getShape(), start, stride, pp );
217 41957 : unsigned nvals=value->getNumberOfValues();
218 41957 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,nvals))
219 : for(unsigned i=0; i<nvals; ++i) {
220 : pp[i*stride] += funit*T(value->getForce(i));
221 : }
222 : }
223 :
224 : template <class T>
225 85784 : void DataPassingObjectTyped<T>::add_force( const std::vector<int>& index, Value* value ) {
226 85784 : plumed_assert( value->getRank()==1 );
227 85784 : std::vector<unsigned> s(1,index.size());
228 : T* pp;
229 85784 : getPointer( f, s, start, stride, pp );
230 85782 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,index.size()))
231 : for(unsigned i=0; i<index.size(); ++i) {
232 : pp[i*stride] += funit*T(value->getForce(index[i]));
233 : }
234 85782 : }
235 :
236 : template <class T>
237 62589 : void DataPassingObjectTyped<T>::add_force( const std::vector<AtomNumber>& index, const std::vector<unsigned>& i, Value* value ) {
238 : plumed_dbg_assert( value->getRank()==1 );
239 62589 : std::vector<unsigned> maxel(1,index.size());
240 : #ifndef NDEBUG
241 : // bounds are only checked in debug mode since they require this extra step that is potentially expensive
242 : maxel[0]=(i.size()>0?*std::max_element(i.begin(),i.end())+1:0);
243 : #else
244 62589 : maxel[0]=0;
245 : #endif
246 : T* pp;
247 62589 : getPointer( f, maxel, start, stride, pp );
248 : unsigned k=0;
249 1856013 : for(const auto & p : index) {
250 1793424 : pp[stride*i[k]] += funit*T(value->getForce(p.index()));
251 1793424 : k++;
252 : }
253 62589 : }
254 :
255 : template <class T>
256 150 : void DataPassingObjectTyped<T>::rescale_force( const unsigned& n, const double& factor, Value* value ) {
257 150 : plumed_assert( value->getRank()>0 );
258 150 : std::vector<unsigned> s( value->getShape() );
259 150 : if( s.size()==1 ) {
260 150 : s[0] = n;
261 : }
262 : T* pp;
263 150 : getPointer( f, s, start, stride, pp );
264 150 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,n))
265 : for(unsigned i=0; i<n; ++i) {
266 : pp[i*stride] *= T(factor);
267 : }
268 150 : }
269 :
270 : }
|
