1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2011-2021 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
dotProduct(const std::vector<T> & A,const std::vector<T> & B)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; for(unsigned i=0; i<A.size(); ++i) { val+=A[i]*B[i]; }
40 return val;
41 }
42
43 /// Calculate the dot product between a vector and itself
norm(const std::vector<T> & A)44 template <typename T> T norm( const std::vector<T>& A ) {
45 T val; for(unsigned i=0; i<A.size(); ++i) { val+=A[i]*A[i]; }
46 return val;
47 }
48
49 /// This class stores a full matrix and allows one to do some simple matrix operations
50 template <typename T>
51 class Matrix:
52 public MatrixSquareBracketsAccess<Matrix<T>,T>
53 {
54 /// Multiply matrix by scalar
55 template <typename U> friend Matrix<U> operator*(U&, const Matrix<U>& );
56 /// Matrix matrix multiply
57 template <typename U> friend void mult( const Matrix<U>&, const Matrix<U>&, Matrix<U>& );
58 /// Matrix times a std::vector
59 template <typename U> friend void mult( const Matrix<U>&, const std::vector<U>&, std::vector<U>& );
60 /// std::vector times a Matrix
61 template <typename U> friend void mult( const std::vector<U>&, const Matrix<U>&, std::vector<U>& );
62 /// Matrix transpose
63 template <typename U> friend void transpose( const Matrix<U>&, Matrix<U>& );
64 /// Output the entire matrix on a single line
65 template <typename U> friend Log& operator<<(Log&, const Matrix<U>& );
66 /// Output the Matrix in matrix form
67 template <typename U> friend void matrixOut( Log&, const Matrix<U>& );
68 /// Diagonalize a symmetric matrix - returns zero if diagonalization worked
69 template <typename U> friend int diagMat( const Matrix<U>&, std::vector<double>&, Matrix<double>& );
70 /// Calculate the Moore-Penrose Pseudoinverse of a matrix
71 template <typename U> friend int pseudoInvert( const Matrix<U>&, Matrix<double>& );
72 /// Calculate the logarithm of the determinant of a symmetric matrix - returns zero if succesfull
73 template <typename U> friend int logdet( const Matrix<U>&, double& );
74 /// Invert a matrix (works for both symmetric and asymmetric matrices) - returns zero if sucesfull
75 template <typename U> friend int Invert( const Matrix<U>&, Matrix<double>& );
76 /// Do a cholesky decomposition of a matrix
77 template <typename U> friend void cholesky( const Matrix<U>&, Matrix<U>& );
78 /// Solve a system of equations using the cholesky decomposition
79 template <typename U> friend void chol_elsolve( const Matrix<U>&, const std::vector<U>&, std::vector<U>& );
80 private:
81 /// Number of elements in matrix (nrows*ncols)
82 unsigned sz;
83 /// Number of rows in matrix
84 unsigned rw;
85 /// Number of columns in matrix
86 unsigned cl;
87 /// The data in the matrix
88 std::vector<T> data;
89 public:
sz(nr * nc)90 Matrix(const unsigned nr=0, const unsigned nc=0 ) : sz(nr*nc), rw(nr), cl(nc), data(nr*nc) {}
Matrix(const Matrix<T> & t)91 Matrix(const Matrix<T>& t) : sz(t.sz), rw(t.rw), cl(t.cl), data(t.data) {}
92 /// Resize the matrix
resize(const unsigned nr,const unsigned nc)93 void resize( const unsigned nr, const unsigned nc ) { rw=nr; cl=nc; sz=nr*nc; data.resize(sz); }
94 /// Return the number of rows
nrows()95 inline unsigned nrows() const { return rw; }
96 /// Return the number of columns
ncols()97 inline unsigned ncols() const { return cl; }
98 /// Return the contents of the matrix as a vector of length rw*cl
getVector()99 inline std::vector<T>& getVector() { return data; }
100 /// Set the matrix from a vector input
setFromVector(const std::vector<T> & vecin)101 inline void setFromVector( const std::vector<T>& vecin ) { plumed_assert( vecin.size()==sz ); for(unsigned i=0; i<sz; ++i) data[i]=vecin[i]; }
102 /// Return element i,j of the matrix
operator()103 inline const T& operator () (const unsigned& i, const unsigned& j) const { return data[j+i*cl]; }
104 /// Return a referenre to element i,j of the matrix
operator()105 inline T& operator () (const unsigned& i, const unsigned& j) { return data[j+i*cl]; }
106 /// Set all elements of the matrix equal to the value of v
107 Matrix<T>& operator=(const T& v) {
108 for(unsigned i=0; i<sz; ++i) { data[i]=v; }
109 return *this;
110 }
111 /// Set the Matrix equal to another Matrix
112 Matrix<T>& operator=(const Matrix<T>& m) {
113 sz=m.sz;
114 rw=m.rw;
115 cl=m.cl;
116 data=m.data;
117 return *this;
118 }
119 /// Set the Matrix equal to the value of a standard vector - used for readin
120 Matrix<T>& operator=(const std::vector<T>& v) {
121 plumed_dbg_assert( v.size()==sz );
122 for(unsigned i=0; i<sz; ++i) { data[i]=v[i]; }
123 return *this;
124 }
125 /// Add v to all elements of the Matrix
126 Matrix<T> operator+=(const T& v) {
127 for(unsigned i=0; i<sz; ++i) { data[i]+=v; }
128 return *this;
129 }
130 /// Multiply all elements by v
131 Matrix<T> operator*=(const T& v) {
132 for(unsigned i=0; i<sz; ++i) { data[i]*=v; }
133 return *this;
134 }
135 /// Matrix addition
136 Matrix<T>& operator+=(const Matrix<T>& m) {
137 plumed_dbg_assert( m.rw==rw && m.cl==cl );
138 data+=m.data;
139 return *this;
140 }
141 /// Subtract v from all elements of the Matrix
142 Matrix<T> operator-=(const T& v) {
143 for(unsigned i=0; i<sz; ++i) { data-=v; }
144 return *this;
145 }
146 /// Matrix subtraction
147 Matrix<T>& operator-=(const Matrix<T>& m) {
148 plumed_dbg_assert( m.rw==rw && m.cl==cl );
149 data-=m.data;
150 return *this;
151 }
152 /// Test if the matrix is symmetric or not
isSymmetric()153 unsigned isSymmetric() const {
154 if (rw!=cl) { return 0; }
155 unsigned sym=1;
156 for(unsigned i=1; i<rw; ++i) for(unsigned j=0; j<i; ++j) if( std::fabs(data[i+j*cl]-data[j+i*cl])>1.e-10 ) { sym=0; break; }
157 return sym;
158 }
159 };
160
161 /// Multiply matrix by scalar
162 template <typename T> Matrix<T> operator*(T& v, const Matrix<T>& m ) {
163 Matrix<T> new_m(m);
164 new_m*=v;
165 return new_m;
166 }
167
mult(const Matrix<T> & A,const Matrix<T> & B,Matrix<T> & C)168 template <typename T> void mult( const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& C ) {
169 plumed_assert(A.cl==B.rw);
170 if( A.rw !=C.rw || B.cl !=C.cl ) { C.resize( A.rw, B.cl ); } C=static_cast<T>( 0 );
171 for(unsigned i=0; i<A.rw; ++i) for(unsigned j=0; j<B.cl; ++j) for (unsigned k=0; k<A.cl; ++k) C(i,j)+=A(i,k)*B(k,j);
172 }
173
mult(const Matrix<T> & A,const std::vector<T> & B,std::vector<T> & C)174 template <typename T> void mult( const Matrix<T>& A, const std::vector<T>& B, std::vector<T>& C) {
175 plumed_assert( A.cl==B.size() );
176 if( C.size()!=A.rw ) { C.resize(A.rw); }
177 for(unsigned i=0; i<A.rw; ++i) { C[i]= static_cast<T>( 0 ); }
178 for(unsigned i=0; i<A.rw; ++i) for(unsigned k=0; k<A.cl; ++k) C[i]+=A(i,k)*B[k] ;
179 }
180
mult(const std::vector<T> & A,const Matrix<T> & B,std::vector<T> & C)181 template <typename T> void mult( const std::vector<T>& A, const Matrix<T>& B, std::vector<T>& C) {
182 plumed_assert( B.rw==A.size() );
183 if( C.size()!=B.cl ) {C.resize( B.cl );}
184 for(unsigned i=0; i<B.cl; ++i) { C[i]=static_cast<T>( 0 ); }
185 for(unsigned i=0; i<B.cl; ++i) for(unsigned k=0; k<B.rw; ++k) C[i]+=A[k]*B(k,i);
186 }
187
transpose(const Matrix<T> & A,Matrix<T> & AT)188 template <typename T> void transpose( const Matrix<T>& A, Matrix<T>& AT ) {
189 if( A.rw!=AT.cl || A.cl!=AT.rw ) AT.resize( A.cl, A.rw );
190 for(unsigned i=0; i<A.cl; ++i) for(unsigned j=0; j<A.rw; ++j) AT(i,j)=A(j,i);
191 }
192
193 template <typename T> Log& operator<<(Log& ostr, const Matrix<T>& mat) {
194 for(unsigned i=0; i<mat.sz; ++i) ostr<<mat.data[i]<<" ";
195 return ostr;
196 }
197
matrixOut(Log & ostr,const Matrix<T> & mat)198 template <typename T> void matrixOut( Log& ostr, const Matrix<T>& mat) {
199 for(unsigned i=0; i<mat.rw; ++i) {
200 for(unsigned j=0; j<mat.cl; ++j) { ostr<<mat(i,j)<<" "; }
201 ostr<<"\n";
202 }
203 return;
204 }
205
diagMat(const Matrix<T> & A,std::vector<double> & eigenvals,Matrix<double> & eigenvecs)206 template <typename T> int diagMat( const Matrix<T>& A, std::vector<double>& eigenvals, Matrix<double>& eigenvecs ) {
207
208 // Check matrix is square and symmetric
209 plumed_assert( A.rw==A.cl ); plumed_assert( A.isSymmetric()==1 );
210 std::vector<double> da(A.sz);
211 unsigned k=0;
212 std::vector<double> evals(A.cl);
213 // Transfer the matrix to the local array
214 for (unsigned i=0; i<A.cl; ++i) for (unsigned j=0; j<A.rw; ++j) da[k++]=static_cast<double>( A(j,i) );
215
216 int n=A.cl; int lwork=-1, liwork=-1, m, info, one=1;
217 std::vector<double> work(A.cl);
218 std::vector<int> iwork(A.cl);
219 double vl, vu, abstol=0.0;
220 std::vector<int> isup(2*A.cl);
221 std::vector<double> evecs(A.sz);
222
223 plumed_lapack_dsyevr("V", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
224 &abstol, &m, evals.data(), evecs.data(), &n,
225 isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
226 if (info!=0) return info;
227
228 // Retrieve correct sizes for work and iwork then reallocate
229 liwork=iwork[0]; iwork.resize(liwork);
230 lwork=static_cast<int>( work[0] ); work.resize(lwork);
231
232 plumed_lapack_dsyevr("V", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
233 &abstol, &m, evals.data(), evecs.data(), &n,
234 isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
235 if (info!=0) return info;
236
237 if( eigenvals.size()!=A.cl ) { eigenvals.resize( A.cl ); }
238 if( eigenvecs.rw!=A.rw || eigenvecs.cl!=A.cl ) { eigenvecs.resize( A.rw, A.cl ); }
239 k=0;
240 for(unsigned i=0; i<A.cl; ++i) {
241 eigenvals[i]=evals[i];
242 // N.B. For ease of producing projectors we store the eigenvectors
243 // ROW-WISE in the eigenvectors matrix. The first index is the
244 // eigenvector number and the second the component
245 for(unsigned j=0; j<A.rw; ++j) { eigenvecs(i,j)=evecs[k++]; }
246 }
247
248 // This changes eigenvectors so that the first non-null element
249 // of each of them is positive
250 // We can do it because the phase is arbitrary, and helps making
251 // the result reproducible
252 for(int i=0; i<n; ++i) {
253 int j;
254 for(j=0; j<n; j++) if(eigenvecs(i,j)*eigenvecs(i,j)>1e-14) break;
255 if(j<n) if(eigenvecs(i,j)<0.0) for(j=0; j<n; j++) eigenvecs(i,j)*=-1;
256 }
257 return 0;
258 }
259
pseudoInvert(const Matrix<T> & A,Matrix<double> & pseudoinverse)260 template <typename T> int pseudoInvert( const Matrix<T>& A, Matrix<double>& pseudoinverse ) {
261 std::vector<double> da(A.sz);
262 unsigned k=0;
263 // Transfer the matrix to the local array
264 for (unsigned i=0; i<A.cl; ++i) for (unsigned j=0; j<A.rw; ++j) da[k++]=static_cast<double>( A(j,i) );
265
266 int nsv, info, nrows=A.rw, ncols=A.cl;
267 if(A.rw>A.cl) {nsv=A.cl;} else {nsv=A.rw;}
268
269 // Create some containers for stuff from single value decomposition
270 std::vector<double> S(nsv);
271 std::vector<double> U(nrows*nrows);
272 std::vector<double> VT(ncols*ncols);
273 std::vector<int> iwork(8*nsv);
274
275 // This optimizes the size of the work array used in lapack singular value decomposition
276 int lwork=-1;
277 std::vector<double> work(1);
278 plumed_lapack_dgesdd( "A", &nrows, &ncols, da.data(), &nrows, S.data(), U.data(), &nrows, VT.data(), &ncols, work.data(), &lwork, iwork.data(), &info );
279 if(info!=0) return info;
280
281 // Retrieve correct sizes for work and rellocate
282 lwork=(int) work[0]; work.resize(lwork);
283
284 // This does the singular value decomposition
285 plumed_lapack_dgesdd( "A", &nrows, &ncols, da.data(), &nrows, S.data(), U.data(), &nrows, VT.data(), &ncols, work.data(), &lwork, iwork.data(), &info );
286 if(info!=0) return info;
287
288 // Compute the tolerance on the singular values ( machine epsilon * number of singular values * maximum singular value )
289 double tol; tol=S[0]; for(int i=1; i<nsv; ++i) { if( S[i]>tol ) { tol=S[i]; } } tol*=nsv*epsilon;
290
291 // Get the inverses of the singlular values
292 Matrix<double> Si( ncols, nrows ); Si=0.0;
293 for(int i=0; i<nsv; ++i) { if( S[i]>tol ) { Si(i,i)=1./S[i]; } else { Si(i,i)=0.0; } }
294
295 // Now extract matrices for pseudoinverse
296 Matrix<double> V( ncols, ncols ), UT( nrows, nrows ), tmp( ncols, nrows );
297 k=0; for(int i=0; i<nrows; ++i) { for(int j=0; j<nrows; ++j) { UT(i,j)=U[k++]; } }
298 k=0; for(int i=0; i<ncols; ++i) { for(int j=0; j<ncols; ++j) { V(i,j)=VT[k++]; } }
299
300 // And do matrix algebra to construct the pseudoinverse
301 if( pseudoinverse.rw!=ncols || pseudoinverse.cl!=nrows ) pseudoinverse.resize( ncols, nrows );
302 mult( V, Si, tmp ); mult( tmp, UT, pseudoinverse );
303
304 return 0;
305 }
306
Invert(const Matrix<T> & A,Matrix<double> & inverse)307 template <typename T> int Invert( const Matrix<T>& A, Matrix<double>& inverse ) {
308
309 if( A.isSymmetric()==1 ) {
310 // GAT -- I only ever use symmetric matrices so I can invert them like this.
311 // I choose to do this as I have had problems with the more general way of doing this that
312 // is implemented below.
313 std::vector<double> eval(A.rw); Matrix<double> evec(A.rw,A.cl), tevec(A.rw,A.cl);
314 int err; err=diagMat( A, eval, evec );
315 if(err!=0) return err;
316 for (unsigned i=0; i<A.rw; ++i) for (unsigned j=0; j<A.cl; ++j) tevec(i,j)=evec(j,i)/eval[j];
317 mult(tevec,evec,inverse);
318 } else {
319 std::vector<double> da(A.sz);
320 std::vector<int> ipiv(A.cl);
321 unsigned k=0; int n=A.rw, info;
322 for(unsigned i=0; i<A.cl; ++i) for(unsigned j=0; j<A.rw; ++j) da[k++]=static_cast<double>( A(j,i) );
323
324 plumed_lapack_dgetrf(&n,&n,da.data(),&n,ipiv.data(),&info);
325 if(info!=0) return info;
326
327 int lwork=-1;
328 std::vector<double> work(A.cl);
329 plumed_lapack_dgetri(&n,da.data(),&n,ipiv.data(),work.data(),&lwork,&info);
330 if(info!=0) return info;
331
332 lwork=static_cast<int>( work[0] ); work.resize(lwork);
333 plumed_lapack_dgetri(&n,da.data(),&n,ipiv.data(),work.data(),&lwork,&info);
334 if(info!=0) return info;
335
336 if( inverse.cl!=A.cl || inverse.rw!=A.rw ) { inverse.resize(A.rw,A.cl); }
337 k=0; for(unsigned i=0; i<A.rw; ++i) for(unsigned j=0; j<A.cl; ++j) inverse(j,i)=da[k++];
338 }
339
340 return 0;
341 }
342
cholesky(const Matrix<T> & A,Matrix<T> & B)343 template <typename T> void cholesky( const Matrix<T>& A, Matrix<T>& B ) {
344
345 plumed_assert( A.rw==A.cl && A.isSymmetric() );
346 Matrix<T> L(A.rw,A.cl); L=0.;
347 std::vector<T> D(A.rw,0.);
348 for(unsigned i=0; i<A.rw; ++i) {
349 L(i,i)=static_cast<T>( 1 );
350 for (unsigned j=0; j<i; ++j) {
351 L(i,j)=A(i,j);
352 for (unsigned k=0; k<j; ++k) L(i,j)-=L(i,k)*L(j,k)*D[k];
353 if (D[j]!=0.) L(i,j)/=D[j]; else L(i,j)=static_cast<T>( 0 );
354 }
355 D[i]=A(i,i);
356 for (unsigned k=0; k<i; ++k) D[i]-=L(i,k)*L(i,k)*D[k];
357 }
358
359 for(unsigned i=0; i<A.rw; ++i) D[i]=(D[i]>0.?sqrt(D[i]):0.);
360 if( B.rw!=A.rw || B.cl!=A.cl ) { B.resize( A.rw, A.cl); }
361 B=0.; for(unsigned i=0; i<A.rw; ++i) for(unsigned j=0; j<=i; ++j) B(i,j)+=L(i,j)*D[j];
362 }
363
chol_elsolve(const Matrix<T> & M,const std::vector<T> & b,std::vector<T> & y)364 template <typename T> void chol_elsolve( const Matrix<T>& M, const std::vector<T>& b, std::vector<T>& y ) {
365
366 plumed_assert( M.rw==M.cl && M(0,1)==0.0 && b.size()==M.rw );
367 if( y.size()!=M.rw ) { y.resize( M.rw ); }
368 for(unsigned i=0; i<M.rw; ++i) {
369 y[i]=b[i];
370 for(unsigned j=0; j<i; ++j) y[i]-=M(i,j)*y[j];
371 y[i]*=1.0/M(i,i);
372 }
373 }
374
logdet(const Matrix<T> & M,double & ldet)375 template <typename T> int logdet( const Matrix<T>& M, double& ldet ) {
376 // Check matrix is square and symmetric
377 plumed_assert( M.rw==M.cl || M.isSymmetric() );
378
379 std::vector<double> da(M.sz);
380 unsigned k=0;
381 std::vector<double> evals(M.cl);
382 // Transfer the matrix to the local array
383 for (unsigned i=0; i<M.rw; ++i) for (unsigned j=0; j<M.cl; ++j) da[k++]=static_cast<double>( M(j,i) );
384
385 int n=M.cl; int lwork=-1, liwork=-1, info, m, one=1;
386 std::vector<double> work(M.rw);
387 std::vector<int> iwork(M.rw);
388 double vl, vu, abstol=0.0;
389 std::vector<int> isup(2*M.rw);
390 std::vector<double> evecs(M.sz);
391 plumed_lapack_dsyevr("N", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
392 &abstol, &m, evals.data(), evecs.data(), &n,
393 isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
394 if (info!=0) return info;
395
396 // Retrieve correct sizes for work and iwork then reallocate
397 lwork=static_cast<int>( work[0] ); work.resize(lwork);
398 liwork=iwork[0]; iwork.resize(liwork);
399
400 plumed_lapack_dsyevr("N", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
401 &abstol, &m, evals.data(), evecs.data(), &n,
402 isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
403 if (info!=0) return info;
404
405 // Transfer the eigenvalues and eigenvectors to the output
406 ldet=0; for(unsigned i=0; i<M.cl; i++) { ldet+=log(evals[i]); }
407
408 return 0;
409 }
410
411
412
413 }
414 #endif
415