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 498390 : static void getPointer(const TypesafePtr & p, const std::vector<unsigned>& shape, const unsigned& start, const unsigned& stride, T*&pp ) {
30 498390 : if( shape.size()==1 && stride==1 ) { pp=p.get<T*>( {shape[0]} ); }
31 461218 : else if( shape.size()==1 ) { auto p_=p.get<T*>( {shape[0],stride} ); pp = p_+start; }
32 106366 : else if( shape.size()==2 ) { pp=p.get<T*>( {shape[0],shape[1]} ); }
33 498372 : }
34 :
35 : template <class T>
36 : class DataPassingObjectTyped : public DataPassingObject {
37 : private:
38 : /// A pointer to the value
39 : TypesafePtr v;
40 : /// A pointer to the force
41 : TypesafePtr f;
42 : public:
43 : /// This convers a number from the MD code into a double
44 : double MD2double(const TypesafePtr &) const override ;
45 : /// This is used when you want to save the passed object to a double variable in PLUMED rather than the pointer
46 : /// this can be used even when you don't pass a pointer from the MD code
47 : void saveValueAsDouble( const TypesafePtr & val ) override;
48 : /// Set the pointer to the value
49 : void setValuePointer( const TypesafePtr & p, const std::vector<unsigned>& shape, const bool& isconst ) override;
50 : /// Set the pointer to the force
51 : void setForcePointer( const TypesafePtr & p, const std::vector<unsigned>& shape ) override;
52 : /// This gets the data in the pointer and passes it to the output value
53 : void share_data( std::vector<double>& values ) const override ;
54 : /// Share the data and put it in the value from sequential data
55 : void share_data( const unsigned& j, const unsigned& k, Value* value ) override;
56 : /// Share the data and put it in the value from a scattered data
57 : void share_data( const std::vector<AtomNumber>&index, const std::vector<unsigned>& i, Value* value ) override;
58 : /// Pass the force from the value to the output value
59 : void add_force( Value* vv ) override;
60 : void add_force( const std::vector<int>& index, Value* value ) override;
61 : void add_force( const std::vector<AtomNumber>& index, const std::vector<unsigned>& i, Value* value ) override;
62 : /// Rescale the force on the output value
63 : void rescale_force( const unsigned& n, const double& factor, Value* value ) override;
64 : /// This transfers everything to the output
65 : void setData( Value* value ) override;
66 : };
67 :
68 8433 : std::unique_ptr<DataPassingObject> DataPassingObject::create(unsigned n) {
69 8433 : if(n==sizeof(double)) {
70 8433 : return std::make_unique<DataPassingObjectTyped<double>>();
71 0 : } else if(n==sizeof(float)) {
72 0 : return std::make_unique<DataPassingObjectTyped<float>>();
73 : }
74 0 : std::string pp; Tools::convert(n,pp);
75 0 : plumed_merror("cannot create an MD interface with sizeof(real)=="+ pp);
76 : return NULL;
77 : }
78 :
79 : template <class T>
80 0 : double DataPassingObjectTyped<T>::MD2double(const TypesafePtr & m) const {
81 0 : double d=double(m.template get<T>()); return d;
82 : }
83 :
84 : template <class T>
85 5076 : void DataPassingObjectTyped<T>::saveValueAsDouble( const TypesafePtr & val ) {
86 5076 : hasbackup=true; bvalue=double(val.template get<T>());
87 5076 : }
88 :
89 : template <class T>
90 447223 : void DataPassingObjectTyped<T>::setValuePointer( const TypesafePtr & val, const std::vector<unsigned>& shape, const bool& isconst ) {
91 447223 : if( shape.size()==0 ) {
92 355269 : if( isconst ) val.get<const T*>(); else val.get<T*>(); // just check type and discard pointer
93 91954 : } else if( shape.size()==1 ) {
94 23470 : if( isconst ) val.get<const T*>({shape[0]}); else val.get<T*>({shape[0]}); // just check type and discard pointer
95 68484 : } else if( shape.size()==2 ) {
96 68484 : if( isconst ) val.get<const T*>({shape[0],shape[1]}); else val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
97 : }
98 447222 : v=val.copy();
99 447222 : }
100 :
101 : template <class T>
102 282703 : void DataPassingObjectTyped<T>::setForcePointer( const TypesafePtr & val, const std::vector<unsigned>& shape ) {
103 282703 : if( shape.size()==0 ) val.get<T*>(); // just check type and discard pointer
104 64795 : else if( shape.size()==1 ) val.get<T*>({shape[0]}); // just check type and discard pointer
105 64795 : else if( shape.size()==2 ) val.get<T*>({shape[0],shape[1]}); // just check type and discard pointer
106 282687 : f=val.copy();
107 282685 : }
108 :
109 : template <class T>
110 12447 : void DataPassingObjectTyped<T>::setData( Value* value ) {
111 12447 : if( value->getRank()==0 ) { *v.template get<T*>() = static_cast<T>(value->get()) / unit; return; }
112 11744 : T* pp; getPointer( v, value->getShape(), start, stride, pp ); unsigned nvals=value->getNumberOfValues();
113 429326 : for(unsigned i=0; i<nvals; ++i) pp[i] = T( value->get(i) );
114 : }
115 :
116 : template <class T>
117 239979 : void DataPassingObjectTyped<T>::share_data( const unsigned& j, const unsigned& k, Value* value ) {
118 239979 : if( value->getRank()==0 ) {
119 7352 : if( hasbackup ) value->set( unit*bvalue );
120 72 : else value->set( unit*double(v.template get<T>()) );
121 7352 : return;
122 : }
123 232627 : std::vector<unsigned> s(value->getShape()); if( s.size()==1 ) s[0]=k-j;
124 232627 : const T* pp; getPointer( v, s, start, stride, pp );
125 232611 : std::vector<double> & d=value->data;
126 232611 : #pragma omp parallel for num_threads(value->getGoodNumThreads(j,k))
127 : for(unsigned i=j; i<k; ++i) d[i]=unit*pp[i*stride];
128 : }
129 :
130 : template <class T>
131 1350 : void DataPassingObjectTyped<T>::share_data( std::vector<double>& values ) const {
132 1350 : std::vector<unsigned> maxel(1,values.size()); const T* pp; getPointer( v, maxel, start, stride, pp );
133 1350 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(values))
134 : for(unsigned i=0; i<values.size(); ++i) values[i]=unit*pp[i*stride];
135 1350 : }
136 :
137 : template <class T>
138 70557 : void DataPassingObjectTyped<T>::share_data( const std::vector<AtomNumber>&index, const std::vector<unsigned>& i, Value* value ) {
139 70557 : plumed_dbg_assert( value->getRank()==1 ); std::vector<unsigned> maxel(1,index.size());
140 : #ifndef NDEBUG
141 : // bounds are only checked in debug mode since they require this extra step that is potentially expensive
142 : maxel[0]=(i.size()>0?*std::max_element(i.begin(),i.end())+1:0);
143 : #else
144 70557 : maxel[0]=0;
145 : #endif
146 70557 : const T* pp; getPointer( v, maxel, start, stride, pp );
147 : // cannot be parallelized with omp because access to data is not ordered
148 2470879 : unsigned k=0; for(const auto & p : index) { value->data[p.index()]=unit*pp[i[k]*stride]; k++; }
149 70557 : }
150 :
151 : template <class T>
152 39905 : void DataPassingObjectTyped<T>::add_force( Value* value ) {
153 39905 : if( value->getRank()==0 ) { *f.template get<T*>() += funit*static_cast<T>(value->getForce(0)); return; }
154 39865 : T* pp; getPointer( f, value->getShape(), start, stride, pp ); unsigned nvals=value->getNumberOfValues();
155 39865 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,nvals))
156 : for(unsigned i=0; i<nvals; ++i) pp[i*stride] += funit*T(value->getForce(i));
157 : }
158 :
159 : template <class T>
160 79589 : void DataPassingObjectTyped<T>::add_force( const std::vector<int>& index, Value* value ) {
161 79589 : plumed_assert( value->getRank()==1 ); std::vector<unsigned> s(1,index.size()); T* pp; getPointer( f, s, start, stride, pp );
162 79587 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,index.size()))
163 : for(unsigned i=0; i<index.size(); ++i) pp[i*stride] += funit*T(value->getForce(index[i]));
164 79587 : }
165 :
166 : template <class T>
167 62508 : void DataPassingObjectTyped<T>::add_force( const std::vector<AtomNumber>& index, const std::vector<unsigned>& i, Value* value ) {
168 62508 : plumed_dbg_assert( value->getRank()==1 ); std::vector<unsigned> maxel(1,index.size());
169 : #ifndef NDEBUG
170 : // bounds are only checked in debug mode since they require this extra step that is potentially expensive
171 : maxel[0]=(i.size()>0?*std::max_element(i.begin(),i.end())+1:0);
172 : #else
173 62508 : maxel[0]=0;
174 : #endif
175 62508 : T* pp; getPointer( f, maxel, start, stride, pp );
176 1844868 : unsigned k=0; for(const auto & p : index) { pp[stride*i[k]] += funit*T(value->getForce(p.index())); k++; }
177 62508 : }
178 :
179 : template <class T>
180 150 : void DataPassingObjectTyped<T>::rescale_force( const unsigned& n, const double& factor, Value* value ) {
181 150 : plumed_assert( value->getRank()>0 ); std::vector<unsigned> s( value->getShape() ); if( s.size()==1 ) s[0] = n;
182 150 : T* pp; getPointer( f, s, start, stride, pp );
183 150 : #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(pp,n))
184 : for(unsigned i=0; i<n; ++i) pp[i*stride] *= T(factor);
185 150 : }
186 :
187 : }
|
