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 <initializer_list>
27 : #include <iosfwd>
28 : #include <array>
29 : #include <type_traits>
30 : #include "LoopUnroller.h"
31 :
32 : namespace PLMD {
33 :
34 : /**
35 : \ingroup TOOLBOX
36 : Class implementing fixed size vectors of doubles
37 :
38 : \tparam n The number of elements of the vector.
39 :
40 : This class implements a vector of doubles with size fixed at
41 : compile time. It is useful for small fixed size objects (e.g.
42 : 3d vectors) as it does not waste space to store the vector size.
43 : Moreover, as the compiler knows the size, it can be completely
44 : optimized inline.
45 : All the methods are inlined for better optimization and
46 : all the loops are explicitly unrolled using PLMD::LoopUnroller class.
47 : Vector elements are initialized to zero by default. Notice that
48 : this means that constructor is a bit slow. This point might change
49 : in future if we find performance issues.
50 : Accepts both [] and () syntax for access.
51 : Several functions are declared as friends even if not necessary so as to
52 : properly appear in Doxygen documentation.
53 :
54 : Aliases are defined to simplify common declarations (Vector, Vector2d, Vector3d, Vector4d).
55 : Also notice that some operations are only available for 3 dimensional vectors.
56 :
57 : Example of usage
58 : \verbatim
59 : #include "Vector.h"
60 :
61 : using namespace PLMD;
62 :
63 : int main(){
64 : VectorGeneric<3> v1;
65 : v1[0]=3.0;
66 : // use equivalently () and [] syntax:
67 : v1(1)=5.0;
68 : // initialize with components
69 : VectorGeneric<3> v2=VectorGeneric<3>(1.0,2.0,3.0);
70 : VectorGeneric<3> v3=crossProduct(v1,v2);
71 : double d=dotProduct(v1,v2);
72 : v3+=v1;
73 : v2=v1+2.0*v3;
74 : }
75 : \endverbatim
76 :
77 : */
78 : template<typename T, unsigned n> class VectorTyped;
79 :
80 : template<typename T, unsigned n>
81 : constexpr VectorTyped<T, n> delta(const VectorTyped<T, n>&,const VectorTyped<T, n>&);
82 :
83 : template<typename T, unsigned n>
84 : constexpr T dotProduct(const VectorTyped<T, n>&,const VectorTyped<T, n>&);
85 :
86 :
87 : template<typename T, unsigned n>
88 : std::ostream & operator<<(std::ostream &os, const VectorTyped<T, n>& v);
89 : template<typename T, unsigned n>
90 : class VectorTyped {
91 : std::array<T,n> d;
92 : /// Auxiliary private function for constructor
93 : void auxiliaryConstructor();
94 : /// Auxiliary private function for constructor
95 : template<typename... Args>
96 : void auxiliaryConstructor(T first,Args... arg);
97 : public:
98 : template<unsigned m=n,typename=std::enable_if_t<m==3,void>>
99 304461615 : constexpr VectorTyped(T first, T second, T third):d{first,second,third} {}
100 : /// Constructor accepting n T parameters.
101 : /// Can be used as Vector<3>(1.0,2.0,3.0) or Vector<2>(2.0,3.0).
102 : /// In case a wrong number of parameters is given, a static assertion will fail.
103 : template<typename... Args,unsigned m=n,typename=std::enable_if_t<m!=3,void>>
104 : VectorTyped(T first,Args... arg);
105 : /// create it null
106 : constexpr VectorTyped();
107 : /// Returns a pointer to the underlying array serving as element storage.
108 : constexpr T* data() noexcept;
109 : /// Returns a pointer to the underlying array serving as element storage.
110 : constexpr const T* data() const noexcept;
111 : /// set it to zero
112 : constexpr void zero();
113 : /// array-like access [i]
114 : constexpr T & operator[](unsigned i);
115 : /// array-like access [i]
116 : constexpr const T & operator[](unsigned i)const;
117 : /// parenthesis access (i)
118 : constexpr T & operator()(unsigned i);
119 : /// parenthesis access (i)
120 : constexpr const T & operator()(unsigned i)const;
121 : /// increment
122 : constexpr VectorTyped& operator +=(const VectorTyped& b);
123 : /// decrement
124 : constexpr VectorTyped& operator -=(const VectorTyped& b);
125 : /// multiply
126 : constexpr VectorTyped& operator *=(T s);
127 : /// divide
128 : constexpr VectorTyped& operator /=(T s);
129 : /// sign +
130 : constexpr VectorTyped operator +()const;
131 : /// sign -
132 : constexpr VectorTyped operator -()const;
133 : /// return v1+v2
134 : template<typename U, unsigned m>
135 : friend constexpr VectorTyped<U, m> operator+(const VectorTyped<U, m>&,const VectorTyped<U, m>&);
136 : /// return v1-v2
137 : template<typename U, unsigned m>
138 : friend constexpr VectorTyped<U, m> operator-(VectorTyped<U, m>,const VectorTyped<U, m>&);
139 : /// return s*v
140 : template<typename U, typename J, unsigned m>
141 : friend constexpr VectorTyped<U, m> operator*(J,VectorTyped<U, m>);
142 : /// return v*s
143 : template<typename U, typename J, unsigned m>
144 : friend constexpr VectorTyped<U, m> operator*(VectorTyped<U, m>,J);
145 : /// return v/s
146 : template<typename U, typename J, unsigned m>
147 : friend constexpr VectorTyped<U, m> operator/(const VectorTyped<U, m>&,J);
148 : /// return v2-v1
149 : friend constexpr VectorTyped delta<>(const VectorTyped&v1,const VectorTyped&v2);
150 : /// return v1 .scalar. v2
151 : friend constexpr T dotProduct<>(const VectorTyped&,const VectorTyped&);
152 : //this bad boy produces a warning (in fact because declrare the crossproduc as a friend for ALL thhe possible combinations of n and T)
153 : /// return v1 .vector. v2
154 : /// Only available for size 3
155 : template<typename U>
156 : friend constexpr VectorTyped<U, 3> crossProduct(const VectorTyped<U, 3>&,const VectorTyped<U, 3>&);
157 : /// compute the squared modulo
158 : constexpr T modulo2()const;
159 : /// Compute the modulo.
160 : /// Shortcut for sqrt(v.modulo2())
161 : constexpr T modulo()const;
162 : /// friend version of modulo2 (to simplify some syntax)
163 : template<typename U, unsigned m>
164 : friend constexpr U modulo2(const VectorTyped<U, m>&);
165 : /// friend version of modulo (to simplify some syntax)
166 : template<typename U, unsigned m>
167 : friend constexpr U modulo(const VectorTyped<U, m>&);
168 : ///conversion to another type
169 : template<typename TT>
170 : constexpr VectorTyped<TT,n> convert()const;
171 : /// << operator.
172 : /// Allows printing vector `v` with `std::cout<<v;`
173 : friend std::ostream & operator<< <>(std::ostream &os, const VectorTyped&);
174 : };
175 :
176 : template<typename T, unsigned n>
177 : void VectorTyped<T, n>::auxiliaryConstructor()
178 : {}
179 :
180 : template<typename T, unsigned n>
181 : template<typename... Args>
182 35691828 : void VectorTyped<T, n>::auxiliaryConstructor(T first,Args... arg) {
183 35691828 : d[n-(sizeof...(Args))-1]=first;
184 26768935 : auxiliaryConstructor(arg...);
185 35691828 : }
186 :
187 : template<typename T, unsigned n>
188 : template<typename... Args, unsigned,typename >
189 8922893 : VectorTyped<T, n>::VectorTyped(T first,Args... arg) {
190 : static_assert((sizeof...(Args))+1==n,"you are trying to initialize a Vector with the wrong number of arguments");
191 8922893 : auxiliaryConstructor(first,arg...);
192 8922893 : }
193 :
194 : template<typename T, unsigned n>
195 : constexpr T* VectorTyped<T, n>::data() noexcept {
196 : return d.data();
197 : }
198 :
199 : template<typename T, unsigned n>
200 : constexpr const T* VectorTyped<T, n>::data() const noexcept {
201 : return d.data();
202 : }
203 :
204 : template<typename T, unsigned n>
205 : constexpr void VectorTyped<T, n>::zero() {
206 : LoopUnroller<n>::_zero(d.data());
207 : }
208 :
209 : template<typename T, unsigned n>
210 79557183 : constexpr VectorTyped<T, n>::VectorTyped() {
211 : LoopUnroller<n>::_zero(d.data());
212 : }
213 :
214 : template<typename T, unsigned n>
215 : constexpr T & VectorTyped<T, n>::operator[](unsigned i) {
216 1849708081 : return d[i];
217 : }
218 :
219 : template<typename T, unsigned n>
220 : constexpr const T & VectorTyped<T, n>::operator[](unsigned i)const {
221 1174954671 : return d[i];
222 : }
223 :
224 : template<typename T, unsigned n>
225 : constexpr T & VectorTyped<T, n>::operator()(unsigned i) {
226 1236395736 : return d[i];
227 : }
228 :
229 : template<typename T, unsigned n>
230 : constexpr const T & VectorTyped<T, n>::operator()(unsigned i)const {
231 1877001820 : return d[i];
232 : }
233 :
234 : template<typename T, unsigned n>
235 : constexpr VectorTyped<T, n>& VectorTyped<T, n>::operator +=(const VectorTyped<T, n>& b) {
236 3106120 : LoopUnroller<n>::_add(d.data(),b.d.data());
237 618 : return *this;
238 : }
239 :
240 : template<typename T, unsigned n>
241 : constexpr VectorTyped<T, n>& VectorTyped<T, n>::operator -=(const VectorTyped<T, n>& b) {
242 1119126 : LoopUnroller<n>::_sub(d.data(),b.d.data());
243 : return *this;
244 : }
245 :
246 : template<typename T, unsigned n>
247 : constexpr VectorTyped<T, n>& VectorTyped<T, n>::operator *=(T s) {
248 128 : LoopUnroller<n>::_mul(d.data(),s);
249 : return *this;
250 : }
251 :
252 : template<typename T, unsigned n>
253 : constexpr VectorTyped<T, n>& VectorTyped<T, n>::operator /=(T s) {
254 862722 : LoopUnroller<n>::_mul(d.data(),1.0/s);
255 60 : return *this;
256 : }
257 :
258 : template<typename T, unsigned n>
259 : constexpr VectorTyped<T, n> VectorTyped<T, n>::operator +()const {
260 : return *this;
261 : }
262 :
263 : template<typename T, unsigned n>
264 : constexpr VectorTyped<T, n> VectorTyped<T, n>::operator -()const {
265 : VectorTyped<T, n> r;
266 : LoopUnroller<n>::_neg(r.d.data(),d.data());
267 144558 : return r;
268 : }
269 :
270 : template<typename T, unsigned n>
271 3106120 : constexpr VectorTyped<T, n> operator+(const VectorTyped<T, n>&v1,const VectorTyped<T, n>&v2) {
272 14483348 : VectorTyped<T, n> v(v1);
273 3190316 : return v+=v2;
274 : }
275 :
276 : template<typename T, typename TT, unsigned n>
277 : constexpr VectorTyped<T, n> sumT(const VectorTyped<T, n>&v1,const VectorTyped<TT, n>&v2) {
278 : VectorTyped<T, n> v(v1);
279 : LoopUnroller<n>::_add(v.data(),v2.data());
280 : return v;
281 : }
282 :
283 : template<typename T, unsigned n>
284 : constexpr VectorTyped<T, n> operator-(VectorTyped<T, n>v1,const VectorTyped<T, n>&v2) {
285 1055258265 : return v1-=v2;
286 : }
287 :
288 : template<typename T, typename J, unsigned n>
289 : constexpr VectorTyped<T, n> operator*(J s,VectorTyped<T, n>v) {
290 418125668 : return v*=s;
291 : }
292 :
293 : template<typename T, typename J, unsigned n>
294 : constexpr VectorTyped<T, n> operator*(VectorTyped<T, n> v,J s) {
295 7324152 : return v*=s;
296 : }
297 :
298 : template<typename T, typename J, unsigned n>
299 : constexpr VectorTyped<T, n> operator/(const VectorTyped<T, n>&v,J s) {
300 244879665 : return v*(T(1.0)/s);
301 : }
302 :
303 : template<typename T, unsigned n>
304 1119126 : constexpr VectorTyped<T, n> delta(const VectorTyped<T, n>&v1,const VectorTyped<T, n>&v2) {
305 1119741 : return v2-v1;
306 : }
307 :
308 : template<typename T, unsigned n>
309 : constexpr T VectorTyped<T, n>::modulo2()const {
310 : return LoopUnroller<n>::_sum2(d.data());
311 : }
312 :
313 : template<typename T, unsigned n>
314 : constexpr T dotProduct(const VectorTyped<T, n>& v1,const VectorTyped<T, n>& v2) {
315 11456290 : return LoopUnroller<n>::_dot(v1.d.data(),v2.d.data());
316 : }
317 :
318 : template<typename T>
319 : constexpr inline
320 4990426 : VectorTyped<T, 3> crossProduct(const VectorTyped<T, 3>& v1,const VectorTyped<T, 3>& v2) {
321 : return VectorTyped<T, 3>(
322 4990426 : v1[1]*v2[2]-v1[2]*v2[1],
323 4990426 : v1[2]*v2[0]-v1[0]*v2[2],
324 4990426 : v1[0]*v2[1]-v1[1]*v2[0]);
325 : }
326 :
327 : template<typename T, unsigned n>
328 440225230 : constexpr T VectorTyped<T, n>::modulo()const {
329 440225230 : return sqrt(modulo2());
330 : }
331 :
332 : template<typename T, unsigned n>
333 : constexpr T modulo2(const VectorTyped<T, n>&v) {
334 : return v.modulo2();
335 : }
336 :
337 : template<typename T, unsigned n>
338 : constexpr T modulo(const VectorTyped<T, n>&v) {
339 776790 : return v.modulo();
340 : }
341 :
342 : template<typename T, unsigned n>
343 : template<typename TT>
344 : constexpr VectorTyped<TT,n> VectorTyped<T,n>::convert()const {
345 : VectorTyped<TT,n> target;
346 : LoopUnroller<n>::_copy(target.data(),d.data());
347 : return target;
348 : }
349 :
350 : template<typename T, unsigned n>
351 54 : std::ostream & operator<<(std::ostream &os, const VectorTyped<T, n>& v) {
352 162 : for(unsigned i=0; i<n-1; i++) {
353 108 : os<<v(i)<<" ";
354 : }
355 54 : os<<v(n-1);
356 54 : return os;
357 : }
358 :
359 : template<unsigned n>
360 : using VectorGeneric=VectorTyped<double, n>;
361 :
362 : /// \ingroup TOOLBOX
363 : /// Alias for one dimensional vectors
364 : typedef VectorGeneric<1> Vector1d;
365 : /// \ingroup TOOLBOX
366 : /// Alias for two dimensional vectors
367 : typedef VectorGeneric<2> Vector2d;
368 : /// \ingroup TOOLBOX
369 : /// Alias for three dimensional vectors
370 : typedef VectorGeneric<3> Vector3d;
371 : /// \ingroup TOOLBOX
372 : /// Alias for four dimensional vectors
373 : typedef VectorGeneric<4> Vector4d;
374 : /// \ingroup TOOLBOX
375 : /// Alias for five dimensional vectors
376 : typedef VectorGeneric<5> Vector5d;
377 : /// \ingroup TOOLBOX
378 : /// Alias for three dimensional vectors
379 : typedef Vector3d Vector;
380 :
381 : /// \ingroup TOOLBOX
382 : /// Alias for one dimensional vectors
383 : typedef VectorTyped<float,1> Vector1f;
384 : /// \ingroup TOOLBOX
385 : /// Alias for two dimensional vectors
386 : typedef VectorTyped<float,2> Vector2f;
387 : /// \ingroup TOOLBOX
388 : /// Alias for three dimensional vectors
389 : typedef VectorTyped<float,3> Vector3f;
390 : /// \ingroup TOOLBOX
391 : /// Alias for four dimensional vectors
392 : typedef VectorTyped<float,4> Vector4f;
393 : /// \ingroup TOOLBOX
394 : /// Alias for five dimensional vectors
395 : typedef VectorTyped<float,5> Vector5f;
396 :
397 : /// \ingroup TOOLBOX
398 : /// Alias for three dimensional vectors of any type:
399 : template <typename T>
400 : using VectorT=VectorTyped<T,3>;
401 :
402 : static_assert(sizeof(VectorGeneric<2>)==2*sizeof(double), "code cannot work if this is not satisfied");
403 : static_assert(sizeof(VectorGeneric<3>)==3*sizeof(double), "code cannot work if this is not satisfied");
404 : static_assert(sizeof(VectorGeneric<4>)==4*sizeof(double), "code cannot work if this is not satisfied");
405 : namespace Versors {
406 : template <typename T>
407 : constexpr auto xp= PLMD::VectorTyped<T,3>( {
408 : T( 1),T( 0),T( 0)
409 : });
410 : template <typename T>
411 : constexpr auto xm= PLMD::VectorTyped<T,3>( {
412 : T(-1),T( 0),T( 0)
413 : });
414 : template <typename T>
415 : constexpr auto yp= PLMD::VectorTyped<T,3>( {
416 : T( 0),T( 1),T( 0)
417 : });
418 : template <typename T>
419 : constexpr auto ym= PLMD::VectorTyped<T,3>( {
420 : T( 0),T(-1),T( 0)
421 : });
422 : template <typename T>
423 : constexpr auto zp= PLMD::VectorTyped<T,3>( {
424 : T( 0),T( 0),T( 1)
425 : });
426 : template <typename T>
427 : constexpr auto zm= PLMD::VectorTyped<T,3>( {
428 : T( 0),T( 0),T(-1)
429 : });
430 : } // namespace Versors
431 : } //PLMD
432 :
433 : #endif //__PLUMED_tools_Vector_h
|
