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 : #include <type_traits>
32 :
33 : namespace PLMD {
34 :
35 : //should I add something for calling ssyevr?
36 : /// Small class to contain local utilities.
37 : /// Should not be used outside of the TensorGeneric class.
38 : struct TensorGenericAux {
39 : // local redefinition, just to avoid including lapack.h here
40 : static void local_xsyevr(const char *jobz, const char *range, const char *uplo, int *n,
41 : double *a, int *lda, double *vl, double *vu, int *
42 : il, int *iu, double *abstol, int *m, double *w,
43 : double *z__, int *ldz, int *isuppz, double *work,
44 : int *lwork, int *iwork, int *liwork, int *info);
45 : static void local_xsyevr(const char *jobz, const char *range, const char *uplo, int *n,
46 : float *a, int *lda, float *vl, float *vu, int *
47 : il, int *iu, float *abstol, int *m, float *w,
48 : float *z__, int *ldz, int *isuppz, float *work,
49 : int *lwork, int *iwork, int *liwork, int *info);
50 : };
51 :
52 : /**
53 : \ingroup TOOLBOX
54 : Class implementing fixed size matrices of doubles
55 :
56 : \tparam n The number rows
57 : \tparam m The number columns
58 :
59 : This class implements a matrix of doubles with size fixed at
60 : compile time. It is useful for small fixed size objects (e.g.
61 : 3x3 tensors) as it does not waste space to store the vector size.
62 : Moreover, as the compiler knows the size, it can be completely
63 : optimized inline.
64 : Most of the loops are explicitly unrolled using PLMD::LoopUnroller class
65 : Matrix elements are initialized to zero by default. Notice that
66 : this means that constructor is a bit slow. This point might change
67 : in future if we find performance issues.
68 : It takes advantage of MatrixSquareBracketsAccess to provide both
69 : () and [] syntax for access.
70 : Several functions are declared as friends even if not necessary so as to
71 : properly appear in Doxygen documentation.
72 :
73 : Aliases are defined to simplify common declarations (Tensor, Tensor2d, Tensor3d, Tensor4d).
74 : Also notice that some operations are only available for 3x3 tensors.
75 :
76 : Example of usage
77 : \verbatim
78 :
79 : #include "Tensor.h"
80 :
81 : using namespace PLMD;
82 :
83 : int main(){
84 : Tensor a;
85 : TensorGeneric<3,2> b;
86 : TensorGeneric<3,2> c=matmul(a,b);
87 : return 0;
88 : }
89 :
90 : \endverbatim
91 : */
92 : template<typename T, unsigned n, unsigned m> class TensorTyped;
93 :
94 : template<typename T, unsigned n, unsigned m>
95 : std::ostream & operator<<(std::ostream &os, const TensorTyped<T,n,m>& t);
96 :
97 : /// Compute the Frobenius norm 2.
98 : template<typename T, unsigned n, unsigned m>
99 : T frobeniusNorm(const TensorTyped<T,m,n>&t);
100 :
101 : /// Diagonalize tensor.
102 : /// Syntax is the same as Matrix::diagMat.
103 : /// In addition, it is possible to call if with m smaller than n. In this case,
104 : /// only the first (smaller) m eigenvalues and eigenvectors are retrieved.
105 : /// If case lapack fails (info!=0) it throws an exception.
106 : /// Notice that tensor is assumed to be symmetric!!!
107 :
108 : template<typename T, unsigned n,unsigned m>
109 : void diagMatSym(const TensorTyped<T,n,n>&,VectorTyped<T,m>&evals,TensorTyped<T,m,n>&evec);
110 : /// Compute lowest eigenvalue and eigenvector, using a branchless iterative implementation.
111 : /// This function is amenable for running on openACC.
112 : /// The accuracy could be controlled increasing the number of iterations. The default value (24)
113 : /// seems sufficient for most applications.
114 : /// In principle, it could be extended to compute m eigenvalues by:
115 : /// - computing the lowest
116 : /// - compute the proper projector and reduce the matrix to <n-1,n-1>
117 : /// - proceed iteratively
118 : /// The interface should then be likely the same as diagMatSym
119 : template<typename T, unsigned n>
120 : T lowestEigenpairSym(const TensorTyped<T,n, n>& K, VectorTyped<T,n>& eigenvector, unsigned niter=24);
121 :
122 : template<typename T, unsigned n, unsigned m>
123 : class TensorTyped:
124 : public MatrixSquareBracketsAccess<TensorTyped<T,n,m>,T> {
125 : std::array<T,n*m> d;
126 : /// Auxiliary private function for constructor
127 : constexpr void auxiliaryConstructor();
128 : /// Auxiliary private function for constructor
129 : template<typename... Args>
130 : constexpr void auxiliaryConstructor(T first,Args... arg);
131 : template <int n_, int m_>
132 : using is3x3 = std::enable_if_t<(n_==3 && m_==3)>;
133 : public:
134 : //upload the tensor to the current acc device
135 : void toACCDevice()const {
136 : #pragma acc enter data copyin(this[0:1], d, d[0:n*m-1])
137 : }
138 : void removeFromACCDevice() const {
139 : // just delete
140 : #pragma acc exit data delete( d[0:n*m-1], d, this[0:1])
141 : }
142 : /// Constructor accepting n*m T parameters.
143 : /// Can be used as Tensor<2,2>(1.0,2.0,3.0,4.0)
144 : /// In case a wrong number of parameters is given, a static assertion will fail.
145 : template<typename... Args>
146 : constexpr TensorTyped(T first,Args... arg);
147 : /// initialize the tensor to zero
148 : constexpr TensorTyped();
149 : /// initialize a tensor as an external product of two Vector
150 : constexpr TensorTyped(const VectorTyped<T,n>&v1,const VectorTyped<T,m>&v2);
151 : /// set it to zero
152 : constexpr void zero();
153 : /// get the underline pointer to data
154 : constexpr T* data();
155 : /// get the underline pointer to data
156 : constexpr const T* data() const;
157 : /// access element
158 : constexpr T & operator() (unsigned i,unsigned j);
159 : /// access element
160 : constexpr const T & operator() (unsigned i,unsigned j)const;
161 : /// increment
162 : constexpr TensorTyped& operator +=(const TensorTyped& b);
163 : /// decrement
164 : constexpr TensorTyped& operator -=(const TensorTyped& b);
165 : /// multiply
166 : constexpr TensorTyped& operator *=(T);
167 : /// divide
168 : constexpr TensorTyped& operator /=(T);
169 : /// return +t
170 : constexpr TensorTyped operator +()const;
171 : /// return -t
172 : constexpr TensorTyped operator -()const;
173 : /// set j-th column
174 : constexpr TensorTyped& setCol(unsigned j,const VectorTyped<T,n> & c);
175 : /// set i-th row
176 : constexpr TensorTyped& setRow(unsigned i,const VectorTyped<T,m> & r);
177 : /// get j-th column
178 : constexpr VectorTyped<T,n> getCol(unsigned j)const;
179 : /// get i-th row
180 : constexpr VectorTyped<T,m> getRow(unsigned i)const;
181 : /// returns the determinant
182 : template<int m_=m, int n_=n,typename = is3x3<m_,n_>>
183 : constexpr T determinant()const;
184 : /// return an identity tensor
185 : static constexpr TensorTyped<T,n,n> identity();
186 : /// return the matrix inverse
187 : template<int m_=m, int n_=n,typename = is3x3<m_,n_>>
188 : constexpr TensorTyped inverse()const;
189 : /// return the transpose matrix
190 : constexpr TensorTyped<T,m,n> transpose()const;
191 : /// << operator.
192 : /// Allows printing tensor `t` with `std::cout<<t;`
193 : friend std::ostream & operator<< <>(std::ostream &os, const TensorTyped&);
194 : };
195 :
196 : template<typename T, unsigned n, unsigned m>
197 : constexpr void TensorTyped<T,n,m>::auxiliaryConstructor()
198 : {}
199 :
200 : template<typename T, unsigned n, unsigned m>
201 : template<typename... Args>
202 : constexpr void TensorTyped<T,n,m>::auxiliaryConstructor(T first,Args... arg) {
203 22369772 : d[n*m-(sizeof...(Args))-1]=first;
204 0 : auxiliaryConstructor(arg...);
205 : }
206 :
207 : template<typename T, unsigned n, unsigned m>
208 : template<typename... Args>
209 3156452 : constexpr TensorTyped<T,n,m>::TensorTyped(T first,Args... arg) {
210 : static_assert((sizeof...(Args))+1==n*m,"you are trying to initialize a Tensor with the wrong number of arguments");
211 : auxiliaryConstructor(first,arg...);
212 : }
213 :
214 : template<typename T, unsigned n, unsigned m>
215 : constexpr T* TensorTyped<T,n,m>::data() {
216 : return d.data();
217 : }
218 : template<typename T, unsigned n, unsigned m>
219 : constexpr const T* TensorTyped<T,n,m>::data() const {
220 : return d.data();
221 : }
222 :
223 : template<typename T, unsigned n, unsigned m>
224 11324708 : constexpr TensorTyped<T,n,m>::TensorTyped() {
225 : LoopUnroller<n*m>::_zero(d.data());
226 : }
227 :
228 : /* between RVO and compile time this should be faster, but slows down openACC, a lot
229 : template <typename T, unsigned i, unsigned j, unsigned m>
230 : void external_rec(T*const out,const T*const v1, const T*const v2){
231 : if constexpr (j>0) {
232 : external_rec<i,j-1,m>(out,v1,v2);
233 : } else if constexpr (i>0) {
234 : external_rec<i-1,m-1,m>(out,v1,v2);
235 : }
236 : out[i*m+j]=v1[i]*v2[j];
237 : }
238 : template<typename T, unsigned n, unsigned m>
239 : std::array<T,n*m> externaProd(const VectorGeneric<T,n>&v1,const VectorGeneric<T,m>&v2){
240 : std::array<T,n*m> toRet;
241 : external_rec<n-1,m-1,m>(toRet.data(),v1.data(),v2.data());
242 : return toRet;
243 : }
244 : template<typename T, unsigned n, unsigned m>
245 : TensorGeneric<T,n,m>::TensorGeneric(const VectorGeneric<T,n>&v1,const VectorGeneric<T,m>&v2)
246 : :d(externaProd(v1,v2)) {}
247 : */
248 :
249 : template<typename T, unsigned n, unsigned m>
250 111639276 : constexpr TensorTyped<T,n,m>::TensorTyped(const VectorTyped<T,n>&v1,const VectorTyped<T,m>&v2) {
251 446557104 : for(unsigned i=0; i<n; i++) {
252 1339671312 : for(unsigned j=0; j<m; j++) {
253 1004753484 : d[i*m+j]=v1[i]*v2[j];
254 : }
255 : }
256 111639276 : }
257 :
258 : template<typename T, unsigned n, unsigned m>
259 : constexpr void TensorTyped<T,n,m>::zero() {
260 : LoopUnroller<n*m>::_zero(d.data());
261 100 : }
262 :
263 : template<typename T, unsigned n, unsigned m>
264 : constexpr T & TensorTyped<T,n,m>::operator() (unsigned i,unsigned j) {
265 : #ifdef _GLIBCXX_DEBUG
266 : // index i is implicitly checked by the std::array class
267 : plumed_assert(j<m);
268 : #endif
269 47479734 : return d[m*i+j];
270 : }
271 :
272 : template<typename T, unsigned n, unsigned m>
273 : constexpr const T & TensorTyped<T,n,m>::operator() (unsigned i,unsigned j)const {
274 : #ifdef _GLIBCXX_DEBUG
275 : // index i is implicitly checked by the std::array class
276 : plumed_assert(j<m);
277 : #endif
278 3772500020 : return d[m*i+j];
279 : }
280 :
281 : template<typename T, unsigned n, unsigned m>
282 : constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::operator +=(const TensorTyped<T,n,m>& b) {
283 : LoopUnroller<n*m>::_add(d.data(),b.d.data());
284 : return *this;
285 : }
286 :
287 : template<typename T, unsigned n, unsigned m>
288 : constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::operator -=(const TensorTyped<T,n,m>& b) {
289 : LoopUnroller<n*m>::_sub(d.data(),b.d.data());
290 : return *this;
291 : }
292 :
293 : template<typename T, unsigned n, unsigned m>
294 : constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::operator *=(T s) {
295 77829008 : LoopUnroller<n*m>::_mul(d.data(),s);
296 : return *this;
297 : }
298 :
299 : template<typename T, unsigned n, unsigned m>
300 : constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::operator /=(T s) {
301 18636720 : LoopUnroller<n*m>::_mul(d.data(),1.0/s);
302 : return *this;
303 : }
304 :
305 : template<typename T, unsigned n, unsigned m>
306 : constexpr TensorTyped<T,n,m> TensorTyped<T,n,m>::operator+()const {
307 449931 : return *this;
308 : }
309 :
310 : template<typename T, unsigned n, unsigned m>
311 : constexpr TensorTyped<T,n,m> TensorTyped<T,n,m>::operator-()const {
312 1945850 : TensorTyped<T,n,m> r;
313 1619251 : LoopUnroller<n*m>::_neg(r.d.data(),d.data());
314 : return r;
315 : }
316 :
317 : template<typename T, unsigned n, unsigned m>
318 : constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::setCol(unsigned j,const VectorTyped<T,n> & c) {
319 485892 : for(unsigned i=0; i<n; ++i) {
320 364419 : (*this)(i,j)=c(i);
321 : }
322 : return *this;
323 : }
324 :
325 : template<typename T, unsigned n, unsigned m>
326 : constexpr TensorTyped<T,n,m>& TensorTyped<T,n,m>::setRow(unsigned i,const VectorTyped<T,m> & r) {
327 12541938 : for(unsigned j=0; j<m; ++j) {
328 9440397 : (*this)(i,j)=r(j);
329 : }
330 : return *this;
331 : }
332 :
333 : template<typename T, unsigned n, unsigned m>
334 : constexpr VectorTyped<T,n> TensorTyped<T,n,m>::getCol(unsigned j)const {
335 10458675 : VectorTyped<T,n> v;
336 51303258 : for(unsigned i=0; i<n; ++i) {
337 43200045 : v(i)=(*this)(i,j);
338 : }
339 : return v;
340 : }
341 :
342 : template<typename T, unsigned n, unsigned m>
343 : constexpr VectorTyped<T,m> TensorTyped<T,n,m>::getRow(unsigned i)const {
344 9840503 : VectorTyped<T,m> v;
345 44363848 : for(unsigned j=0; j<m; ++j) {
346 35733749 : v(j)=(*this)(i,j);
347 : }
348 : return v;
349 : }
350 :
351 : template<typename T, unsigned n, unsigned m>
352 : constexpr TensorTyped<T,n,m> operator+(const TensorTyped<T,n,m>&t1,const TensorTyped<T,n,m>&t2) {
353 4162891 : TensorTyped<T,n,m> t(t1);
354 : t+=t2;
355 : return t;
356 : }
357 :
358 : template<typename T, unsigned n, unsigned m>
359 : constexpr TensorTyped<T,n,m> operator-(const TensorTyped<T,n,m>&t1,const TensorTyped<T,n,m>&t2) {
360 21074803 : TensorTyped<T,n,m> t(t1);
361 : t-=t2;
362 : return t;
363 : }
364 :
365 : template<typename T, typename J, unsigned n, unsigned m>
366 : constexpr TensorTyped<T,n,m> operator*(const TensorTyped<T,n,m>&t1,J s) {
367 73951680 : TensorTyped<T,n,m> t(t1);
368 : t*=s;
369 : return t;
370 : }
371 :
372 : template<typename T, typename J, unsigned n, unsigned m>
373 : constexpr TensorTyped<T,n,m> operator*(J s,const TensorTyped<T,n,m>&t1) {
374 : return t1*s;
375 : }
376 :
377 : template<typename T, typename J, unsigned n, unsigned m>
378 : constexpr TensorTyped<T,n,m> operator/(const TensorTyped<T,n,m>&t1,J s) {
379 21771 : return t1*(T(1.0)/s);
380 : }
381 :
382 : template<typename T, unsigned n, unsigned m>
383 : template<int m_, int n_,typename>
384 568457 : constexpr inline T TensorTyped<T,n,m>::determinant()const {
385 : return
386 568457 : d[0]*d[4]*d[8]
387 568457 : + d[1]*d[5]*d[6]
388 568457 : + d[2]*d[3]*d[7]
389 568457 : - d[0]*d[5]*d[7]
390 568457 : - d[1]*d[3]*d[8]
391 568457 : - d[2]*d[4]*d[6];
392 : }
393 :
394 : //consider to make this a constexpr function
395 : template<typename T, unsigned n, unsigned m>
396 : constexpr inline TensorTyped<T,n,n> TensorTyped<T,n,m>::identity() {
397 : TensorTyped<T,n,n> t;
398 : for(unsigned i=0; i<n; i++) {
399 : t(i,i)=1.0;
400 : }
401 : return t;
402 : }
403 :
404 : template<typename T, unsigned n, unsigned m>
405 6975960 : constexpr TensorTyped<T,m,n> TensorTyped<T,n,m>::transpose()const {
406 6975960 : TensorTyped<T,m,n> t;
407 27903840 : for(unsigned i=0; i<m; i++)
408 83711520 : for(unsigned j=0; j<n; j++) {
409 62783640 : t(i,j)=(*this)(j,i);
410 : }
411 6975960 : return t;
412 : }
413 :
414 : template<typename T, unsigned n, unsigned m>
415 : template<int m_, int n_,typename>
416 364767 : constexpr inline TensorTyped<T,n,m> TensorTyped<T,n,m>::inverse()const {
417 364767 : TensorTyped t;
418 364767 : T invdet=1.0/determinant();
419 1459068 : for(unsigned i=0; i<3; i++)
420 4377204 : for(unsigned j=0; j<3; j++)
421 3282903 : t(j,i)=invdet*( (*this)((i+1)%3,(j+1)%3)*(*this)((i+2)%3,(j+2)%3)
422 3282903 : -(*this)((i+1)%3,(j+2)%3)*(*this)((i+2)%3,(j+1)%3));
423 364767 : return t;
424 : }
425 :
426 : /// matrix-matrix multiplication
427 : template<typename T, unsigned n, unsigned m, unsigned l>
428 23572032 : constexpr TensorTyped<T,n,l> matmul(const TensorTyped<T,n,m>&a,const TensorTyped<T,m,l>&b) {
429 23572032 : TensorTyped<T,n,l> t;
430 112924848 : for(unsigned i=0; i<n; i++)
431 431958144 : for(unsigned j=0; j<l; j++)
432 1668608832 : for(unsigned k=0; k<m; k++) {
433 1326003504 : t(i,j)+=a(i,k)*b(k,j);
434 : }
435 23572032 : return t;
436 : }
437 :
438 : /// matrix-vector multiplicationi, the first argument determines the type of the result
439 : template<typename T,typename TT, unsigned n, unsigned m>
440 80117374 : constexpr VectorTyped<T,n> matmul(const TensorTyped<T,n,m>&a,const VectorTyped<TT,m>&b) {
441 80117374 : VectorTyped<T,n> t;
442 320469496 : for(unsigned i=0; i<n; i++)
443 961408488 : for(unsigned j=0; j<m; j++) {
444 721056366 : t(i)+=a(i,j)*b(j);
445 : }
446 80117374 : return t;
447 : }
448 :
449 : /// vector-matrix multiplication, the first argument determines the type of the result
450 : template<typename T,typename TT, unsigned n, unsigned m>
451 126614546 : constexpr VectorTyped<T,n> matmul(const VectorTyped<T,m>&a,const TensorTyped<TT,m,n>&b) {
452 126614546 : VectorTyped<T,n> t;
453 507234714 : for(unsigned i=0; i<n; i++)
454 1526741690 : for(unsigned j=0; j<m; j++) {
455 1146121522 : t(i)+=a(j)*b(j,i);
456 : }
457 126614546 : return t;
458 : }
459 :
460 : /// vector-vector multiplication (maps to dotProduct)
461 : template<typename T, unsigned n_>
462 : constexpr T matmul(const VectorTyped<T,n_>&a,const VectorTyped<T,n_>&b) {
463 : return dotProduct(a,b);
464 : }
465 :
466 : /// matrix-matrix-matrix multiplication
467 : template<typename T, unsigned n, unsigned m, unsigned l, unsigned i>
468 : constexpr TensorTyped<T,n,i> matmul(const TensorTyped<T,n,m>&a,const TensorTyped<T,m,l>&b,const TensorTyped<T,l,i>&c) {
469 : return matmul(matmul(a,b),c);
470 : }
471 :
472 : /// matrix-matrix-vector multiplication
473 : template<typename T, unsigned n, unsigned m, unsigned l>
474 407268 : constexpr VectorTyped<T,n> matmul(const TensorTyped<T,n,m>&a,const TensorTyped<T,m,l>&b,const VectorTyped<T,l>&c) {
475 407268 : return matmul(matmul(a,b),c);
476 : }
477 :
478 : /// vector-matrix-matrix multiplication
479 : template<typename T, unsigned n, unsigned m, unsigned l>
480 : constexpr VectorTyped<T,l> matmul(const VectorTyped<T,n>&a,const TensorTyped<T,n,m>&b,const TensorTyped<T,m,l>&c) {
481 : return matmul(matmul(a,b),c);
482 : }
483 :
484 : /// vector-matrix-vector multiplication
485 : template<typename T, unsigned n, unsigned m>
486 776530 : constexpr T matmul(const VectorTyped<T,n>&a,const TensorTyped<T,n,m>&b,const VectorTyped<T,m>&c) {
487 776530 : return matmul(matmul(a,b),c);
488 : }
489 :
490 : template <typename T>
491 : inline
492 : constexpr T determinant(const TensorTyped<T,3,3>&t) {
493 817 : return t.determinant();
494 : }
495 :
496 : template <typename T>
497 : inline
498 : constexpr TensorTyped<T,3,3> inverse(const TensorTyped<T,3,3>&t) {
499 182596 : return t.inverse();
500 : }
501 :
502 : /// returns the transpose of a tensor (same as TensorGeneric::transpose())
503 : template<typename T, unsigned n, unsigned m>
504 : constexpr TensorTyped<T,n,m> transpose(const TensorTyped<T,m,n>&t) {
505 944580 : return t.transpose();
506 : }
507 :
508 : /// returns the transpose of a tensor (same as TensorGeneric(const VectorGeneric&,const VectorGeneric&))
509 : template<typename T, unsigned n, unsigned m>
510 : constexpr TensorTyped<T,n,m> extProduct(const VectorTyped<T,n>&v1,const VectorTyped<T,m>&v2) {
511 3566376 : return TensorTyped<T,n,m>(v1,v2);
512 : }
513 :
514 : template <typename T>
515 : inline
516 : constexpr TensorTyped<T,3,3> dcrossDv1(const VectorTyped<T,3>&v1,const VectorTyped<T,3>&v2) {
517 : (void) v1; // this is to avoid warnings. still the syntax of this function is a bit dummy...
518 : return TensorTyped<T,3,3>(
519 : T(0.0), v2[2], -v2[1],
520 : -v2[2], T(0.0), v2[0],
521 1112903 : v2[1], -v2[0], T(0.0));
522 : }
523 :
524 : template <typename T>
525 : inline
526 : constexpr TensorTyped<T,3,3> dcrossDv2(const VectorTyped<T,3>&v1,const VectorTyped<T,3>&v2) {
527 : (void) v2; // this is to avoid warnings. still the syntax of this function is a bit dummy...
528 : return TensorTyped<T,3,3>(
529 : T(0.0),-v1[2], v1[1],
530 : v1[2], T(0.0),-v1[0],
531 1691135 : -v1[1],v1[0], T(0.0));
532 : }
533 :
534 : template<typename T, unsigned n,unsigned m>
535 18636720 : T frobeniusNorm(const TensorTyped<T,m,n>&t) {
536 18636720 : return std::sqrt(LoopUnroller<m*n>::_dot(t.data(),t.data()));
537 : }
538 :
539 : template<typename T, unsigned n, unsigned m>
540 0 : std::ostream & operator<<(std::ostream &os, const TensorTyped<T,n,m>& t) {
541 0 : for(unsigned i=0; i<n; i++)
542 0 : for(unsigned j=0; j<m; j++) {
543 0 : if(i!=(n-1) || j!=(m-1)) {
544 0 : os<<t(i,j)<<" ";
545 : } else {
546 0 : os<<t(i,j);
547 : }
548 : }
549 0 : return os;
550 : }
551 :
552 : template <unsigned n,unsigned m>
553 : using TensorGeneric=TensorTyped<double,n,m>;
554 :
555 : /// \ingroup TOOLBOX
556 : typedef TensorGeneric<1,1> Tensor1d;
557 : /// \ingroup TOOLBOX
558 : typedef TensorGeneric<2,2> Tensor2d;
559 : /// \ingroup TOOLBOX
560 : typedef TensorGeneric<3,3> Tensor3d;
561 : /// \ingroup TOOLBOX
562 : typedef TensorGeneric<4,4> Tensor4d;
563 : /// \ingroup TOOLBOX
564 : typedef TensorGeneric<5,5> Tensor5d;
565 : /// \ingroup TOOLBOX
566 : typedef Tensor3d Tensor;
567 :
568 : /// \ingroup TOOLBOX
569 : typedef TensorTyped<float,1,1> Tensor1f;
570 : /// \ingroup TOOLBOX
571 : typedef TensorTyped<float,2,2> Tensor2f;
572 : /// \ingroup TOOLBOX
573 : typedef TensorTyped<float,3,3> Tensor3f;
574 : /// \ingroup TOOLBOX
575 : typedef TensorTyped<float,4,4> Tensor4f;
576 : /// \ingroup TOOLBOX
577 : typedef TensorTyped<float,5,5> Tensor5f;
578 :
579 : /// \ingroup TOOLBOX
580 : /// Alias for three by three Tensor of any type:
581 : template <typename T>
582 : using TensorT=TensorTyped<T,3,3>;
583 :
584 : template <typename T>
585 : inline
586 144558 : TensorTyped<T,3,3> VcrossTensor(const VectorTyped<T,3>&v1,const TensorTyped<T,3,3>&v2) {
587 144558 : TensorTyped<T,3,3> t;
588 578232 : for(unsigned i=0; i<3; i++) {
589 867348 : t.setRow(i,matmul(dcrossDv2(v1,v1),v2.getRow(i)));
590 : }
591 144558 : return t;
592 : }
593 :
594 : template <typename T>
595 : inline
596 48186 : TensorTyped<T,3,3> VcrossTensor(const TensorTyped<T,3,3>&v2,const VectorTyped<T,3>&v1) {
597 48186 : TensorTyped<T,3,3> t;
598 192744 : for(unsigned i=0; i<3; i++) {
599 289116 : t.setRow(i,-matmul(dcrossDv2(v1,v1),v2.getRow(i)));
600 : }
601 48186 : return t;
602 : }
603 :
604 : template <typename T>
605 : inline
606 48186 : TensorTyped<T,3,3> deriNorm(const VectorTyped<T,3>&v1,const TensorTyped<T,3,3>&v2) {
607 : // delta(v) = delta(v1/v1.norm) = 1/v1.norm*(delta(v1) - (v.delta(v1))cross v;
608 48186 : T over_norm = 1./v1.modulo();
609 48186 : return over_norm*(v2 - over_norm*over_norm*(extProduct(matmul(v2,v1),v1)));
610 : }
611 :
612 : /// Diagonalize tensor.
613 : /// Syntax is the same as Matrix::diagMat.
614 : /// In addition, it is possible to call if with m_ smaller than n_. In this case,
615 : /// only the first (smaller) m_ eigenvalues and eigenvectors are retrieved.
616 : /// If case lapack fails (info!=0) it throws an exception.
617 : /// Notice that tensor is assumed to be symmetric!!!
618 : template<typename precision, unsigned n, unsigned m>
619 150843 : void diagMatSym(const TensorTyped<precision,n,n>&mat,VectorTyped<precision,m>&evals,TensorTyped<precision,m,n>&evec) {
620 : // some guess number to make sure work is large enough.
621 : // for correctness it should be >=20. However, it is recommended to be the block size.
622 : // I put some likely exaggerated number
623 : constexpr int bs=100;
624 : // temporary data, on stack so as to avoid allocations
625 : std::array<int,10*n> iwork;
626 : std::array<precision,(6+bs)*n> work;
627 : std::array<int,2*m> isup;
628 : // documentation says that "matrix is destroyed" !!!
629 150843 : auto mat_copy=mat;
630 : // documentation says this is size n (even if m<n)
631 : std::array<precision,n> evals_tmp;
632 150843 : int nn=n; // dimension of matrix
633 150843 : precision vl=0.0, vu=1.0; // ranges - not used
634 150843 : int one=1,mm=m; // minimum and maximum index
635 150843 : precision abstol=0.0; // tolerance
636 150843 : int mout=0; // number of eigenvalues found (same as mm)
637 150843 : int info=0; // result
638 150843 : int liwork=iwork.size();
639 150843 : int lwork=work.size();
640 150843 : TensorGenericAux::local_xsyevr("V", (n==m?"A":"I"), "U", &nn, mat_copy.data(), &nn, &vl, &vu, &one, &mm,
641 : &abstol, &mout, evals_tmp.data(), evec.data(), &nn,
642 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
643 150843 : if(info!=0)
644 0 : plumed_error()<<"Error diagonalizing matrix\n"
645 : <<"Matrix:\n"<<mat<<"\n"
646 : <<"Info: "<<info<<"\n";
647 150843 : plumed_assert(mout==m);
648 753349 : for(unsigned i=0; i<m; i++) {
649 602506 : evals[i]=evals_tmp[i];
650 : }
651 : // This changes eigenvectors so that the first non-null element
652 : // of each of them is positive
653 : // We can do it because the phase is arbitrary, and helps making
654 : // the result reproducible
655 753349 : for(unsigned i=0; i<m; ++i) {
656 : unsigned j=0;
657 605841 : for(j=0; j<n; j++)
658 605841 : if(evec(i,j)*evec(i,j)>1e-14) {
659 : break;
660 : }
661 602506 : if(j<n)
662 602506 : if(evec(i,j)<0.0)
663 738762 : for(j=0; j<n; j++) {
664 590816 : evec(i,j)*=-1;
665 : }
666 : }
667 150843 : }
668 :
669 : template<typename T, unsigned n>
670 776530 : T lowestEigenpairSym(const TensorTyped<T,n,n>& K, VectorTyped<T,n>& eigenvector, const unsigned niter) {
671 : // Estimate upper bound for largest eigenvalue using Gershgorin disks
672 : T upper_bound = std::numeric_limits<T>::lowest();
673 3882650 : for (unsigned i = 0; i < n; ++i) {
674 : auto row = K.getRow(i);
675 3106120 : T center = row[i];
676 3106120 : row[i] = 0.0; // zero out the diagonal entry
677 :
678 : // Compute sum of absolute values of the off-diagonal elements
679 : T row_sum = 0.0;
680 15530600 : for (unsigned j = 0; j < n; ++j) {
681 12424480 : row_sum += std::fabs(row[j]);
682 : }
683 :
684 3106120 : upper_bound = std::fmax(upper_bound, center + row_sum);
685 : }
686 :
687 : // Build shifted matrix A = -K + lambda_shift * I
688 : TensorTyped<T,n,n> A = -K;
689 3882650 : for (unsigned i = 0; i < n; ++i) {
690 3106120 : A[i][i] += upper_bound;
691 : }
692 :
693 : // Repeated squaring + normalization to project onto dominant eigenspace
694 19413250 : for (unsigned k = 0; k < niter; ++k) {
695 18636720 : A = matmul(A, A);
696 18636720 : auto frob = frobeniusNorm(A);
697 18636720 : A /= (frob > 0.0 ? frob : 1.0);
698 : }
699 :
700 : // Extract dominant eigenvector from projector A
701 776530 : VectorTyped<T,n> v;
702 : T best_norm2 = 0.0;
703 3882650 : for (unsigned j = 0; j < n; ++j) {
704 : auto row = A.getRow(j);
705 : T norm2 = modulo2(row);
706 3106120 : T use = static_cast<T>(norm2 > best_norm2);
707 3106120 : best_norm2 = use * norm2 + (1.0 - use) * best_norm2;
708 3106120 : v = use * row + (1.0 - use) * v;
709 : }
710 :
711 : v /= modulo(v);
712 :
713 776530 : eigenvector = v;
714 :
715 776530 : return matmul(eigenvector, K, eigenvector);
716 : }
717 :
718 :
719 : static_assert(sizeof(Tensor)==9*sizeof(double), "code cannot work if this is not satisfied");
720 : static_assert(sizeof(TensorT<float>)==9*sizeof(float), "code cannot work if this is not satisfied");
721 :
722 : } //PLMD
723 :
724 : #endif
|
