Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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_tools_Vector_h
23 : #define __PLUMED_tools_Vector_h
24 :
25 : #include <cmath>
26 : #include <iosfwd>
27 : #include <array>
28 : #include "LoopUnroller.h"
29 :
30 : namespace PLMD {
31 :
32 : /**
33 : \ingroup TOOLBOX
34 : Class implementing fixed size vectors of doubles
35 :
36 : \tparam n The number of elements of the vector.
37 :
38 : This class implements a vector of doubles with size fixed at
39 : compile time. It is useful for small fixed size objects (e.g.
40 : 3d vectors) as it does not waste space to store the vector size.
41 : Moreover, as the compiler knows the size, it can be completely
42 : opimized inline.
43 : All the methods are inlined for better optimization and
44 : all the loops are explicitly unrolled using PLMD::LoopUnroller class.
45 : Vector elements are initialized to zero by default. Notice that
46 : this means that constructor is a bit slow. This point might change
47 : in future if we find performance issues.
48 : Accepts both [] and () syntax for access.
49 : Several functions are declared as friends even if not necessary so as to
50 : properly appear in Doxygen documentation.
51 :
52 : Aliases are defined to simplify common declarations (Vector, Vector2d, Vector3d, Vector4d).
53 : Also notice that some operations are only available for 3 dimensional vectors.
54 :
55 : Example of usage
56 : \verbatim
57 : #include "Vector.h"
58 :
59 : using namespace PLMD;
60 :
61 : int main(){
62 : VectorGeneric<3> v1;
63 : v1[0]=3.0;
64 : // use equivalently () and [] syntax:
65 : v1(1)=5.0;
66 : // initialize with components
67 : VectorGeneric<3> v2=VectorGeneric<3>(1.0,2.0,3.0);
68 : VectorGeneric<3> v3=crossProduct(v1,v2);
69 : double d=dotProduct(v1,v2);
70 : v3+=v1;
71 : v2=v1+2.0*v3;
72 : }
73 : \endverbatim
74 :
75 : */
76 :
77 :
78 : template <unsigned n>
79 : class VectorGeneric {
80 : std::array<double,n> d;
81 : /// Auxiliary private function for constructor
82 : void auxiliaryConstructor();
83 : /// Auxiliary private function for constructor
84 : template<typename... Args>
85 : void auxiliaryConstructor(double first,Args... arg);
86 : public:
87 : /// Constructor accepting n double parameters.
88 : /// Can be used as Vector<3>(1.0,2.0,3.0) or Vector<2>(2.0,3.0).
89 : /// In case a wrong number of parameters is given, a static assertion will fail.
90 : template<typename... Args>
91 : VectorGeneric(double first,Args... arg);
92 : /// create it null
93 : VectorGeneric();
94 : /// set it to zero
95 : void zero();
96 : /// array-like access [i]
97 : double & operator[](unsigned i);
98 : /// array-like access [i]
99 : const double & operator[](unsigned i)const;
100 : /// parenthesis access (i)
101 : double & operator()(unsigned i);
102 : /// parenthesis access (i)
103 : const double & operator()(unsigned i)const;
104 : /// increment
105 : VectorGeneric& operator +=(const VectorGeneric& b);
106 : /// decrement
107 : VectorGeneric& operator -=(const VectorGeneric& b);
108 : /// multiply
109 : VectorGeneric& operator *=(double s);
110 : /// divide
111 : VectorGeneric& operator /=(double s);
112 : /// sign +
113 : VectorGeneric operator +()const;
114 : /// sign -
115 : VectorGeneric operator -()const;
116 : /// return v1+v2
117 : template<unsigned m>
118 : friend VectorGeneric<m> operator+(const VectorGeneric<m>&,const VectorGeneric<m>&);
119 : /// return v1-v2
120 : template<unsigned m>
121 : friend VectorGeneric<m> operator-(const VectorGeneric<m>&,const VectorGeneric<m>&);
122 : /// return s*v
123 : template<unsigned m>
124 : friend VectorGeneric<m> operator*(double,const VectorGeneric<m>&);
125 : /// return v*s
126 : template<unsigned m>
127 : friend VectorGeneric<m> operator*(const VectorGeneric<m>&,double);
128 : /// return v/s
129 : template<unsigned m>
130 : friend VectorGeneric<m> operator/(const VectorGeneric<m>&,double);
131 : /// return v2-v1
132 : template<unsigned m>
133 : friend VectorGeneric<m> delta(const VectorGeneric<m>&v1,const VectorGeneric<m>&v2);
134 : /// return v1 .scalar. v2
135 : template<unsigned m>
136 : friend double dotProduct(const VectorGeneric<m>&,const VectorGeneric<m>&);
137 : /// return v1 .vector. v2
138 : /// Only available for size 3
139 : friend VectorGeneric<3> crossProduct(const VectorGeneric<3>&,const VectorGeneric<3>&);
140 : /// compute the squared modulo
141 : double modulo2()const;
142 : /// Compute the modulo.
143 : /// Shortcut for sqrt(v.modulo2())
144 : double modulo()const;
145 : /// friend version of modulo2 (to simplify some syntax)
146 : template<unsigned m>
147 : friend double modulo2(const VectorGeneric<m>&);
148 : /// friend version of modulo (to simplify some syntax)
149 : template<unsigned m>
150 : friend double modulo(const VectorGeneric<m>&);
151 : /// << operator.
152 : /// Allows printing vector `v` with `std::cout<<v;`
153 : template<unsigned m>
154 : friend std::ostream & operator<<(std::ostream &os, const VectorGeneric<m>&);
155 : };
156 :
157 : template <unsigned n>
158 : void VectorGeneric<n>::auxiliaryConstructor()
159 : {}
160 :
161 : template <unsigned n>
162 : template<typename... Args>
163 119602528 : void VectorGeneric<n>::auxiliaryConstructor(double first,Args... arg)
164 : {
165 119602528 : d[n-(sizeof...(Args))-1]=first;
166 79920350 : auxiliaryConstructor(arg...);
167 119602528 : }
168 :
169 : template <unsigned n>
170 : template<typename... Args>
171 39682178 : VectorGeneric<n>::VectorGeneric(double first,Args... arg)
172 : {
173 : static_assert((sizeof...(Args))+1==n,"you are trying to initialize a Vector with the wrong number of arguments");
174 39682178 : auxiliaryConstructor(first,arg...);
175 39682178 : }
176 :
177 : template <unsigned n>
178 12803723 : void VectorGeneric<n>::zero() {
179 12803723 : LoopUnroller<n>::_zero(d.data());
180 12803723 : }
181 :
182 : template <unsigned n>
183 475870486 : VectorGeneric<n>::VectorGeneric() {
184 475321566 : LoopUnroller<n>::_zero(d.data());
185 475870486 : }
186 :
187 : template <unsigned n>
188 2899777226 : double & VectorGeneric<n>::operator[](unsigned i) {
189 2899777226 : return d[i];
190 : }
191 :
192 : template <unsigned n>
193 3779351632 : const double & VectorGeneric<n>::operator[](unsigned i)const {
194 3779351632 : return d[i];
195 : }
196 :
197 : template <unsigned n>
198 1261457172 : double & VectorGeneric<n>::operator()(unsigned i) {
199 1261457172 : return d[i];
200 : }
201 :
202 : template <unsigned n>
203 1260767808 : const double & VectorGeneric<n>::operator()(unsigned i)const {
204 1260767808 : return d[i];
205 : }
206 :
207 : template <unsigned n>
208 1622467416 : VectorGeneric<n>& VectorGeneric<n>::operator +=(const VectorGeneric<n>& b) {
209 1622467416 : LoopUnroller<n>::_add(d.data(),b.d.data());
210 1622467416 : return *this;
211 : }
212 :
213 : template <unsigned n>
214 1049811592 : VectorGeneric<n>& VectorGeneric<n>::operator -=(const VectorGeneric<n>& b) {
215 1049811592 : LoopUnroller<n>::_sub(d.data(),b.d.data());
216 1049811592 : return *this;
217 : }
218 :
219 : template <unsigned n>
220 1659864777 : VectorGeneric<n>& VectorGeneric<n>::operator *=(double s) {
221 1659864777 : LoopUnroller<n>::_mul(d.data(),s);
222 1659864777 : return *this;
223 : }
224 :
225 : template <unsigned n>
226 8301 : VectorGeneric<n>& VectorGeneric<n>::operator /=(double s) {
227 8301 : LoopUnroller<n>::_mul(d.data(),1.0/s);
228 8301 : return *this;
229 : }
230 :
231 : template <unsigned n>
232 : VectorGeneric<n> VectorGeneric<n>::operator +()const {
233 : return *this;
234 : }
235 :
236 : template <unsigned n>
237 106091865 : VectorGeneric<n> VectorGeneric<n>::operator -()const {
238 106091865 : VectorGeneric<n> r;
239 106091865 : LoopUnroller<n>::_neg(r.d.data(),d.data());
240 106091865 : return r;
241 : }
242 :
243 : template <unsigned n>
244 1368903652 : VectorGeneric<n> operator+(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
245 1368903652 : VectorGeneric<n> v(v1);
246 1368903652 : return v+=v2;
247 : }
248 :
249 : template <unsigned n>
250 854663986 : VectorGeneric<n> operator-(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
251 854663986 : VectorGeneric<n> v(v1);
252 854663986 : return v-=v2;
253 : }
254 :
255 : template <unsigned n>
256 1659517798 : VectorGeneric<n> operator*(double s,const VectorGeneric<n>&v) {
257 1659517798 : VectorGeneric<n> vv(v);
258 1659517798 : return vv*=s;
259 : }
260 :
261 : template <unsigned n>
262 281949147 : VectorGeneric<n> operator*(const VectorGeneric<n>&v,double s) {
263 281949147 : return s*v;
264 : }
265 :
266 : template <unsigned n>
267 40818496 : VectorGeneric<n> operator/(const VectorGeneric<n>&v,double s) {
268 40818496 : return v*(1.0/s);
269 : }
270 :
271 : template <unsigned n>
272 318536896 : VectorGeneric<n> delta(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
273 318536896 : return v2-v1;
274 : }
275 :
276 : template <unsigned n>
277 739774321 : double VectorGeneric<n>::modulo2()const {
278 739774321 : return LoopUnroller<n>::_sum2(d.data());
279 : }
280 :
281 : template <unsigned n>
282 166313906 : double dotProduct(const VectorGeneric<n>& v1,const VectorGeneric<n>& v2) {
283 166313906 : return LoopUnroller<n>::_dot(v1.d.data(),v2.d.data());
284 : }
285 :
286 : inline
287 4881117 : VectorGeneric<3> crossProduct(const VectorGeneric<3>& v1,const VectorGeneric<3>& v2) {
288 : return VectorGeneric<3>(
289 4881117 : v1[1]*v2[2]-v1[2]*v2[1],
290 4881117 : v1[2]*v2[0]-v1[0]*v2[2],
291 4881117 : v1[0]*v2[1]-v1[1]*v2[0]);
292 : }
293 :
294 : template<unsigned n>
295 84431275 : double VectorGeneric<n>::modulo()const {
296 84431275 : return sqrt(modulo2());
297 : }
298 :
299 : template<unsigned n>
300 191884672 : double modulo2(const VectorGeneric<n>&v) {
301 191884672 : return v.modulo2();
302 : }
303 :
304 : template<unsigned n>
305 260 : double modulo(const VectorGeneric<n>&v) {
306 260 : return v.modulo();
307 : }
308 :
309 : template<unsigned n>
310 : std::ostream & operator<<(std::ostream &os, const VectorGeneric<n>& v) {
311 : for(unsigned i=0; i<n-1; i++) os<<v(i)<<" ";
312 : os<<v(n-1);
313 : return os;
314 : }
315 :
316 :
317 : /// \ingroup TOOLBOX
318 : /// Alias for one dimensional vectors
319 : typedef VectorGeneric<1> Vector1d;
320 : /// \ingroup TOOLBOX
321 : /// Alias for two dimensional vectors
322 : typedef VectorGeneric<2> Vector2d;
323 : /// \ingroup TOOLBOX
324 : /// Alias for three dimensional vectors
325 : typedef VectorGeneric<3> Vector3d;
326 : /// \ingroup TOOLBOX
327 : /// Alias for four dimensional vectors
328 : typedef VectorGeneric<4> Vector4d;
329 : /// \ingroup TOOLBOX
330 : /// Alias for five dimensional vectors
331 : typedef VectorGeneric<5> Vector5d;
332 : /// \ingroup TOOLBOX
333 : /// Alias for three dimensional vectors
334 : typedef Vector3d Vector;
335 :
336 : static_assert(sizeof(VectorGeneric<2>)==2*sizeof(double), "code cannot work if this is not satisfied");
337 : static_assert(sizeof(VectorGeneric<3>)==3*sizeof(double), "code cannot work if this is not satisfied");
338 : static_assert(sizeof(VectorGeneric<4>)==4*sizeof(double), "code cannot work if this is not satisfied");
339 :
340 : }
341 :
342 : #endif
343 :
|
