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 191221918 : void VectorGeneric<n>::auxiliaryConstructor(double first,Args... arg) {
164 191221918 : d[n-(sizeof...(Args))-1]=first;
165 131763012 : auxiliaryConstructor(arg...);
166 191221918 : }
167 :
168 : template <unsigned n>
169 : template<typename... Args>
170 59458906 : VectorGeneric<n>::VectorGeneric(double first,Args... arg) {
171 : static_assert((sizeof...(Args))+1==n,"you are trying to initialize a Vector with the wrong number of arguments");
172 59458906 : auxiliaryConstructor(first,arg...);
173 59458906 : }
174 :
175 : template <unsigned n>
176 25037430 : void VectorGeneric<n>::zero() {
177 25037430 : LoopUnroller<n>::_zero(d.data());
178 25037430 : }
179 :
180 : template <unsigned n>
181 646926254 : VectorGeneric<n>::VectorGeneric() {
182 646343536 : LoopUnroller<n>::_zero(d.data());
183 646926254 : }
184 :
185 : template <unsigned n>
186 2361748037 : double & VectorGeneric<n>::operator[](unsigned i) {
187 2361748037 : return d[i];
188 : }
189 :
190 : template <unsigned n>
191 2551360111 : const double & VectorGeneric<n>::operator[](unsigned i)const {
192 2551360111 : return d[i];
193 : }
194 :
195 : template <unsigned n>
196 1581184674 : double & VectorGeneric<n>::operator()(unsigned i) {
197 1581184674 : return d[i];
198 : }
199 :
200 : template <unsigned n>
201 1531216359 : const double & VectorGeneric<n>::operator()(unsigned i)const {
202 1531216359 : return d[i];
203 : }
204 :
205 : template <unsigned n>
206 1960946512 : VectorGeneric<n>& VectorGeneric<n>::operator +=(const VectorGeneric<n>& b) {
207 1960946512 : LoopUnroller<n>::_add(d.data(),b.d.data());
208 1960946512 : return *this;
209 : }
210 :
211 : template <unsigned n>
212 1914905212 : VectorGeneric<n>& VectorGeneric<n>::operator -=(const VectorGeneric<n>& b) {
213 1914905212 : LoopUnroller<n>::_sub(d.data(),b.d.data());
214 1914905212 : return *this;
215 : }
216 :
217 : template <unsigned n>
218 2439672640 : VectorGeneric<n>& VectorGeneric<n>::operator *=(double s) {
219 2439672640 : LoopUnroller<n>::_mul(d.data(),s);
220 2439672640 : return *this;
221 : }
222 :
223 : template <unsigned n>
224 57841 : VectorGeneric<n>& VectorGeneric<n>::operator /=(double s) {
225 57841 : LoopUnroller<n>::_mul(d.data(),1.0/s);
226 57841 : return *this;
227 : }
228 :
229 : template <unsigned n>
230 : VectorGeneric<n> VectorGeneric<n>::operator +()const {
231 : return *this;
232 : }
233 :
234 : template <unsigned n>
235 206458710 : VectorGeneric<n> VectorGeneric<n>::operator -()const {
236 206458710 : VectorGeneric<n> r;
237 206458710 : LoopUnroller<n>::_neg(r.d.data(),d.data());
238 206458710 : return r;
239 : }
240 :
241 : template <unsigned n>
242 1319412804 : VectorGeneric<n> operator+(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
243 1319412804 : VectorGeneric<n> v(v1);
244 1319412804 : return v+=v2;
245 : }
246 :
247 : template <unsigned n>
248 1371788948 : VectorGeneric<n> operator-(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
249 1371788948 : VectorGeneric<n> v(v1);
250 1371788948 : return v-=v2;
251 : }
252 :
253 : template <unsigned n>
254 2437590445 : VectorGeneric<n> operator*(double s,const VectorGeneric<n>&v) {
255 2437590445 : VectorGeneric<n> vv(v);
256 2437590445 : return vv*=s;
257 : }
258 :
259 : template <unsigned n>
260 947466513 : VectorGeneric<n> operator*(const VectorGeneric<n>&v,double s) {
261 947466513 : return s*v;
262 : }
263 :
264 : template <unsigned n>
265 359170206 : VectorGeneric<n> operator/(const VectorGeneric<n>&v,double s) {
266 359170206 : return v*(1.0/s);
267 : }
268 :
269 : template <unsigned n>
270 993455922 : VectorGeneric<n> delta(const VectorGeneric<n>&v1,const VectorGeneric<n>&v2) {
271 993455922 : return v2-v1;
272 : }
273 :
274 : template <unsigned n>
275 1109907003 : double VectorGeneric<n>::modulo2()const {
276 1109907003 : return LoopUnroller<n>::_sum2(d.data());
277 : }
278 :
279 : template <unsigned n>
280 265375907 : double dotProduct(const VectorGeneric<n>& v1,const VectorGeneric<n>& v2) {
281 265375907 : return LoopUnroller<n>::_dot(v1.d.data(),v2.d.data());
282 : }
283 :
284 : inline
285 5459803 : VectorGeneric<3> crossProduct(const VectorGeneric<3>& v1,const VectorGeneric<3>& v2) {
286 : return VectorGeneric<3>(
287 5459803 : v1[1]*v2[2]-v1[2]*v2[1],
288 5459803 : v1[2]*v2[0]-v1[0]*v2[2],
289 5459803 : v1[0]*v2[1]-v1[1]*v2[0]);
290 : }
291 :
292 : template<unsigned n>
293 408021831 : double VectorGeneric<n>::modulo()const {
294 408021831 : return sqrt(modulo2());
295 : }
296 :
297 : template<unsigned n>
298 203397173 : double modulo2(const VectorGeneric<n>&v) {
299 203397173 : return v.modulo2();
300 : }
301 :
302 : template<unsigned n>
303 34471352 : double modulo(const VectorGeneric<n>&v) {
304 34471352 : return v.modulo();
305 : }
306 :
307 : template<unsigned n>
308 0 : std::ostream & operator<<(std::ostream &os, const VectorGeneric<n>& v) {
309 0 : for(unsigned i=0; i<n-1; i++) {
310 0 : os<<v(i)<<" ";
311 : }
312 0 : os<<v(n-1);
313 0 : 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 :
|
