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