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 : {
90 : std::array<double,n*m> d;
91 : /// Auxiliary private function for constructor
92 : void auxiliaryConstructor();
93 : /// Auxiliary private function for constructor
94 : template<typename... Args>
95 : void auxiliaryConstructor(double first,Args... arg);
96 : public:
97 : /// Constructor accepting n*m double parameters.
98 : /// Can be used as Tensor<2,2>(1.0,2.0,3.0,4.0)
99 : /// In case a wrong number of parameters is given, a static assertion will fail.
100 : template<typename... Args>
101 : TensorGeneric(double first,Args... arg);
102 : /// initialize the tensor to zero
103 : TensorGeneric();
104 : /// initialize a tensor as an external product of two Vector
105 : TensorGeneric(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2);
106 : /// set it to zero
107 : void zero();
108 : /// access element
109 : double & operator() (unsigned i,unsigned j);
110 : /// access element
111 : const double & operator() (unsigned i,unsigned j)const;
112 : /// increment
113 : TensorGeneric& operator +=(const TensorGeneric<n,m>& b);
114 : /// decrement
115 : TensorGeneric& operator -=(const TensorGeneric<n,m>& b);
116 : /// multiply
117 : TensorGeneric& operator *=(double);
118 : /// divide
119 : TensorGeneric& operator /=(double);
120 : /// return +t
121 : TensorGeneric operator +()const;
122 : /// return -t
123 : TensorGeneric operator -()const;
124 : /// set j-th column
125 : TensorGeneric& setCol(unsigned j,const VectorGeneric<n> & c);
126 : /// set i-th row
127 : TensorGeneric& setRow(unsigned i,const VectorGeneric<m> & r);
128 : /// get j-th column
129 : VectorGeneric<n> getCol(unsigned j)const;
130 : /// get i-th row
131 : VectorGeneric<m> getRow(unsigned i)const;
132 : /// return t1+t2
133 : template<unsigned n_,unsigned m_>
134 : friend TensorGeneric<n_,m_> operator+(const TensorGeneric<n_,m_>&,const TensorGeneric<n_,m_>&);
135 : /// return t1+t2
136 : template<unsigned n_,unsigned m_>
137 : friend TensorGeneric<n_,m_> operator-(const TensorGeneric<n_,m_>&,const TensorGeneric<n_,m_>&);
138 : /// scale the tensor by a factor s
139 : template<unsigned n_,unsigned m_>
140 : friend TensorGeneric<n_,m_> operator*(double,const TensorGeneric<n_,m_>&);
141 : /// scale the tensor by a factor s
142 : template<unsigned n_,unsigned m_>
143 : friend TensorGeneric<n_,m_> operator*(const TensorGeneric<n_,m_>&,double s);
144 : /// scale the tensor by a factor 1/s
145 : template<unsigned n_,unsigned m_>
146 : friend TensorGeneric<n_,m_> operator/(const TensorGeneric<n_,m_>&,double s);
147 : /// returns the determinant
148 : double determinant()const;
149 : /// return an identity tensor
150 : static TensorGeneric<n,n> identity();
151 : /// return the matrix inverse
152 : TensorGeneric inverse()const;
153 : /// return the transpose matrix
154 : TensorGeneric<m,n> transpose()const;
155 : /// matrix-matrix multiplication
156 : template<unsigned n_,unsigned m_,unsigned l_>
157 : friend TensorGeneric<n_,l_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&);
158 : /// matrix-vector multiplication
159 : template<unsigned n_,unsigned m_>
160 : friend VectorGeneric<n_> matmul(const TensorGeneric<n_,m_>&,const VectorGeneric<m_>&);
161 : /// vector-matrix multiplication
162 : template<unsigned n_,unsigned m_>
163 : friend VectorGeneric<n_> matmul(const VectorGeneric<m_>&,const TensorGeneric<m_,n_>&);
164 : /// vector-vector multiplication (maps to dotProduct)
165 : template<unsigned n_>
166 : friend double matmul(const VectorGeneric<n_>&,const VectorGeneric<n_>&);
167 : /// matrix-matrix-matrix multiplication
168 : template<unsigned n_,unsigned m_,unsigned l_,unsigned i_>
169 : friend TensorGeneric<n_,i_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&,const TensorGeneric<l_,i_>&);
170 : /// matrix-matrix-vector multiplication
171 : template<unsigned n_,unsigned m_,unsigned l_>
172 : friend VectorGeneric<n_> matmul(const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&,const VectorGeneric<l_>&);
173 : /// vector-matrix-matrix multiplication
174 : template<unsigned n_,unsigned m_,unsigned l_>
175 : friend VectorGeneric<l_> matmul(const VectorGeneric<n_>&,const TensorGeneric<n_,m_>&,const TensorGeneric<m_,l_>&);
176 : /// vector-matrix-vector multiplication
177 : template<unsigned n_,unsigned m_>
178 : friend double matmul(const VectorGeneric<n_>&,const TensorGeneric<n_,m_>&,const VectorGeneric<m_>&);
179 : /// returns the determinant of a tensor
180 : friend double determinant(const TensorGeneric<3,3>&);
181 : /// returns the inverse of a tensor (same as inverse())
182 : friend TensorGeneric<3,3> inverse(const TensorGeneric<3,3>&);
183 : /// returns the transpose of a tensor (same as transpose())
184 : template<unsigned n_,unsigned m_>
185 : friend TensorGeneric<n_,m_> transpose(const TensorGeneric<m_,n_>&);
186 : /// returns the transpose of a tensor (same as TensorGeneric(const VectorGeneric&,const VectorGeneric&))
187 : template<unsigned n_,unsigned m_>
188 : friend TensorGeneric<n_,m_> extProduct(const VectorGeneric<n>&,const VectorGeneric<m>&);
189 : friend TensorGeneric<3,3> dcrossDv1(const VectorGeneric<3>&,const VectorGeneric<3>&);
190 : friend TensorGeneric<3,3> dcrossDv2(const VectorGeneric<3>&,const VectorGeneric<3>&);
191 : friend TensorGeneric<3,3> VcrossTensor(const VectorGeneric<3>&,const TensorGeneric<3,3>&);
192 : friend TensorGeneric<3,3> VcrossTensor(const TensorGeneric<3,3>&,const VectorGeneric<3>&);
193 : /// Derivative of a normalized vector
194 : friend TensorGeneric<3,3> deriNorm(const VectorGeneric<3>&,const TensorGeneric<3,3>&);
195 : /// << operator.
196 : /// Allows printing tensor `t` with `std::cout<<t;`
197 : template<unsigned n_,unsigned m_>
198 : friend std::ostream & operator<<(std::ostream &os, const TensorGeneric<n_,m_>&);
199 : /// Diagonalize tensor.
200 : /// Syntax is the same as Matrix::diagMat.
201 : /// In addition, it is possible to call if with m_ smaller than n_. In this case,
202 : /// only the first (smaller) m_ eigenvalues and eigenvectors are retrieved.
203 : /// If case lapack fails (info!=0) it throws an exception.
204 : /// Notice that tensor is assumed to be symmetric!!!
205 : template<unsigned n_,unsigned m_>
206 : friend void diagMatSym(const TensorGeneric<n_,n_>&,VectorGeneric<m_>&evals,TensorGeneric<m_,n_>&evec);
207 : };
208 :
209 : template <unsigned n,unsigned m>
210 : void TensorGeneric<n,m>::auxiliaryConstructor()
211 : {}
212 :
213 : template <unsigned n,unsigned m>
214 : template<typename... Args>
215 103767498 : void TensorGeneric<n,m>::auxiliaryConstructor(double first,Args... arg)
216 : {
217 103767498 : d[n*m-(sizeof...(Args))-1]=first;
218 92237776 : auxiliaryConstructor(arg...);
219 103767498 : }
220 :
221 : template <unsigned n,unsigned m>
222 : template<typename... Args>
223 11529722 : TensorGeneric<n,m>::TensorGeneric(double first,Args... arg)
224 : {
225 : static_assert((sizeof...(Args))+1==n*m,"you are trying to initialize a Tensor with the wrong number of arguments");
226 11529722 : auxiliaryConstructor(first,arg...);
227 11529722 : }
228 :
229 : template<unsigned n,unsigned m>
230 59806532 : TensorGeneric<n,m>::TensorGeneric() {
231 59806532 : LoopUnroller<n*m>::_zero(d.data());
232 59806532 : }
233 :
234 : template<unsigned n,unsigned m>
235 101538290 : TensorGeneric<n,m>::TensorGeneric(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2) {
236 1319997770 : for(unsigned i=0; i<n; i++)for(unsigned j=0; j<m; j++)d[i*m+j]=v1[i]*v2[j];
237 101538290 : }
238 :
239 : template<unsigned n,unsigned m>
240 932023 : void TensorGeneric<n,m>::zero() {
241 932023 : LoopUnroller<n*m>::_zero(d.data());
242 932023 : }
243 :
244 : template<unsigned n,unsigned m>
245 351794501 : 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 351794501 : return d[m*i+j];
251 : }
252 :
253 : template<unsigned n,unsigned m>
254 5054268506 : 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 5054268506 : return d[m*i+j];
260 : }
261 :
262 : template<unsigned n,unsigned m>
263 60103808 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator +=(const TensorGeneric<n,m>& b) {
264 60103808 : LoopUnroller<n*m>::_add(d.data(),b.d.data());
265 60103808 : return *this;
266 : }
267 :
268 : template<unsigned n,unsigned m>
269 49162145 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator -=(const TensorGeneric<n,m>& b) {
270 49162145 : LoopUnroller<n*m>::_sub(d.data(),b.d.data());
271 49162145 : return *this;
272 : }
273 :
274 : template<unsigned n,unsigned m>
275 81597648 : TensorGeneric<n,m>& TensorGeneric<n,m>::operator *=(double s) {
276 81597648 : LoopUnroller<n*m>::_mul(d.data(),s);
277 81597648 : 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 250728 : TensorGeneric<n,m> TensorGeneric<n,m>::operator-()const {
293 250728 : TensorGeneric<n,m> r;
294 250728 : LoopUnroller<n*m>::_neg(r.d.data(),d.data());
295 250728 : 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) (*this)(i,j)=c(i);
301 79542 : return *this;
302 : }
303 :
304 : template<unsigned n,unsigned m>
305 3075597 : TensorGeneric<n,m>& TensorGeneric<n,m>::setRow(unsigned i,const VectorGeneric<m> & r) {
306 12302388 : for(unsigned j=0; j<m; ++j) (*this)(i,j)=r(j);
307 3075597 : return *this;
308 : }
309 :
310 : template<unsigned n,unsigned m>
311 12312882 : VectorGeneric<n> TensorGeneric<n,m>::getCol(unsigned j)const {
312 12312882 : VectorGeneric<n> v;
313 61455240 : for(unsigned i=0; i<n; ++i) v(i)=(*this)(i,j);
314 12312882 : return v;
315 : }
316 :
317 : template<unsigned n,unsigned m>
318 3536124 : VectorGeneric<m> TensorGeneric<n,m>::getRow(unsigned i)const {
319 3536124 : VectorGeneric<m> v;
320 14144496 : for(unsigned j=0; j<m; ++j) v(j)=(*this)(i,j);
321 3536124 : return v;
322 : }
323 :
324 : template<unsigned n,unsigned m>
325 2932161 : TensorGeneric<n,m> operator+(const TensorGeneric<n,m>&t1,const TensorGeneric<n,m>&t2) {
326 2932161 : TensorGeneric<n,m> t(t1);
327 2932161 : t+=t2;
328 2932161 : return t;
329 : }
330 :
331 : template<unsigned n,unsigned m>
332 1821830 : TensorGeneric<n,m> operator-(const TensorGeneric<n,m>&t1,const TensorGeneric<n,m>&t2) {
333 1821830 : TensorGeneric<n,m> t(t1);
334 1821830 : t-=t2;
335 1821830 : return t;
336 : }
337 :
338 : template<unsigned n,unsigned m>
339 81592623 : TensorGeneric<n,m> operator*(const TensorGeneric<n,m>&t1,double s) {
340 81592623 : TensorGeneric<n,m> t(t1);
341 81592623 : t*=s;
342 81592623 : return t;
343 : }
344 :
345 : template<unsigned n,unsigned m>
346 55998693 : TensorGeneric<n,m> operator*(double s,const TensorGeneric<n,m>&t1) {
347 55998693 : return t1*s;
348 : }
349 :
350 : template<unsigned n,unsigned m>
351 37375 : TensorGeneric<n,m> operator/(const TensorGeneric<n,m>&t1,double s) {
352 37375 : return t1*(1.0/s);
353 : }
354 :
355 : template<>
356 : inline
357 272145 : double TensorGeneric<3,3>::determinant()const {
358 : return
359 272145 : d[0]*d[4]*d[8]
360 272145 : + d[1]*d[5]*d[6]
361 272145 : + d[2]*d[3]*d[7]
362 272145 : - d[0]*d[5]*d[7]
363 272145 : - d[1]*d[3]*d[8]
364 272145 : - d[2]*d[4]*d[6];
365 : }
366 :
367 : template<unsigned n,unsigned m>
368 : inline
369 13275572 : TensorGeneric<n,n> TensorGeneric<n,m>::identity() {
370 13275572 : TensorGeneric<n,n> t;
371 53102288 : for(unsigned i=0; i<n; i++) t(i,i)=1.0;
372 13275572 : return t;
373 : }
374 :
375 : template<unsigned n,unsigned m>
376 4882753 : TensorGeneric<m,n> TensorGeneric<n,m>::transpose()const {
377 4882753 : TensorGeneric<m,n> t;
378 63475789 : for(unsigned i=0; i<m; i++)for(unsigned j=0; j<n; j++) t(i,j)=(*this)(j,i);
379 4882753 : return t;
380 : }
381 :
382 : template<>
383 : inline
384 178307 : TensorGeneric<3,3> TensorGeneric<3,3>::inverse()const {
385 178307 : TensorGeneric t;
386 178307 : double invdet=1.0/determinant();
387 2317991 : for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++)
388 1604763 : t(j,i)=invdet*( (*this)((i+1)%3,(j+1)%3)*(*this)((i+2)%3,(j+2)%3)
389 1604763 : -(*this)((i+1)%3,(j+2)%3)*(*this)((i+2)%3,(j+1)%3));
390 178307 : return t;
391 : }
392 :
393 : template<unsigned n,unsigned m,unsigned l>
394 5616481 : TensorGeneric<n,l> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b) {
395 5616481 : TensorGeneric<n,l> t;
396 224659240 : for(unsigned i=0; i<n; i++) for(unsigned j=0; j<l; j++) for(unsigned k=0; k<m; k++) {
397 151644987 : t(i,j)+=a(i,k)*b(k,j);
398 : }
399 5616481 : return t;
400 : }
401 :
402 : template<unsigned n,unsigned m>
403 60636639 : VectorGeneric<n> matmul(const TensorGeneric<n,m>&a,const VectorGeneric<m>&b) {
404 60636639 : VectorGeneric<n> t;
405 788276307 : for(unsigned i=0; i<n; i++) for(unsigned j=0; j<m; j++) t(i)+=a(i,j)*b(j);
406 60636639 : return t;
407 : }
408 :
409 : template<unsigned n,unsigned m>
410 107695733 : VectorGeneric<n> matmul(const VectorGeneric<m>&a,const TensorGeneric<m,n>&b) {
411 107695733 : VectorGeneric<n> t;
412 1401199211 : for(unsigned i=0; i<n; i++) for(unsigned j=0; j<m; j++) t(i)+=a(j)*b(j,i);
413 107695733 : return t;
414 : }
415 :
416 : template<unsigned n_>
417 : double matmul(const VectorGeneric<n_>&a,const VectorGeneric<n_>&b) {
418 : return dotProduct(a,b);
419 : }
420 :
421 : template<unsigned n,unsigned m,unsigned l,unsigned i>
422 : TensorGeneric<n,i> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b,const TensorGeneric<l,i>&c) {
423 : return matmul(matmul(a,b),c);
424 : }
425 :
426 : template<unsigned n,unsigned m,unsigned l>
427 363240 : VectorGeneric<n> matmul(const TensorGeneric<n,m>&a,const TensorGeneric<m,l>&b,const VectorGeneric<l>&c) {
428 363240 : return matmul(matmul(a,b),c);
429 : }
430 :
431 : template<unsigned n,unsigned m,unsigned l>
432 : VectorGeneric<l> matmul(const VectorGeneric<n>&a,const TensorGeneric<n,m>&b,const TensorGeneric<m,l>&c) {
433 : return matmul(matmul(a,b),c);
434 : }
435 :
436 : template<unsigned n,unsigned m>
437 : double matmul(const VectorGeneric<n>&a,const TensorGeneric<n,m>&b,const VectorGeneric<m>&c) {
438 : return matmul(matmul(a,b),c);
439 : }
440 :
441 : inline
442 : double determinant(const TensorGeneric<3,3>&t) {
443 417 : return t.determinant();
444 : }
445 :
446 : inline
447 : TensorGeneric<3,3> inverse(const TensorGeneric<3,3>&t) {
448 89354 : return t.inverse();
449 : }
450 :
451 : template<unsigned n,unsigned m>
452 867561 : TensorGeneric<n,m> transpose(const TensorGeneric<m,n>&t) {
453 867561 : return t.transpose();
454 : }
455 :
456 : template<unsigned n,unsigned m>
457 2611634 : TensorGeneric<n,m> extProduct(const VectorGeneric<n>&v1,const VectorGeneric<m>&v2) {
458 2611634 : return TensorGeneric<n,m>(v1,v2);
459 : }
460 :
461 : inline
462 5182135 : TensorGeneric<3,3> dcrossDv1(const VectorGeneric<3>&v1,const VectorGeneric<3>&v2) {
463 : (void) v1; // this is to avoid warnings. still the syntax of this function is a bit dummy...
464 : return TensorGeneric<3,3>(
465 : 0.0, v2[2],-v2[1],
466 5182135 : -v2[2], 0.0, v2[0],
467 5182135 : v2[1],-v2[0], 0.0);
468 : }
469 :
470 : inline
471 5759791 : TensorGeneric<3,3> dcrossDv2(const VectorGeneric<3>&v1,const VectorGeneric<3>&v2) {
472 : (void) v2; // this is to avoid warnings. still the syntax of this function is a bit dummy...
473 : return TensorGeneric<3,3>(
474 : 0.0,-v1[2],v1[1],
475 5759791 : v1[2],0.0,-v1[0],
476 5759791 : -v1[1],v1[0],0.0);
477 : }
478 :
479 : template<unsigned n,unsigned m>
480 0 : std::ostream & operator<<(std::ostream &os, const TensorGeneric<n,m>& t) {
481 0 : for(unsigned i=0; i<n; i++)for(unsigned j=0; j<m; j++) {
482 0 : if(i!=(n-1) || j!=(m-1)) os<<t(i,j)<<" ";
483 0 : else os<<t(i,j);
484 : }
485 0 : return os;
486 : }
487 :
488 : /// \ingroup TOOLBOX
489 : typedef TensorGeneric<1,1> Tensor1d;
490 : /// \ingroup TOOLBOX
491 : typedef TensorGeneric<2,2> Tensor2d;
492 : /// \ingroup TOOLBOX
493 : typedef TensorGeneric<3,3> Tensor3d;
494 : /// \ingroup TOOLBOX
495 : typedef TensorGeneric<4,4> Tensor4d;
496 : /// \ingroup TOOLBOX
497 : typedef TensorGeneric<5,5> Tensor5d;
498 : /// \ingroup TOOLBOX
499 : typedef Tensor3d Tensor;
500 :
501 : inline
502 144414 : TensorGeneric<3,3> VcrossTensor(const VectorGeneric<3>&v1,const TensorGeneric<3,3>&v2) {
503 :
504 144414 : TensorGeneric<3,3> t;
505 577656 : for(unsigned i=0; i<3; i++) {
506 433242 : t.setRow(i,matmul(dcrossDv2(v1,v1),v2.getRow(i)));
507 : }
508 144414 : return t;
509 : }
510 :
511 : inline
512 48138 : TensorGeneric<3,3> VcrossTensor(const TensorGeneric<3,3>&v2,const VectorGeneric<3>&v1) {
513 48138 : TensorGeneric<3,3> t;
514 192552 : for(unsigned i=0; i<3; i++) {
515 144414 : t.setRow(i,-matmul(dcrossDv2(v1,v1),v2.getRow(i)));
516 : }
517 48138 : return t;
518 : }
519 :
520 :
521 : inline
522 48138 : TensorGeneric<3,3> deriNorm(const VectorGeneric<3>&v1,const TensorGeneric<3,3>&v2) {
523 : // delta(v) = delta(v1/v1.norm) = 1/v1.norm*(delta(v1) - (v.delta(v1))cross v;
524 48138 : double over_norm = 1./v1.modulo();
525 48138 : return over_norm*(v2 - over_norm*over_norm*(extProduct(matmul(v2,v1),v1)));
526 : }
527 :
528 : template<unsigned n,unsigned m>
529 623758 : void diagMatSym(const TensorGeneric<n,n>&mat,VectorGeneric<m>&evals,TensorGeneric<m,n>&evec) {
530 : // some guess number to make sure work is large enough.
531 : // for correctness it should be >=20. However, it is recommended to be the block size.
532 : // I put some likely exaggerated number
533 : constexpr int bs=100;
534 : // temporary data, on stack so as to avoid allocations
535 : std::array<int,10*n> iwork;
536 : std::array<double,(6+bs)*n> work;
537 : std::array<int,2*m> isup;
538 : // documentation says that "matrix is destroyed" !!!
539 623758 : auto mat_copy=mat;
540 : // documentation says this is size n (even if m<n)
541 : std::array<double,n> evals_tmp;
542 623758 : int nn=n; // dimension of matrix
543 623758 : double vl=0.0, vu=1.0; // ranges - not used
544 623758 : int one=1,mm=m; // minimum and maximum index
545 623758 : double abstol=0.0; // tolerance
546 623758 : int mout=0; // number of eigenvalues found (same as mm)
547 623758 : int info=0; // result
548 623758 : int liwork=iwork.size();
549 623758 : int lwork=work.size();
550 623758 : TensorGenericAux::local_dsyevr("V", (n==m?"A":"I"), "U", &nn, const_cast<double*>(&mat_copy[0][0]), &nn, &vl, &vu, &one, &mm,
551 623758 : &abstol, &mout, &evals_tmp[0], &evec[0][0], &nn,
552 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
553 623758 : if(info!=0) plumed_error()<<"Error diagonalizing matrix\n"
554 : <<"Matrix:\n"<<mat<<"\n"
555 : <<"Info: "<<info<<"\n";
556 623758 : plumed_assert(mout==m);
557 1424758 : for(unsigned i=0; i<m; i++) evals[i]=evals_tmp[i];
558 : // This changes eigenvectors so that the first non-null element
559 : // of each of them is positive
560 : // We can do it because the phase is arbitrary, and helps making
561 : // the result reproducible
562 1424758 : for(unsigned i=0; i<m; ++i) {
563 : unsigned j=0;
564 801593 : for(j=0; j<n; j++) if(evec(i,j)*evec(i,j)>1e-14) break;
565 1395931 : if(j<n) if(evec(i,j)<0.0) for(j=0; j<n; j++) evec(i,j)*=-1;
566 : }
567 623758 : }
568 :
569 : static_assert(sizeof(Tensor)==9*sizeof(double), "code cannot work if this is not satisfied");
570 :
571 : }
572 :
573 : #endif
574 :
|
