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 546236 : static void getPointer(const TypesafePtr & p, const std::vector<std::size_t>& shape, const unsigned& start, const unsigned& stride, T*&pp ) {
30 546236 : if( shape.size()==1 && stride==1 ) {
31 40331 : pp=p.get<T*>( {shape[0]} );
32 505905 : } else if( shape.size()==1 ) {
33 388368 : auto p_=p.get<T*>( {shape[0],stride} );
34 388350 : pp = p_+start;
35 117537 : } else if( shape.size()==2 ) {
36 117537 : pp=p.get<T*>( {shape[0],shape[1]} );
37 : }
38 546218 : }
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<std::size_t>& shape, const bool& isconst ) override;
55 : /// Set the pointer to the force
56 : void setForcePointer( const TypesafePtr & p, const std::vector<std::size_t>& 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 8783 : std::unique_ptr<DataPassingObject> DataPassingObject::create(unsigned n) {
74 8783 : if(n==sizeof(double)) {
75 8783 : 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 5123 : void DataPassingObjectTyped<T>::saveValueAsDouble( const TypesafePtr & val ) {
93 5123 : hasbackup=true;
94 5123 : bvalue=double(val.template get<T>());
95 5123 : }
96 :
97 : template <class T>
98 482700 : void DataPassingObjectTyped<T>::setValuePointer( const TypesafePtr & val, const std::vector<std::size_t>& shape, const bool& isconst ) {
99 482700 : if( shape.size()==0 ) {
100 383012 : if( isconst ) {
101 382884 : val.get<const T*>();
102 : } else {
103 128 : val.get<T*>(); // just check type and discard pointer
104 : }
105 99688 : } else if( shape.size()==1 ) {
106 25655 : if( isconst )
107 25640 : val.get<const T*>({shape[0]});
108 : else
109 15 : val.get<T*>({shape[0]}); // just check type and discard pointer
110 74033 : } else if( shape.size()==2 ) {
111 74033 : if( isconst )
112 74030 : 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 482699 : v=val.copy();
117 482699 : }
118 :
119 : template <class T>
120 304899 : void DataPassingObjectTyped<T>::setForcePointer( const TypesafePtr & val, const std::vector<std::size_t>& shape ) {
121 304899 : if( shape.size()==0 ) {
122 234555 : val.get<T*>(); // just check type and discard pointer
123 70344 : } else if( shape.size()==1 )
124 0 : val.get<T*>({shape[0]}); // just check type and discard pointer
125 70344 : else if( shape.size()==2 )
126 70344 : val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
127 304883 : f=val.copy();
128 304881 : }
129 :
130 : template <class T>
131 13539 : void DataPassingObjectTyped<T>::setData( Value* value ) {
132 13539 : if( value->getRank()==0 ) {
133 703 : *v.template get<T*>() = static_cast<T>(value->get()) / unit;
134 703 : return;
135 : }
136 : T* pp;
137 12836 : getPointer( v, value->getShape(), start, stride, pp );
138 12836 : unsigned nvals=value->getNumberOfValues();
139 473006 : for(unsigned i=0; i<nvals; ++i) {
140 460170 : pp[i] = T( value->get(i) );
141 : }
142 : }
143 :
144 : template <class T>
145 264510 : void DataPassingObjectTyped<T>::share_data( const unsigned& j, const unsigned& k, Value* value ) {
146 264510 : if( value->getRank()==0 ) {
147 7704 : if( hasbackup ) {
148 7421 : value->set( unit*bvalue );
149 : } else {
150 283 : value->set( unit*double(v.template get<T>()) );
151 : }
152 7704 : return;
153 : }
154 256806 : std::vector<std::size_t> s(value->getShape());
155 256806 : if( s.size()==1 ) {
156 184709 : s[0]=k-j;
157 : }
158 : const T* pp;
159 256806 : getPointer( v, s, start, stride, pp );
160 256790 : std::vector<double> & d=value->data;
161 256790 : #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<std::size_t> 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 70880 : 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 70880 : std::vector<std::size_t> 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 70880 : maxel[0]=0;
187 : #endif
188 : const T* pp;
189 70880 : getPointer( v, maxel, start, stride, pp );
190 : // cannot be parallelized with omp because access to data is not ordered
191 : unsigned k=0;
192 3003146 : for(const auto & p : index) {
193 2932266 : value->data[p.index()]=unit*pp[i[k]*stride];
194 2932266 : k++;
195 : }
196 70880 : }
197 :
198 : template <class T>
199 45468 : void DataPassingObjectTyped<T>::add_force( Value* value ) {
200 45468 : 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 45428 : getPointer( f, value->getShape(), start, stride, pp );
206 45428 : unsigned nvals=value->getNumberOfValues();
207 45428 : #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 96131 : void DataPassingObjectTyped<T>::add_force( const std::vector<int>& index, Value* value ) {
215 96131 : plumed_assert( value->getRank()==1 );
216 96131 : std::vector<std::size_t> s(1,index.size());
217 : T* pp;
218 96131 : getPointer( f, s, start, stride, pp );
219 96129 : #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 96129 : }
224 :
225 : template <class T>
226 62655 : 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 62655 : std::vector<std::size_t> 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 62655 : maxel[0]=0;
234 : #endif
235 : T* pp;
236 62655 : getPointer( f, maxel, start, stride, pp );
237 : unsigned k=0;
238 1861383 : for(const auto & p : index) {
239 1798728 : pp[stride*i[k]] += funit*T(value->getForce(p.index()));
240 1798728 : k++;
241 : }
242 62655 : }
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<std::size_t> 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 : }
|
