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_Matrix_h
23 : #define __PLUMED_tools_Matrix_h
24 : #include <vector>
25 : #include <string>
26 : #include <set>
27 : #include <cmath>
28 : #include "Exception.h"
29 : #include "MatrixSquareBracketsAccess.h"
30 : #include "Tools.h"
31 : #include "Log.h"
32 : #include "lapack/lapack.h"
33 :
34 : namespace PLMD {
35 :
36 : /// Calculate the dot product between two vectors
37 : template <typename T> T dotProduct( const std::vector<T>& A, const std::vector<T>& B ) {
38 : plumed_assert( A.size()==B.size() );
39 : T val{};
40 : for(unsigned i=0; i<A.size(); ++i) {
41 : val+=A[i]*B[i];
42 : }
43 : return val;
44 : }
45 :
46 : /// Calculate the dot product between a vector and itself
47 86 : template <typename T> T norm( const std::vector<T>& A ) {
48 : T val{};
49 608 : for(unsigned i=0; i<A.size(); ++i) {
50 522 : val+=A[i]*A[i];
51 : }
52 86 : return val;
53 : }
54 :
55 : /// This class stores a full matrix and allows one to do some simple matrix operations
56 : template <typename T>
57 1455812 : class Matrix:
58 : public MatrixSquareBracketsAccess<Matrix<T>,T> {
59 : /// Multiply matrix by scalar
60 : template <typename U> friend Matrix<U> operator*(U&, const Matrix<U>& );
61 : /// Matrix matrix multiply
62 : template <typename U> friend void mult( const Matrix<U>&, const Matrix<U>&, Matrix<U>& );
63 : /// Matrix times a std::vector
64 : template <typename U> friend void mult( const Matrix<U>&, const std::vector<U>&, std::vector<U>& );
65 : /// std::vector times a Matrix
66 : template <typename U> friend void mult( const std::vector<U>&, const Matrix<U>&, std::vector<U>& );
67 : /// Matrix transpose
68 : template <typename U> friend void transpose( const Matrix<U>&, Matrix<U>& );
69 : /// Output the entire matrix on a single line
70 : template <typename U> friend Log& operator<<(Log&, const Matrix<U>& );
71 : /// Output the Matrix in matrix form
72 : template <typename U> friend void matrixOut( Log&, const Matrix<U>& );
73 : /// Diagonalize a symmetric matrix - returns zero if diagonalization worked
74 : template <typename U> friend int diagMat( const Matrix<U>&, std::vector<U>&, Matrix<U>& );
75 : /// Diagonalize a symmetric matrix computing only the eigenpairs from first to last (0-based, ascending eigenvalues) - returns zero if diagonalization worked
76 : template <typename U> friend int diagMatRange( const Matrix<U>&, unsigned, unsigned, std::vector<U>&, Matrix<U>& );
77 : /// Calculate the Moore-Penrose Pseudoinverse of a matrix
78 : template <typename U> friend int pseudoInvert( const Matrix<U>&, Matrix<U>& );
79 : /// Calculate the logarithm of the determinant of a symmetric matrix - returns zero if successful
80 : template <typename U> friend int logdet( const Matrix<U>&, U& );
81 : /// Invert a matrix (works for both symmetric and asymmetric matrices) - returns zero if sucesfull
82 : template <typename U> friend int Invert( const Matrix<U>&, Matrix<U>& );
83 : /// Do a cholesky decomposition of a matrix
84 : template <typename U> friend void cholesky( const Matrix<U>&, Matrix<U>& );
85 : /// Solve a system of equations using the cholesky decomposition
86 : template <typename U> friend void chol_elsolve( const Matrix<U>&, const std::vector<U>&, std::vector<U>& );
87 : private:
88 : /// Number of elements in matrix (nrows*ncols)
89 : unsigned sz;
90 : /// Number of rows in matrix
91 : unsigned rw;
92 : /// Number of columns in matrix
93 : unsigned cl;
94 : /// The data in the matrix
95 : std::vector<T> data;
96 : template <typename Y>
97 : struct openAccVectorHelper {
98 : using pointer = Y*;
99 : size_t sz=0;
100 : pointer data=nullptr;
101 : void toACCDevice() const {
102 : #pragma acc enter data copyin(this[0:1],sz,data[0:sz])
103 : }
104 : void removeFromACCDevice() const {
105 : #pragma acc exit data delete(data[0:sz],sz,this[0:1])
106 : }
107 : void update(std::vector<Y>&x) {
108 1786104 : data = x.data();
109 1786104 : sz = x.size();
110 : }
111 : };
112 : openAccVectorHelper<T> helper{};
113 :
114 : public:
115 1635225 : explicit Matrix(const unsigned nr=0, const unsigned nc=0 ) : sz(nr*nc), rw(nr), cl(nc), data(nr*nc) {
116 : helper.update(data);
117 1635225 : }
118 : void toACCDevice() const {
119 : #pragma acc enter data copyin(this[0:1])
120 : helper.toACCDevice();
121 : }
122 : void removeFromACCDevice() const {
123 : helper.removeFromACCDevice();
124 : // and delete
125 : #pragma acc exit data delete( this[0:1])
126 : }
127 612 : Matrix(const Matrix<T>& t) : sz(t.sz), rw(t.rw), cl(t.cl), data(t.data) {
128 : helper.update(data);
129 612 : }
130 : /// Resize the matrix
131 1542 : void resize( const unsigned nr, const unsigned nc ) {
132 1542 : rw=nr;
133 1542 : cl=nc;
134 1542 : sz=nr*nc;
135 1542 : data.resize(sz);
136 : helper.update(data);
137 1542 : }
138 : /// Return the number of rows
139 : inline unsigned nrows() const {
140 16901187 : return rw;
141 : }
142 : /// Return the number of columns
143 : inline unsigned ncols() const {
144 11483960 : return cl;
145 : }
146 : /// Return the contents of the matrix as a vector of length rw*cl
147 : inline std::vector<T>& getVector() {
148 : return data;
149 : }
150 : /// Set the matrix from a vector input
151 : inline void setFromVector( const std::vector<T>& vecin ) {
152 : plumed_assert( vecin.size()==sz );
153 : for(unsigned i=0; i<sz; ++i) {
154 : data[i]=vecin[i];
155 : }
156 : }
157 : /// Return element i,j of the matrix
158 : inline const T& operator () (const unsigned& i, const unsigned& j) const {
159 18246495 : return helper.data[j+i*cl];
160 : }
161 : /// Return a referenre to element i,j of the matrix
162 : inline T& operator () (const unsigned& i, const unsigned& j) {
163 68390578 : return helper.data[j+i*cl];
164 : }
165 : /// Set all elements of the matrix equal to the value of v
166 : Matrix<T>& operator=(const T& v) {
167 3010243975 : for(unsigned i=0; i<sz; ++i) {
168 3010178103 : data[i]=v;
169 : }
170 : return *this;
171 : }
172 : /// Set the Matrix equal to another Matrix
173 148725 : Matrix<T>& operator=(const Matrix<T>& m) {
174 148725 : sz=m.sz;
175 148725 : rw=m.rw;
176 148725 : cl=m.cl;
177 148725 : data=m.data;
178 : helper.update(data);
179 148725 : return *this;
180 : }
181 : /// Set the Matrix equal to the value of a standard vector - used for readin
182 : Matrix<T>& operator=(const std::vector<T>& v) {
183 : plumed_dbg_assert( v.size()==sz );
184 : for(unsigned i=0; i<sz; ++i) {
185 : data[i]=v[i];
186 : }
187 : return *this;
188 : }
189 : /// Add v to all elements of the Matrix
190 : Matrix<T> operator+=(const T& v) {
191 : for(unsigned i=0; i<sz; ++i) {
192 : data[i]+=v;
193 : }
194 : return *this;
195 : }
196 : /// Multiply all elements by v
197 125 : Matrix<T> operator*=(const T& v) {
198 55038 : for(unsigned i=0; i<sz; ++i) {
199 54913 : data[i]*=v;
200 : }
201 125 : return *this;
202 : }
203 : /// Matrix addition
204 : Matrix<T>& operator+=(const Matrix<T>& m) {
205 : plumed_dbg_assert( m.rw==rw && m.cl==cl );
206 : data+=m.data;
207 : return *this;
208 : }
209 : /// Subtract v from all elements of the Matrix
210 : Matrix<T> operator-=(const T& v) {
211 : for(unsigned i=0; i<sz; ++i) {
212 : data-=v;
213 : }
214 : return *this;
215 : }
216 : /// Matrix subtraction
217 : Matrix<T>& operator-=(const Matrix<T>& m) {
218 : plumed_dbg_assert( m.rw==rw && m.cl==cl );
219 : data-=m.data;
220 : return *this;
221 : }
222 : /// Test if the matrix is symmetric or not
223 160688 : unsigned isSymmetric() const {
224 160688 : if (rw!=cl) {
225 : return 0;
226 : }
227 : unsigned sym=1;
228 319266 : for(unsigned i=1; i<rw; ++i)
229 374415 : for(unsigned j=0; j<i; ++j)
230 350119 : if( std::fabs(data[i+j*cl]-data[j+i*cl])>1.e-10 ) {
231 : sym=0;
232 : break;
233 : }
234 : return sym;
235 : }
236 : };
237 :
238 : /// Multiply matrix by scalar
239 : template <typename T> Matrix<T> operator*(T& v, const Matrix<T>& m ) {
240 : Matrix<T> new_m(m);
241 : new_m*=v;
242 : return new_m;
243 : }
244 :
245 8572 : template <typename T> void mult( const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& C ) {
246 8572 : plumed_assert(A.cl==B.rw);
247 8572 : if( A.rw !=C.rw || B.cl !=C.cl ) {
248 36 : C.resize( A.rw, B.cl );
249 : }
250 : C=static_cast<T>( 0 );
251 : // i-k-j loop order: the innermost loop runs over contiguous elements of B and C,
252 : // so it is cache friendly and can be vectorized. For each element of C the
253 : // products are summed in the same order as with the i-j-k loop order.
254 : const T* a=A.data.data();
255 : const T* b=B.data.data();
256 : T* c=C.data.data();
257 71524 : for(unsigned i=0; i<A.rw; ++i) {
258 62952 : T* ci=c+i*B.cl;
259 7594384 : for(unsigned k=0; k<A.cl; ++k) {
260 7531432 : const T aik=a[i*A.cl+k];
261 7531432 : const T* bk=b+k*B.cl;
262 24111134 : for(unsigned j=0; j<B.cl; ++j) {
263 16579702 : ci[j]+=aik*bk[j];
264 : }
265 : }
266 : }
267 8572 : }
268 :
269 136364 : template <typename T> void mult( const Matrix<T>& A, const std::vector<T>& B, std::vector<T>& C) {
270 136364 : plumed_assert( A.cl==B.size() );
271 136364 : if( C.size()!=A.rw ) {
272 94 : C.resize(A.rw);
273 : }
274 439992 : for(unsigned i=0; i<A.rw; ++i) {
275 303628 : C[i]= static_cast<T>( 0 );
276 : }
277 439992 : for(unsigned i=0; i<A.rw; ++i)
278 1484656 : for(unsigned k=0; k<A.cl; ++k) {
279 1181028 : C[i]+=A(i,k)*B[k] ;
280 : }
281 136364 : }
282 :
283 : template <typename T> void mult( const std::vector<T>& A, const Matrix<T>& B, std::vector<T>& C) {
284 : plumed_assert( B.rw==A.size() );
285 : if( C.size()!=B.cl ) {
286 : C.resize( B.cl );
287 : }
288 : for(unsigned i=0; i<B.cl; ++i) {
289 : C[i]=static_cast<T>( 0 );
290 : }
291 : // loop on the rows of B in the outer loop, so that the inner loop is on contiguous elements
292 : for(unsigned k=0; k<B.rw; ++k)
293 : for(unsigned i=0; i<B.cl; ++i) {
294 : C[i]+=A[k]*B(k,i);
295 : }
296 : }
297 :
298 766 : template <typename T> void transpose( const Matrix<T>& A, Matrix<T>& AT ) {
299 766 : if( A.rw!=AT.cl || A.cl!=AT.rw ) {
300 36 : AT.resize( A.cl, A.rw );
301 : }
302 3482 : for(unsigned i=0; i<A.cl; ++i)
303 70854 : for(unsigned j=0; j<A.rw; ++j) {
304 68138 : AT(i,j)=A(j,i);
305 : }
306 766 : }
307 :
308 : template <typename T> Log& operator<<(Log& ostr, const Matrix<T>& mat) {
309 : for(unsigned i=0; i<mat.sz; ++i) {
310 : ostr<<mat.data[i]<<" ";
311 : }
312 : return ostr;
313 : }
314 :
315 33 : template <typename T> void matrixOut( Log& ostr, const Matrix<T>& mat) {
316 111 : for(unsigned i=0; i<mat.rw; ++i) {
317 270 : for(unsigned j=0; j<mat.cl; ++j) {
318 192 : ostr<<mat(i,j)<<" ";
319 : }
320 78 : ostr<<"\n";
321 : }
322 33 : return;
323 : }
324 :
325 : /// Cholesky factorization A=L*L^T of a symmetric positive definite n x n matrix.
326 : /// Only the upper triangle of A is read (as done by diagMat). L is returned as a lower
327 : /// triangular matrix stored row by row in a vector of size n*n.
328 : /// Returns false if A is not positive definite (or contains NaNs).
329 15606 : template <typename T> bool choleskyFactorization( const Matrix<T>& A, std::vector<T>& L ) {
330 : const unsigned n=A.nrows();
331 15606 : plumed_assert( A.ncols()==n );
332 : // Right-looking factorization A=U^T*U working on the rows of the upper triangle,
333 : // so that the innermost loop is a vectorizable update of contiguous elements
334 15606 : std::vector<T> U(static_cast<std::size_t>(n)*n,static_cast<T>(0));
335 43170 : for(unsigned j=0; j<n; ++j)
336 67087 : for(unsigned i=j; i<n; ++i) {
337 39523 : U[static_cast<std::size_t>(j)*n+i]=A(j,i);
338 : }
339 43170 : for(unsigned j=0; j<n; ++j) {
340 27564 : T* Uj=U.data()+static_cast<std::size_t>(j)*n;
341 27564 : if( !(Uj[j]>0) ) {
342 : return false;
343 : }
344 27564 : Uj[j]=std::sqrt(Uj[j]);
345 27564 : const T invd=static_cast<T>(1)/Uj[j];
346 39523 : for(unsigned i=j+1; i<n; ++i) {
347 11959 : Uj[i]*=invd;
348 : }
349 39523 : for(unsigned i=j+1; i<n; ++i) {
350 11959 : const T uji=Uj[i];
351 11959 : T* Ui=U.data()+static_cast<std::size_t>(i)*n;
352 23919 : for(unsigned k=i; k<n; ++k) {
353 11960 : Ui[k]-=uji*Uj[k];
354 : }
355 : }
356 : }
357 15606 : L.assign(static_cast<std::size_t>(n)*n,static_cast<T>(0));
358 43170 : for(unsigned i=0; i<n; ++i)
359 67087 : for(unsigned j=0; j<=i; ++j) {
360 39523 : L[static_cast<std::size_t>(i)*n+j]=U[static_cast<std::size_t>(j)*n+i];
361 : }
362 : return true;
363 : }
364 :
365 : /// Inverse of a symmetric positive definite matrix from its Cholesky factorization.
366 : /// Returns false (and leaves inverse untouched) if A is not positive definite.
367 15606 : template <typename T> bool choleskyInvert( const Matrix<T>& A, Matrix<T>& inverse ) {
368 : const unsigned n=A.nrows();
369 15606 : std::vector<T> L;
370 15606 : if( !choleskyFactorization( A, L ) ) {
371 : return false;
372 : }
373 : // Compute M=L^{-1} (lower triangular) row by row: M_i = (e_i - sum_{k<i} L_ik M_k) / L_ii
374 15606 : std::vector<T> M(static_cast<std::size_t>(n)*n,static_cast<T>(0));
375 43170 : for(unsigned i=0; i<n; ++i) {
376 27564 : const T* Li=L.data()+static_cast<std::size_t>(i)*n;
377 27564 : T* Mi=M.data()+static_cast<std::size_t>(i)*n;
378 27564 : Mi[i]=static_cast<T>(1);
379 39523 : for(unsigned k=0; k<i; ++k) {
380 11959 : const T lik=Li[k];
381 11959 : const T* Mk=M.data()+static_cast<std::size_t>(k)*n;
382 23919 : for(unsigned j=0; j<=k; ++j) {
383 11960 : Mi[j]-=lik*Mk[j];
384 : }
385 : }
386 27564 : const T invd=static_cast<T>(1)/Li[i];
387 67087 : for(unsigned j=0; j<=i; ++j) {
388 39523 : Mi[j]*=invd;
389 : }
390 : }
391 : // A^{-1} = M^T M = sum_k (row k of M)^T (row k of M); compute the lower triangle
392 15606 : if( inverse.nrows()!=n || inverse.ncols()!=n ) {
393 0 : inverse.resize(n,n);
394 : }
395 : inverse=static_cast<T>(0);
396 43170 : for(unsigned k=0; k<n; ++k) {
397 27564 : const T* Mk=M.data()+static_cast<std::size_t>(k)*n;
398 67087 : for(unsigned p=0; p<=k; ++p) {
399 39523 : const T mkp=Mk[p];
400 91006 : for(unsigned q=0; q<=p; ++q) {
401 51483 : inverse(p,q)+=mkp*Mk[q];
402 : }
403 : }
404 : }
405 43170 : for(unsigned p=0; p<n; ++p)
406 39523 : for(unsigned q=0; q<p; ++q) {
407 11959 : inverse(q,p)=inverse(p,q);
408 : }
409 : return true;
410 : }
411 :
412 7146 : template <typename T> int diagMatRange( const Matrix<T>& A, unsigned first, unsigned last, std::vector<T>& eigenvals, Matrix<T>& eigenvecs ) {
413 :
414 : // Check matrix is square and symmetric
415 7146 : plumed_assert( A.rw==A.cl );
416 7146 : plumed_assert( A.isSymmetric()==1 );
417 7146 : if( eigenvals.size()!=A.cl ) {
418 9 : eigenvals.resize( A.cl );
419 : }
420 7146 : if( eigenvecs.rw!=A.rw || eigenvecs.cl!=A.cl ) {
421 9 : eigenvecs.resize( A.rw, A.cl );
422 : }
423 7146 : if( A.cl==0 ) {
424 : return 0;
425 : }
426 7146 : plumed_assert( first<=last && last<A.cl );
427 7146 : std::vector<T> da(A.sz);
428 : unsigned k=0;
429 7146 : std::vector<T> evals(A.cl);
430 : // Transfer the matrix to the local array
431 24969 : for (unsigned i=0; i<A.cl; ++i)
432 440076 : for (unsigned j=0; j<A.rw; ++j) {
433 422253 : da[k++]=static_cast<T>( A(j,i) );
434 : }
435 :
436 7146 : int n=A.cl;
437 7146 : int lwork=-1;
438 7146 : int liwork=-1;
439 : int m;
440 7146 : int info=0;
441 : #ifdef __PLUMED_HAS_EXTERNAL_LAPACK
442 : const bool partial=true;
443 : #else
444 : // The internal dsyevr reads out of bounds when only a subset of the eigenpairs is
445 : // requested, so all of them are computed and then the requested ones are copied
446 : const bool partial=false;
447 : #endif
448 : // dsyevr numbers the eigenvalues from 1, in ascending order
449 : const unsigned offset=( partial ? first : 0 );
450 7146 : int il=offset+1;
451 7146 : int iu=( partial ? last+1 : A.cl );
452 : const unsigned nev=iu-il+1;
453 7146 : std::vector<T> work(A.cl);
454 7146 : std::vector<int> iwork(A.cl);
455 7146 : T vl, vu, abstol=0.0;
456 7146 : std::vector<int> isup(2*nev);
457 7146 : std::vector<T> evecs(static_cast<std::size_t>(A.cl)*nev);
458 : if constexpr (std::is_same_v<T,double>) {
459 7146 : plumed_lapack_dsyevr("V", "I", "U", &n, da.data(), &n, &vl, &vu, &il, &iu,
460 : &abstol, &m, evals.data(), evecs.data(), &n,
461 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
462 : } else if constexpr (std::is_same_v<T,float>) {
463 : plumed_lapack_ssyevr("V", "I", "U", &n, da.data(), &n, &vl, &vu, &il, &iu,
464 : &abstol, &m, evals.data(), evecs.data(), &n,
465 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
466 : } else {
467 : static_assert(std::is_same_v<T,double> || std::is_same_v<T,float>, "diagMat is only compatible with float and double");
468 : }
469 7146 : if (info!=0) {
470 : return info;
471 : }
472 :
473 : // Retrieve correct sizes for work and iwork then reallocate
474 7146 : liwork=iwork[0];
475 7146 : iwork.resize(liwork);
476 7146 : lwork=static_cast<int>( work[0] );
477 7146 : work.resize(lwork);
478 :
479 : if constexpr (std::is_same_v<T,double>) {
480 7146 : plumed_lapack_dsyevr("V", "I", "U", &n, da.data(), &n, &vl, &vu, &il, &iu,
481 : &abstol, &m, evals.data(), evecs.data(), &n,
482 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
483 : } else if constexpr (std::is_same_v<T,float>) {
484 : plumed_lapack_ssyevr("V", "I", "U", &n, da.data(), &n, &vl, &vu, &il, &iu,
485 : &abstol, &m, evals.data(), evecs.data(), &n,
486 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
487 : } else {
488 : static_assert(std::is_same_v<T,double> || std::is_same_v<T,float>, "diagMat is only compatible with float and double");
489 : }
490 7146 : if (info!=0) {
491 : return info;
492 : }
493 :
494 7146 : k=(first-offset)*A.rw;
495 23653 : for(unsigned i=first; i<=last; ++i) {
496 16507 : eigenvals[i]=evals[i-offset];
497 : // N.B. For ease of producing projectors we store the eigenvectors
498 : // ROW-WISE in the eigenvectors matrix. The first index is the
499 : // eigenvector number and the second the component
500 79211 : for(unsigned j=0; j<A.rw; ++j) {
501 62704 : eigenvecs(i,j)=evecs[k++];
502 : }
503 : }
504 :
505 : // This changes eigenvectors so that the first non-null element
506 : // of each of them is positive
507 : // We can do it because the phase is arbitrary, and helps making
508 : // the result reproducible
509 23653 : for(unsigned i=first; i<=last; ++i) {
510 : int j;
511 16561 : for(j=0; j<n; j++)
512 16561 : if(eigenvecs(i,j)*eigenvecs(i,j)>1e-14) {
513 : break;
514 : }
515 16507 : if(j<n)
516 16507 : if(eigenvecs(i,j)<0.0)
517 22853 : for(j=0; j<n; j++) {
518 19682 : eigenvecs(i,j)*=-1;
519 : }
520 : }
521 : return 0;
522 0 : }
523 :
524 7127 : template <typename T> int diagMat( const Matrix<T>& A, std::vector<T>& eigenvals, Matrix<T>& eigenvecs ) {
525 7127 : if( A.cl==0 ) {
526 0 : return diagMatRange( A, 0, 0, eigenvals, eigenvecs );
527 : }
528 7127 : return diagMatRange( A, 0, A.cl-1, eigenvals, eigenvecs );
529 : }
530 :
531 1 : template <typename T> int pseudoInvert( const Matrix<T>& A, Matrix<T>& pseudoinverse ) {
532 1 : const unsigned nrows=A.rw, ncols=A.cl;
533 1 : const unsigned nsv=std::min(nrows,ncols);
534 1 : if( pseudoinverse.rw!=ncols || pseudoinverse.cl!=nrows ) {
535 0 : pseudoinverse.resize( ncols, nrows );
536 : }
537 : pseudoinverse=static_cast<T>( 0 );
538 1 : if( nsv==0 ) {
539 : return 0;
540 : }
541 :
542 1 : std::vector<T> da(A.sz);
543 : unsigned k=0;
544 : // Transfer the matrix to the local array
545 501 : for (unsigned i=0; i<A.cl; ++i)
546 250500 : for (unsigned j=0; j<A.rw; ++j) {
547 250000 : da[k++]=static_cast<T>( A(j,i) );
548 : }
549 :
550 1 : int info=0;
551 1 : int nr=nrows, nc=ncols;
552 :
553 : // Create some containers for stuff from the singular value decomposition A=U*S*VT
554 : // (U is nrows x nrows and VT is ncols x ncols, both stored column by column).
555 : // The thin decomposition (JOBZ='S') would be enough, but the internal dgesdd
556 : // writes out of bounds with JOBZ='S', so the full one is used.
557 1 : std::vector<T> S(nsv);
558 1 : std::vector<T> U(static_cast<std::size_t>(nrows)*nrows);
559 1 : std::vector<T> VT(static_cast<std::size_t>(ncols)*ncols);
560 1 : std::vector<int> iwork(8*nsv);
561 :
562 : // This optimizes the size of the work array used in lapack singular value decomposition
563 1 : int lwork=-1;
564 1 : std::vector<T> work(1);
565 : if constexpr (std::is_same_v<T,double> ) {
566 1 : plumed_lapack_dgesdd( "A", &nr, &nc, da.data(), &nr, S.data(),
567 : U.data(), &nr, VT.data(), &nc, work.data(),
568 : &lwork, iwork.data(), &info );
569 : } else if constexpr (std::is_same_v<T,float>) {
570 : plumed_lapack_sgesdd( "A", &nr, &nc, da.data(), &nr, S.data(),
571 : U.data(), &nr, VT.data(), &nc, work.data(),
572 : &lwork, iwork.data(), &info );
573 : } else {
574 : static_assert(std::is_same_v<T,double> || std::is_same_v<T,float>, "pseudoInvert works only for float and double");
575 : }
576 1 : if(info!=0) {
577 : return info;
578 : }
579 :
580 : // Retrieve correct sizes for work and rellocate
581 1 : lwork=(int) work[0];
582 1 : work.resize(lwork);
583 :
584 : // This does the singular value decomposition
585 : if constexpr (std::is_same_v<T,double> ) {
586 1 : plumed_lapack_dgesdd( "A", &nr, &nc, da.data(), &nr, S.data(),
587 : U.data(), &nr, VT.data(), &nc, work.data(),
588 : &lwork, iwork.data(), &info );
589 : } else if constexpr (std::is_same_v<T,float>) {
590 : plumed_lapack_sgesdd( "A", &nr, &nc, da.data(), &nr, S.data(),
591 : U.data(), &nr, VT.data(), &nc, work.data(),
592 : &lwork, iwork.data(), &info );
593 : } else {
594 : static_assert(std::is_same_v<T,double> || std::is_same_v<T,float>, "pseudoInvert works only for float and double");
595 : }
596 1 : if(info!=0) {
597 : return info;
598 : }
599 :
600 : // Compute the tolerance on the singular values ( machine epsilon * number of singular values * maximum singular value )
601 : T tol;
602 1 : tol=S[0];
603 500 : for(unsigned i=1; i<nsv; ++i) {
604 499 : if( S[i]>tol ) {
605 : tol=S[i];
606 : }
607 : }
608 1 : tol*=nsv*epsilon;
609 :
610 : // pseudoinverse = V * S^{-1} * U^T = sum_k (column k of V) (column k of U)^T / S_k,
611 : // summing only over the singular values larger than the tolerance
612 501 : for(unsigned isv=0; isv<nsv; ++isv) {
613 500 : if( !(S[isv]>tol) ) {
614 1 : continue;
615 : }
616 499 : const T inv=static_cast<T>(1)/S[isv];
617 499 : const T* Uk=U.data()+static_cast<std::size_t>(isv)*nrows;
618 249999 : for(unsigned i=0; i<ncols; ++i) {
619 249500 : const T a=inv*VT[isv+static_cast<std::size_t>(i)*ncols];
620 124999500 : for(unsigned j=0; j<nrows; ++j) {
621 124750000 : pseudoinverse(i,j)+=a*Uk[j];
622 : }
623 : }
624 : }
625 :
626 : return 0;
627 : }
628 :
629 149887 : template <typename T> int Invert( const Matrix<T>& A, Matrix<T>& inverse ) {
630 :
631 149887 : if( A.isSymmetric()==1 ) {
632 : // Symmetric positive definite matrices (e.g. covariance and metric matrices) are inverted
633 : // with a Cholesky factorization, which is several times faster than diagonalizing them.
634 15606 : if( choleskyInvert( A, inverse ) ) {
635 15606 : return 0;
636 : }
637 : // GAT -- I only ever use symmetric matrices so I can invert them like this.
638 : // I choose to do this as I have had problems with the more general way of doing this that
639 : // is implemented below.
640 0 : std::vector<T> eval(A.rw);
641 0 : Matrix<T> evec(A.rw,A.cl), tevec(A.rw,A.cl);
642 : int err;
643 0 : err=diagMat( A, eval, evec );
644 0 : if(err!=0) {
645 : return err;
646 : }
647 0 : for (unsigned i=0; i<A.rw; ++i)
648 0 : for (unsigned j=0; j<A.cl; ++j) {
649 0 : tevec(i,j)=evec(j,i)/eval[j];
650 : }
651 0 : mult(tevec,evec,inverse);
652 : } else {
653 134281 : std::vector<T> da(A.sz);
654 134281 : std::vector<int> ipiv(A.cl);
655 : unsigned k=0;
656 134281 : int n=A.rw;
657 134281 : int info=0;
658 402846 : for(unsigned i=0; i<A.cl; ++i)
659 805710 : for(unsigned j=0; j<A.rw; ++j) {
660 537145 : da[k++]=static_cast<T>( A(j,i) );
661 : }
662 : if constexpr (std::is_same_v<T,double> ) {
663 134281 : plumed_lapack_dgetrf(&n,&n,da.data(),&n,ipiv.data(),&info);
664 : } else if constexpr (std::is_same_v<T,float>) {
665 : plumed_lapack_sgetrf(&n,&n,da.data(),&n,ipiv.data(),&info);
666 : } else {
667 : static_assert(std::is_same_v<T,double> || std::is_same_v<T,float>, "Invert works only for float and double");
668 : }
669 134281 : if(info!=0) {
670 : return info;
671 : }
672 :
673 134281 : int lwork=-1;
674 134281 : std::vector<T> work(A.cl);
675 : if constexpr (std::is_same_v<T,double> ) {
676 134281 : plumed_lapack_dgetri(&n,da.data(),&n,ipiv.data(),work.data(),&lwork,&info);
677 : } else if constexpr (std::is_same_v<T,float>) {
678 : plumed_lapack_sgetri(&n,da.data(),&n,ipiv.data(),work.data(),&lwork,&info);
679 : } else {
680 : static_assert(std::is_same_v<T,double> || std::is_same_v<T,float>, "Invert works only for float and double");
681 : }
682 134281 : if(info!=0) {
683 : return info;
684 : }
685 :
686 134281 : lwork=static_cast<int>( work[0] );
687 134281 : work.resize(lwork);
688 : if constexpr (std::is_same_v<T,double> ) {
689 134281 : plumed_lapack_dgetri(&n,da.data(),&n,ipiv.data(),work.data(),&lwork,&info);
690 : } else if constexpr (std::is_same_v<T,float>) {
691 : plumed_lapack_sgetri(&n,da.data(),&n,ipiv.data(),work.data(),&lwork,&info);
692 : } else {
693 : static_assert(std::is_same_v<T,double> || std::is_same_v<T,float>, "Invert works only for float and double");
694 : }
695 134281 : if(info!=0) {
696 : return info;
697 : }
698 :
699 134281 : if( inverse.cl!=A.cl || inverse.rw!=A.rw ) {
700 0 : inverse.resize(A.rw,A.cl);
701 : }
702 : k=0;
703 402846 : for(unsigned i=0; i<A.rw; ++i)
704 805710 : for(unsigned j=0; j<A.cl; ++j) {
705 537145 : inverse(j,i)=da[k++];
706 : }
707 0 : }
708 :
709 : return 0;
710 : }
711 :
712 3655 : template <typename T> void cholesky( const Matrix<T>& A, Matrix<T>& B ) {
713 :
714 3655 : plumed_assert( A.rw==A.cl && A.isSymmetric() );
715 3655 : Matrix<T> L(A.rw,A.cl);
716 : L=0.;
717 3655 : std::vector<T> D(A.rw,0.);
718 8969 : for(unsigned i=0; i<A.rw; ++i) {
719 5314 : L(i,i)=static_cast<T>( 1 );
720 6973 : for (unsigned j=0; j<i; ++j) {
721 1659 : L(i,j)=A(i,j);
722 1659 : for (unsigned k=0; k<j; ++k) {
723 0 : L(i,j)-=L(i,k)*L(j,k)*D[k];
724 : }
725 1659 : if (D[j]!=0.) {
726 1659 : L(i,j)/=D[j];
727 : } else {
728 0 : L(i,j)=static_cast<T>( 0 );
729 : }
730 : }
731 5314 : D[i]=A(i,i);
732 6973 : for (unsigned k=0; k<i; ++k) {
733 1659 : D[i]-=L(i,k)*L(i,k)*D[k];
734 : }
735 : }
736 :
737 8969 : for(unsigned i=0; i<A.rw; ++i) {
738 5314 : D[i]=(D[i]>0.?std::sqrt(D[i]):0.);
739 : }
740 3655 : if( B.rw!=A.rw || B.cl!=A.cl ) {
741 0 : B.resize( A.rw, A.cl);
742 : }
743 : B=0.;
744 8969 : for(unsigned i=0; i<A.rw; ++i)
745 12287 : for(unsigned j=0; j<=i; ++j) {
746 6973 : B(i,j)+=L(i,j)*D[j];
747 : }
748 3655 : }
749 :
750 : template <typename T> void chol_elsolve( const Matrix<T>& M, const std::vector<T>& b, std::vector<T>& y ) {
751 :
752 : plumed_assert( M.rw==M.cl && M(0,1)==0.0 && b.size()==M.rw );
753 : if( y.size()!=M.rw ) {
754 : y.resize( M.rw );
755 : }
756 : for(unsigned i=0; i<M.rw; ++i) {
757 : y[i]=b[i];
758 : for(unsigned j=0; j<i; ++j) {
759 : y[i]-=M(i,j)*y[j];
760 : }
761 : y[i]*=1.0/M(i,i);
762 : }
763 : }
764 :
765 0 : template <typename T> int logdet( const Matrix<T>& M, T& ldet ) {
766 : // Check matrix is square and symmetric
767 0 : plumed_assert( M.rw==M.cl || M.isSymmetric() );
768 :
769 : // For a positive definite matrix, log(det(M)) = 2 sum_i log(L_ii) with M=L*L^T
770 : {
771 0 : std::vector<T> L;
772 0 : if( choleskyFactorization( M, L ) ) {
773 0 : ldet=0;
774 0 : for(unsigned i=0; i<M.rw; ++i) {
775 0 : ldet+=2*std::log(L[static_cast<std::size_t>(i)*M.rw+i]);
776 : }
777 : return 0;
778 : }
779 : }
780 :
781 0 : std::vector<T> da(M.sz);
782 : unsigned k=0;
783 0 : std::vector<T> evals(M.cl);
784 : // Transfer the matrix to the local array
785 0 : for (unsigned i=0; i<M.rw; ++i)
786 0 : for (unsigned j=0; j<M.cl; ++j) {
787 0 : da[k++]=static_cast<T>( M(j,i) );
788 : }
789 :
790 0 : int n=M.cl;
791 0 : int lwork=-1;
792 0 : int liwork=-1;
793 0 : int info=0;
794 : int m;
795 0 : int one=1;
796 0 : std::vector<T> work(M.rw);
797 0 : std::vector<int> iwork(M.rw);
798 0 : T vl, vu, abstol=0.0;
799 0 : std::vector<int> isup(2*M.rw);
800 0 : std::vector<T> evecs(M.sz);
801 : if constexpr (std::is_same_v<T,double> ) {
802 0 : plumed_lapack_dsyevr("N", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
803 : &abstol, &m, evals.data(), evecs.data(), &n,
804 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
805 : } else if constexpr (std::is_same_v<T,float>) {
806 : plumed_lapack_ssyevr("N", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
807 : &abstol, &m, evals.data(), evecs.data(), &n,
808 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
809 : } else {
810 : static_assert(std::is_same_v<T,double> || std::is_same_v<T,float>, "logdet works only for float and double");
811 : }
812 0 : if (info!=0) {
813 : return info;
814 : }
815 :
816 : // Retrieve correct sizes for work and iwork then reallocate
817 0 : lwork=static_cast<int>( work[0] );
818 0 : work.resize(lwork);
819 0 : liwork=iwork[0];
820 0 : iwork.resize(liwork);
821 : if constexpr (std::is_same_v<T,double> ) {
822 0 : plumed_lapack_dsyevr("N", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
823 : &abstol, &m, evals.data(), evecs.data(), &n,
824 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
825 : } else if constexpr (std::is_same_v<T,float>) {
826 : plumed_lapack_ssyevr("N", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
827 : &abstol, &m, evals.data(), evecs.data(), &n,
828 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
829 : } else {
830 : static_assert(std::is_same_v<T,double> || std::is_same_v<T,float>, "logdet works only for float and double");
831 : }
832 0 : if (info!=0) {
833 : return info;
834 : }
835 :
836 : // Transfer the eigenvalues and eigenvectors to the output
837 0 : ldet=0;
838 0 : for(unsigned i=0; i<M.cl; i++) {
839 0 : ldet+=log(evals[i]);
840 : }
841 :
842 : return 0;
843 : }
844 :
845 :
846 :
847 : }
848 : #endif
|
