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_Tensor_h
23 : #define __PLUMED_tools_Tensor_h
24 :
25 : #include "MatrixSquareBracketsAccess.h"
26 : #include "Vector.h"
27 : #include "LoopUnroller.h"
28 : #include "Exception.h"
29 :
30 : #include <array>
31 :
32 : namespace PLMD {
33 :
34 : /// Small class to contain local utilities.
35 : /// Should not be used outside of the TensorGeneric class.
36 : class TensorGenericAux {
37 : public:
38 : // local redefinition, just to avoid including lapack.h here
39 : static void local_dsyevr(const char *jobz, const char *range, const char *uplo, int *n,
40 : double *a, int *lda, double *vl, double *vu, int *
41 : il, int *iu, double *abstol, int *m, double *w,
42 : double *z__, int *ldz, int *isuppz, double *work,
43 : int *lwork, int *iwork, int *liwork, int *info);
44 : };
45 :
46 : /**
47 : \ingroup TOOLBOX
48 : Class implementing fixed size matrices of doubles
49 :
50 : \tparam n The number rows
51 : \tparam m The number columns
52 :
53 : This class implements a matrix of doubles with size fixed at
54 : compile time. It is useful for small fixed size objects (e.g.
55 : 3x3 tensors) as it does not waste space to store the vector size.
56 : Moreover, as the compiler knows the size, it can be completely
57 : opimized inline.
58 : Most of the loops are explicitly unrolled using PLMD::LoopUnroller class
59 : Matrix elements are initialized to zero by default. Notice that
60 : this means that constructor is a bit slow. This point might change
61 : in future if we find performance issues.
62 : It takes advantage of MatrixSquareBracketsAccess to provide both
63 : () and [] syntax for access.
64 : Several functions are declared as friends even if not necessary so as to
65 : properly appear in Doxygen documentation.
66 :
67 : Aliases are defined to simplify common declarations (Tensor, Tensor2d, Tensor3d, Tensor4d).
68 : Also notice that some operations are only available for 3x3 tensors.
69 :
70 : Example of usage
71 : \verbatim
72 :
73 : #include "Tensor.h"
74 :
75 : using namespace PLMD;
76 :
77 : int main(){
78 : Tensor a;
79 : TensorGeneric<3,2> b;
80 : TensorGeneric<3,2> c=matmul(a,b);
81 : return 0;
82 : }
83 :
84 : \endverbatim
85 : */
86 : template <unsigned n,unsigned m>
87 : class TensorGeneric:
88 : public MatrixSquareBracketsAccess<TensorGeneric<n,m>,double> {
89 : std::array<double,n*m> d;
90 : /// Auxiliary private function for constructor
91 : void auxiliaryConstructor();
92 : /// Auxiliary private function for constructor
93 : template<typename... Args>
94 : void auxiliaryConstructor(double first,Args... arg);
95 : public:
96 : /// Constructor accepting n*m double parameters.
97 : /// Can be used as Tensor<2,2>(1.0,2.0,3.0,4.0)
98 : /// In case a wrong number of parameters is given, a static assertion will fail.
99 : template<typename... Args>
100 : TensorGeneric(double first,Args... arg);
101 : /// initialize the tensor to zero
102 : TensorGeneric();
103 : /// initialize a tensor as an external product of two Vector
104 : TensorGeneric(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2);
105 : /// set it to zero
106 : void zero();
107 : /// access element
108 : double & operator() (unsigned i,unsigned j);
109 : /// access element
110 : const double & operator() (unsigned i,unsigned j)const;
111 : /// increment
112 : TensorGeneric& operator +=(const TensorGeneric<n,m>& b);
113 : /// decrement
114 : TensorGeneric& operator -=(const TensorGeneric<n,m>& b);
115 : /// multiply
116 : TensorGeneric& operator *=(double);
117 : /// divide
118 : TensorGeneric& operator /=(double);
119 : /// return +t
120 : TensorGeneric operator +()const;
121 : /// return -t
122 : TensorGeneric operator -()const;
123 : /// set j-th column
124 : TensorGeneric& setCol(unsigned j,const VectorGeneric<n> & c);
125 : /// set i-th row
126 : TensorGeneric& setRow(unsigned i,const VectorGeneric<m> & r);
127 : /// get j-th column
128 : VectorGeneric<n> getCol(unsigned j)const;
129 : /// get i-th row
130 : VectorGeneric<m> getRow(unsigned i)const;
131 : /// return t1+t2
132 : template<unsigned n_,unsigned m_>
133 : friend TensorGeneric<n_,m_> operator+(const TensorGeneric<n_,m_>&,const TensorGeneric<n_,m_>&);
134 : /// return t1+t2
135 : template<unsigned n_,unsigned m_>
136 : friend TensorGeneric<n_,m_> operator-(const TensorGeneric<n_,m_>&,const TensorGeneric<n_,m_>&);
137 : /// scale the tensor by a factor s
138 : template<unsigned n_,unsigned m_>
139 : friend TensorGeneric<n_,m_> operator*(double,const TensorGeneric<n_,m_>&);
140 : /// scale the tensor by a factor s
141 : template<unsigned n_,unsigned m_>
142 : friend TensorGeneric<n_,m_> operator*(const TensorGeneric<n_,m_>&,double s);
143 : /// scale the tensor by a factor 1/s
144 : template<unsigned n_,unsigned m_>
145 : friend TensorGeneric<n_,m_> operator/(const TensorGeneric<n_,m_>&,double s);
146 : /// returns the determinant
147 : double determinant()const;
148 : /// return an identity tensor
149 : static TensorGeneric<n,n> identity();
150 : /// return the matrix inverse
151 : TensorGeneric inverse()const;
152 : /// return the transpose matrix
153 : TensorGeneric<m,n> transpose()const;
154 : /// matrix-matrix multiplication
155 : template<unsigned n_,unsigned m_,unsigned l_>
156 : friend TensorGeneric<n_,l_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&);
157 : /// matrix-vector multiplication
158 : template<unsigned n_,unsigned m_>
159 : friend VectorGeneric<n_> matmul(const TensorGeneric<n_,m_>&,const VectorGeneric<m_>&);
160 : /// vector-matrix multiplication
161 : template<unsigned n_,unsigned m_>
162 : friend VectorGeneric<n_> matmul(const VectorGeneric<m_>&,const TensorGeneric<m_,n_>&);
163 : /// vector-vector multiplication (maps to dotProduct)
164 : template<unsigned n_>
165 : friend double matmul(const VectorGeneric<n_>&,const VectorGeneric<n_>&);
166 : /// matrix-matrix-matrix multiplication
167 : template<unsigned n_,unsigned m_,unsigned l_,unsigned i_>
168 : friend TensorGeneric<n_,i_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&,const TensorGeneric<l_,i_>&);
169 : /// matrix-matrix-vector multiplication
170 : template<unsigned n_,unsigned m_,unsigned l_>
171 : friend VectorGeneric<n_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&,const VectorGeneric<l_>&);
172 : /// vector-matrix-matrix multiplication
173 : template<unsigned n_,unsigned m_,unsigned l_>
174 : friend VectorGeneric<l_> matmul(const VectorGeneric<n_>&,const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&);
175 : /// vector-matrix-vector multiplication
176 : template<unsigned n_,unsigned m_>
177 : friend double matmul(const VectorGeneric<n_>&,const TensorGeneric<n_,m_>&,const VectorGeneric<m_>&);
178 : /// returns the determinant of a tensor
179 : friend double determinant(const TensorGeneric<3,3>&);
180 : /// returns the inverse of a tensor (same as inverse())
181 : friend TensorGeneric<3,3> inverse(const TensorGeneric<3,3>&);
182 : /// returns the transpose of a tensor (same as transpose())
183 : template<unsigned n_,unsigned m_>
184 : friend TensorGeneric<n_,m_> transpose(const TensorGeneric<m_,n_>&);
185 : /// returns the transpose of a tensor (same as TensorGeneric(const VectorGeneric&,const VectorGeneric&))
186 : template<unsigned n_,unsigned m_>
187 : friend TensorGeneric<n_,m_> extProduct(const VectorGeneric<n>&,const VectorGeneric<m>&);
188 : friend TensorGeneric<3,3> dcrossDv1(const VectorGeneric<3>&,const VectorGeneric<3>&);
189 : friend TensorGeneric<3,3> dcrossDv2(const VectorGeneric<3>&,const VectorGeneric<3>&);
190 : friend TensorGeneric<3,3> VcrossTensor(const VectorGeneric<3>&,const TensorGeneric<3,3>&);
191 : friend TensorGeneric<3,3> VcrossTensor(const TensorGeneric<3,3>&,const VectorGeneric<3>&);
192 : /// Derivative of a normalized vector
193 : friend TensorGeneric<3,3> deriNorm(const VectorGeneric<3>&,const TensorGeneric<3,3>&);
194 : /// << operator.
195 : /// Allows printing tensor `t` with `std::cout<<t;`
196 : template<unsigned n_,unsigned m_>
197 : friend std::ostream & operator<<(std::ostream &os, const TensorGeneric<n_,m_>&);
198 : /// Diagonalize tensor.
199 : /// Syntax is the same as Matrix::diagMat.
200 : /// In addition, it is possible to call if with m_ smaller than n_. In this case,
201 : /// only the first (smaller) m_ eigenvalues and eigenvectors are retrieved.
202 : /// If case lapack fails (info!=0) it throws an exception.
203 : /// Notice that tensor is assumed to be symmetric!!!
204 : template<unsigned n_,unsigned m_>
205 : friend void diagMatSym(const TensorGeneric<n_,n_>&,VectorGeneric<m_>&evals,TensorGeneric<m_,n_>&evec);
206 : };
207 :
208 : template <unsigned n,unsigned m>
209 : void TensorGeneric<n,m>::auxiliaryConstructor()
210 : {}
211 :
212 : template <unsigned n,unsigned m>
213 : template<typename... Args>
214 103821930 : void TensorGeneric<n,m>::auxiliaryConstructor(double first,Args... arg) {
215 103821930 : d[n*m-(sizeof...(Args))-1]=first;
216 92286160 : auxiliaryConstructor(arg...);
217 103821930 : }
218 :
219 : template <unsigned n,unsigned m>
220 : template<typename... Args>
221 11535770 : TensorGeneric<n,m>::TensorGeneric(double first,Args... arg) {
222 : static_assert((sizeof...(Args))+1==n*m,"you are trying to initialize a Tensor with the wrong number of arguments");
223 11535770 : auxiliaryConstructor(first,arg...);
224 11535770 : }
225 :
226 : template<unsigned n,unsigned m>
227 60486621 : TensorGeneric<n,m>::TensorGeneric() {
228 60486621 : LoopUnroller<n*m>::_zero(d.data());
229 60486621 : }
230 :
231 : template<unsigned n,unsigned m>
232 105255401 : TensorGeneric<n,m>::TensorGeneric(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2) {
233 421021604 : for(unsigned i=0; i<n; i++)
234 1263064812 : for(unsigned j=0; j<m; j++) {
235 947298609 : d[i*m+j]=v1[i]*v2[j];
236 : }
237 105255401 : }
238 :
239 : template<unsigned n,unsigned m>
240 957800 : void TensorGeneric<n,m>::zero() {
241 957800 : LoopUnroller<n*m>::_zero(d.data());
242 957800 : }
243 :
244 : template<unsigned n,unsigned m>
245 353622406 : double & TensorGeneric<n,m>::operator() (unsigned i,unsigned j) {
246 : #ifdef _GLIBCXX_DEBUG
247 : // index i is implicitly checked by the std::array class
248 : plumed_assert(j<m);
249 : #endif
250 353622406 : return d[m*i+j];
251 : }
252 :
253 : template<unsigned n,unsigned m>
254 5063732474 : const double & TensorGeneric<n,m>::operator() (unsigned i,unsigned j)const {
255 : #ifdef _GLIBCXX_DEBUG
256 : // index i is implicitly checked by the std::array class
257 : plumed_assert(j<m);
258 : #endif
259 5063732474 : return d[m*i+j];
260 : }
261 :
262 : template<unsigned n,unsigned m>
263 63669845 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator +=(const TensorGeneric<n,m>& b) {
264 63669845 : LoopUnroller<n*m>::_add(d.data(),b.d.data());
265 63669845 : return *this;
266 : }
267 :
268 : template<unsigned n,unsigned m>
269 49313975 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator -=(const TensorGeneric<n,m>& b) {
270 49313975 : LoopUnroller<n*m>::_sub(d.data(),b.d.data());
271 49313975 : return *this;
272 : }
273 :
274 : template<unsigned n,unsigned m>
275 85306132 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator *=(double s) {
276 85306132 : LoopUnroller<n*m>::_mul(d.data(),s);
277 85306132 : return *this;
278 : }
279 :
280 : template<unsigned n,unsigned m>
281 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator /=(double s) {
282 : LoopUnroller<n*m>::_mul(d.data(),1.0/s);
283 : return *this;
284 : }
285 :
286 : template<unsigned n,unsigned m>
287 177723 : TensorGeneric<n,m> TensorGeneric<n,m>::operator+()const {
288 177723 : return *this;
289 : }
290 :
291 : template<unsigned n,unsigned m>
292 250729 : TensorGeneric<n,m> TensorGeneric<n,m>::operator-()const {
293 250729 : TensorGeneric<n,m> r;
294 250729 : LoopUnroller<n*m>::_neg(r.d.data(),d.data());
295 250729 : return r;
296 : }
297 :
298 : template<unsigned n,unsigned m>
299 79542 : TensorGeneric<n,m>& TensorGeneric<n,m>::setCol(unsigned j,const VectorGeneric<n> & c) {
300 318168 : for(unsigned i=0; i<n; ++i) {
301 238626 : (*this)(i,j)=c(i);
302 : }
303 79542 : return *this;
304 : }
305 :
306 : template<unsigned n,unsigned m>
307 3075777 : TensorGeneric<n,m>& TensorGeneric<n,m>::setRow(unsigned i,const VectorGeneric<m> & r) {
308 12303108 : for(unsigned j=0; j<m; ++j) {
309 9227331 : (*this)(i,j)=r(j);
310 : }
311 3075777 : return *this;
312 : }
313 :
314 : template<unsigned n,unsigned m>
315 12312882 : VectorGeneric<n> TensorGeneric<n,m>::getCol(unsigned j)const {
316 12312882 : VectorGeneric<n> v;
317 61455240 : for(unsigned i=0; i<n; ++i) {
318 49142358 : v(i)=(*this)(i,j);
319 : }
320 12312882 : return v;
321 : }
322 :
323 : template<unsigned n,unsigned m>
324 3536679 : VectorGeneric<m> TensorGeneric<n,m>::getRow(unsigned i)const {
325 3536679 : VectorGeneric<m> v;
326 14146716 : for(unsigned j=0; j<m; ++j) {
327 10610037 : v(j)=(*this)(i,j);
328 : }
329 3536679 : return v;
330 : }
331 :
332 : template<unsigned n,unsigned m>
333 2932917 : TensorGeneric<n,m> operator+(const TensorGeneric<n,m>&t1,const TensorGeneric<n,m>&t2) {
334 2932917 : TensorGeneric<n,m> t(t1);
335 2932917 : t+=t2;
336 2932917 : return t;
337 : }
338 :
339 : template<unsigned n,unsigned m>
340 1822586 : TensorGeneric<n,m> operator-(const TensorGeneric<n,m>&t1,const TensorGeneric<n,m>&t2) {
341 1822586 : TensorGeneric<n,m> t(t1);
342 1822586 : t-=t2;
343 1822586 : return t;
344 : }
345 :
346 : template<unsigned n,unsigned m>
347 85301107 : TensorGeneric<n,m> operator*(const TensorGeneric<n,m>&t1,double s) {
348 85301107 : TensorGeneric<n,m> t(t1);
349 85301107 : t*=s;
350 85301107 : return t;
351 : }
352 :
353 : template<unsigned n,unsigned m>
354 59160541 : TensorGeneric<n,m> operator*(double s,const TensorGeneric<n,m>&t1) {
355 59160541 : return t1*s;
356 : }
357 :
358 : template<unsigned n,unsigned m>
359 37375 : TensorGeneric<n,m> operator/(const TensorGeneric<n,m>&t1,double s) {
360 37375 : return t1*(1.0/s);
361 : }
362 :
363 : template<>
364 : inline
365 290868 : double TensorGeneric<3,3>::determinant()const {
366 : return
367 290868 : d[0]*d[4]*d[8]
368 290868 : + d[1]*d[5]*d[6]
369 290868 : + d[2]*d[3]*d[7]
370 290868 : - d[0]*d[5]*d[7]
371 290868 : - d[1]*d[3]*d[8]
372 290868 : - d[2]*d[4]*d[6];
373 : }
374 :
375 : template<unsigned n,unsigned m>
376 : inline
377 13276328 : TensorGeneric<n,n> TensorGeneric<n,m>::identity() {
378 13276328 : TensorGeneric<n,n> t;
379 53105312 : for(unsigned i=0; i<n; i++) {
380 39828984 : t(i,i)=1.0;
381 : }
382 13276328 : return t;
383 : }
384 :
385 : template<unsigned n,unsigned m>
386 4906432 : TensorGeneric<m,n> TensorGeneric<n,m>::transpose()const {
387 4906432 : TensorGeneric<m,n> t;
388 19625728 : for(unsigned i=0; i<m; i++)
389 58877184 : for(unsigned j=0; j<n; j++) {
390 44157888 : t(i,j)=(*this)(j,i);
391 : }
392 4906432 : return t;
393 : }
394 :
395 : template<>
396 : inline
397 182815 : TensorGeneric<3,3> TensorGeneric<3,3>::inverse()const {
398 182815 : TensorGeneric t;
399 182815 : double invdet=1.0/determinant();
400 731260 : for(unsigned i=0; i<3; i++)
401 2193780 : for(unsigned j=0; j<3; j++)
402 1645335 : t(j,i)=invdet*( (*this)((i+1)%3,(j+1)%3)*(*this)((i+2)%3,(j+2)%3)
403 1645335 : -(*this)((i+1)%3,(j+2)%3)*(*this)((i+2)%3,(j+1)%3));
404 182815 : return t;
405 : }
406 :
407 : template<unsigned n,unsigned m,unsigned l>
408 5620586 : TensorGeneric<n,l> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b) {
409 5620586 : TensorGeneric<n,l> t;
410 22482344 : for(unsigned i=0; i<n; i++)
411 67447032 : for(unsigned j=0; j<l; j++)
412 202341096 : for(unsigned k=0; k<m; k++) {
413 151755822 : t(i,j)+=a(i,k)*b(k,j);
414 : }
415 5620586 : return t;
416 : }
417 :
418 : template<unsigned n,unsigned m>
419 61206072 : VectorGeneric<n> matmul(const TensorGeneric<n,m>&a,const VectorGeneric<m>&b) {
420 61206072 : VectorGeneric<n> t;
421 244824288 : for(unsigned i=0; i<n; i++)
422 734472864 : for(unsigned j=0; j<m; j++) {
423 550854648 : t(i)+=a(i,j)*b(j);
424 : }
425 61206072 : return t;
426 : }
427 :
428 : template<unsigned n,unsigned m>
429 107713661 : VectorGeneric<n> matmul(const VectorGeneric<m>&a,const TensorGeneric<m,n>&b) {
430 107713661 : VectorGeneric<n> t;
431 430854644 : for(unsigned i=0; i<n; i++)
432 1293718614 : for(unsigned j=0; j<m; j++) {
433 970577631 : t(i)+=a(j)*b(j,i);
434 : }
435 107713661 : return t;
436 : }
437 :
438 : template<unsigned n_>
439 : double matmul(const VectorGeneric<n_>&a,const VectorGeneric<n_>&b) {
440 : return dotProduct(a,b);
441 : }
442 :
443 : template<unsigned n,unsigned m,unsigned l,unsigned i>
444 : TensorGeneric<n,i> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b,const TensorGeneric<l,i>&c) {
445 : return matmul(matmul(a,b),c);
446 : }
447 :
448 : template<unsigned n,unsigned m,unsigned l>
449 364320 : VectorGeneric<n> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b,const VectorGeneric<l>&c) {
450 364320 : return matmul(matmul(a,b),c);
451 : }
452 :
453 : template<unsigned n,unsigned m,unsigned l>
454 : VectorGeneric<l> matmul(const VectorGeneric<n>&a,const TensorGeneric<n,m>&b,const TensorGeneric<m,l>&c) {
455 : return matmul(matmul(a,b),c);
456 : }
457 :
458 : template<unsigned n,unsigned m>
459 : double matmul(const VectorGeneric<n>&a,const TensorGeneric<n,m>&b,const VectorGeneric<m>&c) {
460 : return matmul(matmul(a,b),c);
461 : }
462 :
463 : inline
464 : double determinant(const TensorGeneric<3,3>&t) {
465 417 : return t.determinant();
466 : }
467 :
468 : inline
469 : TensorGeneric<3,3> inverse(const TensorGeneric<3,3>&t) {
470 91608 : return t.inverse();
471 : }
472 :
473 : template<unsigned n,unsigned m>
474 870278 : TensorGeneric<n,m> transpose(const TensorGeneric<m,n>&t) {
475 870278 : return t.transpose();
476 : }
477 :
478 : template<unsigned n,unsigned m>
479 2612390 : TensorGeneric<n,m> extProduct(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2) {
480 2612390 : return TensorGeneric<n,m>(v1,v2);
481 : }
482 :
483 : inline
484 5185159 : TensorGeneric<3,3> dcrossDv1(const VectorGeneric<3>&v1,const VectorGeneric<3>&v2) {
485 : (void) v1; // this is to avoid warnings. still the syntax of this function is a bit dummy...
486 : return TensorGeneric<3,3>(
487 : 0.0, v2[2],-v2[1],
488 5185159 : -v2[2], 0.0, v2[0],
489 5185159 : v2[1],-v2[0], 0.0);
490 : }
491 :
492 : inline
493 5762815 : TensorGeneric<3,3> dcrossDv2(const VectorGeneric<3>&v1,const VectorGeneric<3>&v2) {
494 : (void) v2; // this is to avoid warnings. still the syntax of this function is a bit dummy...
495 : return TensorGeneric<3,3>(
496 : 0.0,-v1[2],v1[1],
497 5762815 : v1[2],0.0,-v1[0],
498 5762815 : -v1[1],v1[0],0.0);
499 : }
500 :
501 : template<unsigned n,unsigned m>
502 0 : std::ostream & operator<<(std::ostream &os, const TensorGeneric<n,m>& t) {
503 0 : for(unsigned i=0; i<n; i++)
504 0 : for(unsigned j=0; j<m; j++) {
505 0 : if(i!=(n-1) || j!=(m-1)) {
506 0 : os<<t(i,j)<<" ";
507 : } else {
508 0 : os<<t(i,j);
509 : }
510 : }
511 0 : return os;
512 : }
513 :
514 : /// \ingroup TOOLBOX
515 : typedef TensorGeneric<1,1> Tensor1d;
516 : /// \ingroup TOOLBOX
517 : typedef TensorGeneric<2,2> Tensor2d;
518 : /// \ingroup TOOLBOX
519 : typedef TensorGeneric<3,3> Tensor3d;
520 : /// \ingroup TOOLBOX
521 : typedef TensorGeneric<4,4> Tensor4d;
522 : /// \ingroup TOOLBOX
523 : typedef TensorGeneric<5,5> Tensor5d;
524 : /// \ingroup TOOLBOX
525 : typedef Tensor3d Tensor;
526 :
527 : inline
528 144414 : TensorGeneric<3,3> VcrossTensor(const VectorGeneric<3>&v1,const TensorGeneric<3,3>&v2) {
529 :
530 144414 : TensorGeneric<3,3> t;
531 577656 : for(unsigned i=0; i<3; i++) {
532 433242 : t.setRow(i,matmul(dcrossDv2(v1,v1),v2.getRow(i)));
533 : }
534 144414 : return t;
535 : }
536 :
537 : inline
538 48138 : TensorGeneric<3,3> VcrossTensor(const TensorGeneric<3,3>&v2,const VectorGeneric<3>&v1) {
539 48138 : TensorGeneric<3,3> t;
540 192552 : for(unsigned i=0; i<3; i++) {
541 144414 : t.setRow(i,-matmul(dcrossDv2(v1,v1),v2.getRow(i)));
542 : }
543 48138 : return t;
544 : }
545 :
546 :
547 : inline
548 48138 : TensorGeneric<3,3> deriNorm(const VectorGeneric<3>&v1,const TensorGeneric<3,3>&v2) {
549 : // delta(v) = delta(v1/v1.norm) = 1/v1.norm*(delta(v1) - (v.delta(v1))cross v;
550 48138 : double over_norm = 1./v1.modulo();
551 48138 : return over_norm*(v2 - over_norm*over_norm*(extProduct(matmul(v2,v1),v1)));
552 : }
553 :
554 : template<unsigned n,unsigned m>
555 641954 : void diagMatSym(const TensorGeneric<n,n>&mat,VectorGeneric<m>&evals,TensorGeneric<m,n>&evec) {
556 : // some guess number to make sure work is large enough.
557 : // for correctness it should be >=20. However, it is recommended to be the block size.
558 : // I put some likely exaggerated number
559 : constexpr int bs=100;
560 : // temporary data, on stack so as to avoid allocations
561 : std::array<int,10*n> iwork;
562 : std::array<double,(6+bs)*n> work;
563 : std::array<int,2*m> isup;
564 : // documentation says that "matrix is destroyed" !!!
565 641954 : auto mat_copy=mat;
566 : // documentation says this is size n (even if m<n)
567 : std::array<double,n> evals_tmp;
568 641954 : int nn=n; // dimension of matrix
569 641954 : double vl=0.0, vu=1.0; // ranges - not used
570 641954 : int one=1,mm=m; // minimum and maximum index
571 641954 : double abstol=0.0; // tolerance
572 641954 : int mout=0; // number of eigenvalues found (same as mm)
573 641954 : int info=0; // result
574 641954 : int liwork=iwork.size();
575 641954 : int lwork=work.size();
576 641954 : TensorGenericAux::local_dsyevr("V", (n==m?"A":"I"), "U", &nn, const_cast<double*>(&mat_copy[0][0]), &nn, &vl, &vu, &one, &mm,
577 641954 : &abstol, &mout, &evals_tmp[0], &evec[0][0], &nn,
578 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
579 641954 : if(info!=0)
580 0 : plumed_error()<<"Error diagonalizing matrix\n"
581 : <<"Matrix:\n"<<mat<<"\n"
582 : <<"Info: "<<info<<"\n";
583 641954 : plumed_assert(mout==m);
584 1461150 : for(unsigned i=0; i<m; i++) {
585 819196 : evals[i]=evals_tmp[i];
586 : }
587 : // This changes eigenvectors so that the first non-null element
588 : // of each of them is positive
589 : // We can do it because the phase is arbitrary, and helps making
590 : // the result reproducible
591 1461150 : for(unsigned i=0; i<m; ++i) {
592 : unsigned j=0;
593 819789 : for(j=0; j<n; j++)
594 819789 : if(evec(i,j)*evec(i,j)>1e-14) {
595 : break;
596 : }
597 819196 : if(j<n)
598 819196 : if(evec(i,j)<0.0)
599 780639 : for(j=0; j<n; j++) {
600 624427 : evec(i,j)*=-1;
601 : }
602 : }
603 641954 : }
604 :
605 : static_assert(sizeof(Tensor)==9*sizeof(double), "code cannot work if this is not satisfied");
606 :
607 : }
608 :
609 : #endif
610 :
|
