1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2011-2021 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>
auxiliaryConstructor()158 void VectorGeneric<n>::auxiliaryConstructor()
159 {}
160
161 template <unsigned n>
162 template<typename... Args>
auxiliaryConstructor(double first,Args...arg)163 void VectorGeneric<n>::auxiliaryConstructor(double first,Args... arg)
164 {
165 d[n-(sizeof...(Args))-1]=first;
166 auxiliaryConstructor(arg...);
167 }
168
169 template <unsigned n>
170 template<typename... Args>
VectorGeneric(double first,Args...arg)171 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 auxiliaryConstructor(first,arg...);
175 }
176
177 template <unsigned n>
zero()178 void VectorGeneric<n>::zero() {
179 LoopUnroller<n>::_zero(d.data());
180 }
181
182 template <unsigned n>
VectorGeneric()183 VectorGeneric<n>::VectorGeneric() {
184 LoopUnroller<n>::_zero(d.data());
185 }
186
187 template <unsigned n>
188 double & VectorGeneric<n>::operator[](unsigned i) {
189 return d[i];
190 }
191
192 template <unsigned n>
193 const double & VectorGeneric<n>::operator[](unsigned i)const {
194 return d[i];
195 }
196
197 template <unsigned n>
operator()198 double & VectorGeneric<n>::operator()(unsigned i) {
199 return d[i];
200 }
201
202 template <unsigned n>
operator()203 const double & VectorGeneric<n>::operator()(unsigned i)const {
204 return d[i];
205 }
206
207 template <unsigned n>
208 VectorGeneric<n>& VectorGeneric<n>::operator +=(const VectorGeneric<n>& b) {
209 LoopUnroller<n>::_add(d.data(),b.d.data());
210 return *this;
211 }
212
213 template <unsigned n>
214 VectorGeneric<n>& VectorGeneric<n>::operator -=(const VectorGeneric<n>& b) {
215 LoopUnroller<n>::_sub(d.data(),b.d.data());
216 return *this;
217 }
218
219 template <unsigned n>
220 VectorGeneric<n>& VectorGeneric<n>::operator *=(double s) {
221 LoopUnroller<n>::_mul(d.data(),s);
222 return *this;
223 }
224
225 template <unsigned n>
226 VectorGeneric<n>& VectorGeneric<n>::operator /=(double s) {
227 LoopUnroller<n>::_mul(d.data(),1.0/s);
228 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 VectorGeneric<n> VectorGeneric<n>::operator -()const {
238 VectorGeneric<n> r;
239 LoopUnroller<n>::_neg(r.d.data(),d.data());
240 return r;
241 }
242
243 template <unsigned n>
244 VectorGeneric<n> operator+(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
245 VectorGeneric<n> v(v1);
246 return v+=v2;
247 }
248
249 template <unsigned n>
250 VectorGeneric<n> operator-(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
251 VectorGeneric<n> v(v1);
252 return v-=v2;
253 }
254
255 template <unsigned n>
256 VectorGeneric<n> operator*(double s,const VectorGeneric<n>&v) {
257 VectorGeneric<n> vv(v);
258 return vv*=s;
259 }
260
261 template <unsigned n>
262 VectorGeneric<n> operator*(const VectorGeneric<n>&v,double s) {
263 return s*v;
264 }
265
266 template <unsigned n>
267 VectorGeneric<n> operator/(const VectorGeneric<n>&v,double s) {
268 return v*(1.0/s);
269 }
270
271 template <unsigned n>
delta(const VectorGeneric<n> & v1,const VectorGeneric<n> & v2)272 VectorGeneric<n> delta(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
273 return v2-v1;
274 }
275
276 template <unsigned n>
modulo2()277 double VectorGeneric<n>::modulo2()const {
278 return LoopUnroller<n>::_sum2(d.data());
279 }
280
281 template <unsigned n>
dotProduct(const VectorGeneric<n> & v1,const VectorGeneric<n> & v2)282 double dotProduct(const VectorGeneric<n>& v1,const VectorGeneric<n>& v2) {
283 return LoopUnroller<n>::_dot(v1.d.data(),v2.d.data());
284 }
285
286 inline
crossProduct(const VectorGeneric<3> & v1,const VectorGeneric<3> & v2)287 VectorGeneric<3> crossProduct(const VectorGeneric<3>& v1,const VectorGeneric<3>& v2) {
288 return VectorGeneric<3>(
289 v1[1]*v2[2]-v1[2]*v2[1],
290 v1[2]*v2[0]-v1[0]*v2[2],
291 v1[0]*v2[1]-v1[1]*v2[0]);
292 }
293
294 template<unsigned n>
modulo()295 double VectorGeneric<n>::modulo()const {
296 return sqrt(modulo2());
297 }
298
299 template<unsigned n>
modulo2(const VectorGeneric<n> & v)300 double modulo2(const VectorGeneric<n>&v) {
301 return v.modulo2();
302 }
303
304 template<unsigned n>
modulo(const VectorGeneric<n> & v)305 double modulo(const VectorGeneric<n>&v) {
306 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 }
337
338 #endif
339
340