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_Tensor_h
23 : #define __PLUMED_tools_Tensor_h
24 :
25 : #include "MatrixSquareBracketsAccess.h"
26 : #include "Vector.h"
27 : #include "LoopUnroller.h"
28 :
29 : #ifdef _GLIBCXX_DEBUG
30 : #include "Exception.h"
31 : #endif
32 :
33 : namespace PLMD {
34 :
35 : /**
36 : \ingroup TOOLBOX
37 : Class implementing fixed size matrices of doubles
38 :
39 : \tparam n The number rows
40 : \tparam m The number columns
41 :
42 : This class implements a matrix of doubles with size fixed at
43 : compile time. It is useful for small fixed size objects (e.g.
44 : 3x3 tensors) 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 : Most of the loops are explicitly unrolled using PLMD::LoopUnroller class
48 : Matrix elements are initialized to zero by default. Notice that
49 : this means that constructor is a bit slow. This point might change
50 : in future if we find performance issues.
51 : It takes advantage of MatrixSquareBracketsAccess to provide both
52 : () 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 (Tensor, Tensor2d, Tensor3d, Tensor4d).
57 : Also notice that some operations are only available for 3x3 tensors.
58 :
59 : Example of usage
60 : \verbatim
61 :
62 : #include "Tensor.h"
63 :
64 : using namespace PLMD;
65 :
66 : int main(){
67 : Tensor a;
68 : TensorGeneric<3,2> b;
69 : TensorGeneric<3,2> c=matmul(a,b);
70 : return 0;
71 : }
72 :
73 : \endverbatim
74 : */
75 : template <unsigned n,unsigned m>
76 : class TensorGeneric:
77 : public MatrixSquareBracketsAccess<TensorGeneric<n,m>,double>
78 : {
79 : double d[n*m];
80 : public:
81 : /// initialize the tensor to zero
82 : TensorGeneric();
83 : /// initialize a tensor as an external product of two Vector
84 : TensorGeneric(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2);
85 : /// initialize a tensor with 4 values, in standard C order
86 : TensorGeneric(double,double,double,double);
87 : /// initialize a tensor with 9 values, in standard C order
88 : TensorGeneric(double,double,double,double,double,double,double,double,double);
89 : /// set it to zero
90 : void zero();
91 : /// access element
92 : double & operator() (unsigned i,unsigned j);
93 : /// access element
94 : const double & operator() (unsigned i,unsigned j)const;
95 : /// increment
96 : TensorGeneric& operator +=(const TensorGeneric<n,m>& b);
97 : /// decrement
98 : TensorGeneric& operator -=(const TensorGeneric<n,m>& b);
99 : /// multiply
100 : TensorGeneric& operator *=(double);
101 : /// divide
102 : TensorGeneric& operator /=(double);
103 : /// return +t
104 : TensorGeneric operator +()const;
105 : /// return -t
106 : TensorGeneric operator -()const;
107 : /// set j-th column
108 : TensorGeneric& setCol(unsigned j,const VectorGeneric<n> & c);
109 : /// set i-th row
110 : TensorGeneric& setRow(unsigned i,const VectorGeneric<m> & r);
111 : /// get j-th column
112 : VectorGeneric<n> getCol(unsigned j)const;
113 : /// get i-th row
114 : VectorGeneric<m> getRow(unsigned i)const;
115 : /// return t1+t2
116 : template<unsigned n_,unsigned m_>
117 : friend TensorGeneric<n_,m_> operator+(const TensorGeneric<n_,m_>&,const TensorGeneric<n_,m_>&);
118 : /// return t1+t2
119 : template<unsigned n_,unsigned m_>
120 : friend TensorGeneric<n_,m_> operator-(const TensorGeneric<n_,m_>&,const TensorGeneric<n_,m_>&);
121 : /// scale the tensor by a factor s
122 : template<unsigned n_,unsigned m_>
123 : friend TensorGeneric<n_,m_> operator*(double,const TensorGeneric<n_,m_>&);
124 : /// scale the tensor by a factor s
125 : template<unsigned n_,unsigned m_>
126 : friend TensorGeneric<n_,m_> operator*(const TensorGeneric<n_,m_>&,double s);
127 : /// scale the tensor by a factor 1/s
128 : template<unsigned n_,unsigned m_>
129 : friend TensorGeneric<n_,m_> operator/(const TensorGeneric<n_,m_>&,double s);
130 : /// returns the determinant
131 : double determinant()const;
132 : /// return an identity tensor
133 : static TensorGeneric<n,n> identity();
134 : /// return the matrix inverse
135 : TensorGeneric inverse()const;
136 : /// return the transpose matrix
137 : TensorGeneric<m,n> transpose()const;
138 : /// matrix-matrix multiplication
139 : template<unsigned n_,unsigned m_,unsigned l_>
140 : friend TensorGeneric<n_,l_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&);
141 : /// matrix-vector multiplication
142 : template<unsigned n_,unsigned m_>
143 : friend VectorGeneric<n_> matmul(const TensorGeneric<n_,m_>&,const VectorGeneric<m_>&);
144 : /// vector-matrix multiplication
145 : template<unsigned n_,unsigned m_>
146 : friend VectorGeneric<n_> matmul(const VectorGeneric<m_>&,const TensorGeneric<m_,n_>&);
147 : /// vector-vector multiplication (maps to dotProduct)
148 : template<unsigned n_>
149 : friend double matmul(const VectorGeneric<n_>&,const VectorGeneric<n_>&);
150 : /// matrix-matrix-matrix multiplication
151 : template<unsigned n_,unsigned m_,unsigned l_,unsigned i_>
152 : friend TensorGeneric<n_,i_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&,const TensorGeneric<l_,i_>&);
153 : /// matrix-matrix-vector multiplication
154 : template<unsigned n_,unsigned m_,unsigned l_>
155 : friend VectorGeneric<n_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&,const VectorGeneric<l_>&);
156 : /// vector-matrix-matrix multiplication
157 : template<unsigned n_,unsigned m_,unsigned l_>
158 : friend VectorGeneric<l_> matmul(const VectorGeneric<n_>&,const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&);
159 : /// vector-matrix-vector multiplication
160 : template<unsigned n_,unsigned m_>
161 : friend double matmul(const VectorGeneric<n_>&,const TensorGeneric<n_,m_>&,const VectorGeneric<m_>&);
162 : /// returns the determinant of a tensor
163 : friend double determinant(const TensorGeneric<3,3>&);
164 : /// returns the inverse of a tensor (same as inverse())
165 : friend TensorGeneric<3,3> inverse(const TensorGeneric<3,3>&);
166 : /// returns the transpose of a tensor (same as transpose())
167 : template<unsigned n_,unsigned m_>
168 : friend TensorGeneric<n_,m_> transpose(const TensorGeneric<m_,n_>&);
169 : /// returns the transpose of a tensor (same as TensorGeneric(const VectorGeneric&,const VectorGeneric&))
170 : template<unsigned n_,unsigned m_>
171 : friend TensorGeneric<n_,m_> extProduct(const VectorGeneric<n>&,const VectorGeneric<m>&);
172 4257443 : friend TensorGeneric<3,3> dcrossDv1(const VectorGeneric<3>&,const VectorGeneric<3>&);
173 4835099 : friend TensorGeneric<3,3> dcrossDv2(const VectorGeneric<3>&,const VectorGeneric<3>&);
174 144414 : friend TensorGeneric<3,3> VcrossTensor(const VectorGeneric<3>&,const TensorGeneric<3,3>&);
175 48138 : friend TensorGeneric<3,3> VcrossTensor(const TensorGeneric<3,3>&,const VectorGeneric<3>&);
176 : /// Derivative of a normalized vector
177 48138 : friend TensorGeneric<3,3> deriNorm(const VectorGeneric<3>&,const TensorGeneric<3,3>&);
178 : /// << operator.
179 : /// Allows printing tensor `t` with `std::cout<<t;`
180 : template<unsigned n_,unsigned m_>
181 : friend std::ostream & operator<<(std::ostream &os, const TensorGeneric<n_,m_>&);
182 : };
183 :
184 : template<unsigned n,unsigned m>
185 39660503 : TensorGeneric<n,m>::TensorGeneric() {
186 39660503 : LoopUnroller<n*m>::_zero(d);
187 39660503 : }
188 :
189 : template<unsigned n,unsigned m>
190 129874959 : TensorGeneric<n,m>::TensorGeneric(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2) {
191 519499836 : for(unsigned i=0; i<n; i++)for(unsigned j=0; j<m; j++)d[i*m+j]=v1[i]*v2[j];
192 129874959 : }
193 :
194 : template<>
195 : inline
196 : TensorGeneric<2,2>::TensorGeneric(double d00,double d01,double d10,double d11) {
197 : d[0]=d00;
198 : d[1]=d01;
199 : d[2]=d10;
200 : d[3]=d11;
201 : }
202 :
203 : template<>
204 : inline
205 9156237 : TensorGeneric<3,3>::TensorGeneric(double d00,double d01,double d02,double d10,double d11,double d12,double d20,double d21,double d22) {
206 9156237 : d[0]=d00;
207 9156237 : d[1]=d01;
208 9156237 : d[2]=d02;
209 9156237 : d[3]=d10;
210 9156237 : d[4]=d11;
211 9156237 : d[5]=d12;
212 9156237 : d[6]=d20;
213 9156237 : d[7]=d21;
214 9156237 : d[8]=d22;
215 : }
216 :
217 : template<unsigned n,unsigned m>
218 355191 : void TensorGeneric<n,m>::zero() {
219 355191 : LoopUnroller<n*m>::_zero(d);
220 355191 : }
221 :
222 : template<unsigned n,unsigned m>
223 178727603 : double & TensorGeneric<n,m>::operator() (unsigned i,unsigned j) {
224 : #ifdef _GLIBCXX_DEBUG
225 : plumed_assert(i<n && j<m);
226 : #endif
227 178727603 : return d[m*i+j];
228 : }
229 :
230 : template<unsigned n,unsigned m>
231 4454750629 : const double & TensorGeneric<n,m>::operator() (unsigned i,unsigned j)const {
232 : #ifdef _GLIBCXX_DEBUG
233 : plumed_assert(i<n && j<m);
234 : #endif
235 4454750629 : return d[m*i+j];
236 : }
237 :
238 : template<unsigned n,unsigned m>
239 38149660 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator +=(const TensorGeneric<n,m>& b) {
240 38149660 : LoopUnroller<n*m>::_add(d,b.d);
241 38149660 : return *this;
242 : }
243 :
244 : template<unsigned n,unsigned m>
245 24605920 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator -=(const TensorGeneric<n,m>& b) {
246 24605920 : LoopUnroller<n*m>::_sub(d,b.d);
247 24605920 : return *this;
248 : }
249 :
250 : template<unsigned n,unsigned m>
251 122591275 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator *=(double s) {
252 122591275 : LoopUnroller<n*m>::_mul(d,s);
253 122591275 : return *this;
254 : }
255 :
256 : template<unsigned n,unsigned m>
257 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator /=(double s) {
258 : LoopUnroller<n*m>::_mul(d,1.0/s);
259 : return *this;
260 : }
261 :
262 : template<unsigned n,unsigned m>
263 19111 : TensorGeneric<n,m> TensorGeneric<n,m>::operator+()const {
264 19111 : return *this;
265 : }
266 :
267 : template<unsigned n,unsigned m>
268 61898 : TensorGeneric<n,m> TensorGeneric<n,m>::operator-()const {
269 61898 : TensorGeneric<n,m> r;
270 61898 : LoopUnroller<n*m>::_neg(r.d,d);
271 61898 : return r;
272 : }
273 :
274 : template<unsigned n,unsigned m>
275 29160 : TensorGeneric<n,m>& TensorGeneric<n,m>::setCol(unsigned j,const VectorGeneric<n> & c) {
276 116640 : for(unsigned i=0; i<n; ++i) (*this)(i,j)=c(i);
277 29160 : return *this;
278 : }
279 :
280 : template<unsigned n,unsigned m>
281 2856339 : TensorGeneric<n,m>& TensorGeneric<n,m>::setRow(unsigned i,const VectorGeneric<m> & r) {
282 11425356 : for(unsigned j=0; j<m; ++j) (*this)(i,j)=r(j);
283 2856339 : return *this;
284 : }
285 :
286 : template<unsigned n,unsigned m>
287 24300 : VectorGeneric<n> TensorGeneric<n,m>::getCol(unsigned j)const {
288 24300 : VectorGeneric<n> v;
289 97200 : for(unsigned i=0; i<n; ++i) v(i)=(*this)(i,j);
290 24300 : return v;
291 : }
292 :
293 : template<unsigned n,unsigned m>
294 3171426 : VectorGeneric<m> TensorGeneric<n,m>::getRow(unsigned i)const {
295 3171426 : VectorGeneric<m> v;
296 12685704 : for(unsigned j=0; j<m; ++j) v(j)=(*this)(i,j);
297 3171426 : return v;
298 : }
299 :
300 : template<unsigned n,unsigned m>
301 2194702 : TensorGeneric<n,m> operator+(const TensorGeneric<n,m>&t1,const TensorGeneric<n,m>&t2) {
302 2194702 : TensorGeneric<n,m> t(t1);
303 2194702 : t+=t2;
304 2194702 : return t;
305 : }
306 :
307 : template<unsigned n,unsigned m>
308 1236908 : TensorGeneric<n,m> operator-(const TensorGeneric<n,m>&t1,const TensorGeneric<n,m>&t2) {
309 1236908 : TensorGeneric<n,m> t(t1);
310 1236908 : t-=t2;
311 1236908 : return t;
312 : }
313 :
314 : template<unsigned n,unsigned m>
315 122586248 : TensorGeneric<n,m> operator*(const TensorGeneric<n,m>&t1,double s) {
316 122586248 : TensorGeneric<n,m> t(t1);
317 122586248 : t*=s;
318 122586248 : return t;
319 : }
320 :
321 : template<unsigned n,unsigned m>
322 99559883 : TensorGeneric<n,m> operator*(double s,const TensorGeneric<n,m>&t1) {
323 99559883 : return t1*s;
324 : }
325 :
326 : template<unsigned n,unsigned m>
327 14 : TensorGeneric<n,m> operator/(const TensorGeneric<n,m>&t1,double s) {
328 14 : return t1*(1.0/s);
329 : }
330 :
331 : template<>
332 : inline
333 129758 : double TensorGeneric<3,3>::determinant()const {
334 : return
335 129758 : d[0]*d[4]*d[8]
336 129758 : + d[1]*d[5]*d[6]
337 129758 : + d[2]*d[3]*d[7]
338 129758 : - d[0]*d[5]*d[7]
339 129758 : - d[1]*d[3]*d[8]
340 129758 : - d[2]*d[4]*d[6];
341 : }
342 :
343 : template<unsigned n,unsigned m>
344 : inline
345 1748638 : TensorGeneric<n,n> TensorGeneric<n,m>::identity() {
346 1748638 : TensorGeneric<n,n> t;
347 6994552 : for(unsigned i=0; i<n; i++) t(i,i)=1.0;
348 1748638 : return t;
349 : }
350 :
351 : template<unsigned n,unsigned m>
352 2548408 : TensorGeneric<m,n> TensorGeneric<n,m>::transpose()const {
353 2548408 : TensorGeneric<m,n> t;
354 10193632 : for(unsigned i=0; i<m; i++)for(unsigned j=0; j<n; j++) t(i,j)=(*this)(j,i);
355 2548408 : return t;
356 : }
357 :
358 : template<>
359 : inline
360 81911 : TensorGeneric<3,3> TensorGeneric<3,3>::inverse()const {
361 81911 : TensorGeneric t;
362 81911 : double invdet=1.0/determinant();
363 1064843 : for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++)
364 1474398 : t(j,i)=invdet*( (*this)((i+1)%3,(j+1)%3)*(*this)((i+2)%3,(j+2)%3)
365 737199 : -(*this)((i+1)%3,(j+2)%3)*(*this)((i+2)%3,(j+1)%3));
366 81911 : return t;
367 : }
368 :
369 : template<unsigned n,unsigned m,unsigned l>
370 4484906 : TensorGeneric<n,l> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b) {
371 4484906 : TensorGeneric<n,l> t;
372 139032086 : for(unsigned i=0; i<n; i++) for(unsigned j=0; j<l; j++) for(unsigned k=0; k<m; k++) {
373 121092462 : t(i,j)+=a(i,k)*b(k,j);
374 : }
375 4484906 : return t;
376 : }
377 :
378 : template<unsigned n,unsigned m>
379 27216235 : VectorGeneric<n> matmul(const TensorGeneric<n,m>&a,const VectorGeneric<m>&b) {
380 27216235 : VectorGeneric<n> t;
381 108864940 : for(unsigned i=0; i<n; i++) for(unsigned j=0; j<m; j++) t(i)+=a(i,j)*b(j);
382 27216235 : return t;
383 : }
384 :
385 : template<unsigned n,unsigned m>
386 19518494 : VectorGeneric<n> matmul(const VectorGeneric<m>&a,const TensorGeneric<m,n>&b) {
387 19518494 : VectorGeneric<n> t;
388 78073976 : for(unsigned i=0; i<n; i++) for(unsigned j=0; j<m; j++) t(i)+=a(j)*b(j,i);
389 19518494 : return t;
390 : }
391 :
392 : template<unsigned n_>
393 : double matmul(const VectorGeneric<n_>&a,const VectorGeneric<n_>&b) {
394 : return dotProduct(a,b);
395 : }
396 :
397 : template<unsigned n,unsigned m,unsigned l,unsigned i>
398 : TensorGeneric<n,i> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b,const TensorGeneric<l,i>&c) {
399 : return matmul(matmul(a,b),c);
400 : }
401 :
402 : template<unsigned n,unsigned m,unsigned l>
403 179928 : VectorGeneric<n> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b,const VectorGeneric<l>&c) {
404 179928 : return matmul(matmul(a,b),c);
405 : }
406 :
407 : template<unsigned n,unsigned m,unsigned l>
408 : VectorGeneric<l> matmul(const VectorGeneric<n>&a,const TensorGeneric<n,m>&b,const TensorGeneric<m,l>&c) {
409 : return matmul(matmul(a,b),c);
410 : }
411 :
412 : template<unsigned n,unsigned m>
413 : double matmul(const VectorGeneric<n>&a,const TensorGeneric<n,m>&b,const VectorGeneric<m>&c) {
414 : return matmul(matmul(a,b),c);
415 : }
416 :
417 : inline
418 : double determinant(const TensorGeneric<3,3>&t) {
419 17 : return t.determinant();
420 : }
421 :
422 : inline
423 : TensorGeneric<3,3> inverse(const TensorGeneric<3,3>&t) {
424 41189 : return t.inverse();
425 : }
426 :
427 : template<unsigned n,unsigned m>
428 375783 : TensorGeneric<n,m> transpose(const TensorGeneric<m,n>&t) {
429 375783 : return t.transpose();
430 : }
431 :
432 : template<unsigned n,unsigned m>
433 1235235 : TensorGeneric<n,m> extProduct(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2) {
434 1235235 : return TensorGeneric<n,m>(v1,v2);
435 : }
436 :
437 : inline
438 : TensorGeneric<3,3> dcrossDv1(const VectorGeneric<3>&v1,const VectorGeneric<3>&v2) {
439 : (void) v1; // this is to avoid warnings. still the syntax of this function is a bit dummy...
440 : return TensorGeneric<3,3>(
441 : 0.0, v2[2],-v2[1],
442 4257443 : -v2[2], 0.0, v2[0],
443 12772329 : v2[1],-v2[0], 0.0);
444 : }
445 :
446 : inline
447 : TensorGeneric<3,3> dcrossDv2(const VectorGeneric<3>&v1,const VectorGeneric<3>&v2) {
448 : (void) v2; // this is to avoid warnings. still the syntax of this function is a bit dummy...
449 : return TensorGeneric<3,3>(
450 : 0.0,-v1[2],v1[1],
451 4835099 : v1[2],0.0,-v1[0],
452 14505297 : -v1[1],v1[0],0.0);
453 : }
454 :
455 : template<unsigned n,unsigned m>
456 : std::ostream & operator<<(std::ostream &os, const TensorGeneric<n,m>& t) {
457 : for(unsigned i=0; i<n; i++)for(unsigned j=0; j<m; j++) {
458 : if(i!=(n-1) || j!=(m-1)) os<<t(i,j)<<" ";
459 : else os<<t(i,j);
460 : }
461 : return os;
462 : }
463 :
464 : /// \ingroup TOOLBOX
465 : typedef TensorGeneric<2,2> Tensor2d;
466 : /// \ingroup TOOLBOX
467 : typedef TensorGeneric<3,3> Tensor3d;
468 : /// \ingroup TOOLBOX
469 : typedef TensorGeneric<4,4> Tensor4d;
470 : /// \ingroup TOOLBOX
471 : typedef Tensor3d Tensor;
472 :
473 : inline
474 : TensorGeneric<3,3> VcrossTensor(const VectorGeneric<3>&v1,const TensorGeneric<3,3>&v2) {
475 :
476 144414 : TensorGeneric<3,3> t;
477 1010898 : for(unsigned i=0; i<3; i++) {
478 433242 : t.setRow(i,matmul(dcrossDv2(v1,v1),v2.getRow(i)));
479 : }
480 144414 : return t;
481 : }
482 :
483 : inline
484 : TensorGeneric<3,3> VcrossTensor(const TensorGeneric<3,3>&v2,const VectorGeneric<3>&v1) {
485 48138 : TensorGeneric<3,3> t;
486 336966 : for(unsigned i=0; i<3; i++) {
487 144414 : t.setRow(i,-matmul(dcrossDv2(v1,v1),v2.getRow(i)));
488 : }
489 48138 : return t;
490 : }
491 :
492 :
493 : inline
494 : TensorGeneric<3,3> deriNorm(const VectorGeneric<3>&v1,const TensorGeneric<3,3>&v2) {
495 : // delta(v) = delta(v1/v1.norm) = 1/v1.norm*(delta(v1) - (v.delta(v1))cross v;
496 48138 : double over_norm = 1./v1.modulo();
497 48138 : return over_norm*(v2 - over_norm*over_norm*(extProduct(matmul(v2,v1),v1)));
498 : }
499 :
500 :
501 :
502 :
503 :
504 : }
505 :
506 : #endif
507 :
|
