Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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; 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
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 3180546 : 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 assymetric 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:
90 3179209 : Matrix(const unsigned nr=0, const unsigned nc=0 ) : sz(nr*nc), rw(nr), cl(nc), data(nr*nc) {}
91 473 : Matrix(const Matrix<T>& t) : sz(t.sz), rw(t.rw), cl(t.cl), data(t.data) {}
92 : /// Resize the matrix
93 552885 : void resize( const unsigned nr, const unsigned nc ) { rw=nr; cl=nc; sz=nr*nc; data.resize(sz); }
94 : /// Return the number of rows
95 2913105 : inline unsigned nrows() const { return rw; }
96 : /// Return the number of columns
97 6307067 : inline unsigned ncols() const { return cl; }
98 : /// Return element i,j of the matrix
99 257707019 : inline T operator () (const unsigned& i, const unsigned& j) const { return data[j+i*cl]; }
100 : /// Return a referenre to element i,j of the matrix
101 86618557 : inline T& operator () (const unsigned& i, const unsigned& j) { return data[j+i*cl]; }
102 : /// Set all elements of the matrix equal to the value of v
103 14695 : Matrix<T>& operator=(const T& v) {
104 14695 : for(unsigned i=0; i<sz; ++i) { data[i]=v; }
105 14695 : return *this;
106 : }
107 : /// Set the Matrix equal to another Matrix
108 4489 : Matrix<T>& operator=(const Matrix<T>& m) {
109 4489 : sz=m.sz;
110 4489 : rw=m.rw;
111 4489 : cl=m.cl;
112 4489 : data=m.data;
113 4489 : return *this;
114 : }
115 : /// Set the Matrix equal to the value of a standard vector - used for readin
116 : Matrix<T>& operator=(const std::vector<T>& v) {
117 : plumed_dbg_assert( v.size()==sz );
118 : for(unsigned i=0; i<sz; ++i) { data[i]=v[i]; }
119 : return *this;
120 : }
121 : /// Add v to all elements of the Matrix
122 : Matrix<T> operator+=(const T& v) {
123 : for(unsigned i=0; i<sz; ++i) { data[i]+=v; }
124 : return *this;
125 : }
126 : /// Multiply all elements by v
127 2 : Matrix<T> operator*=(const T& v) {
128 2 : for(unsigned i=0; i<sz; ++i) { data[i]*=v; }
129 2 : return *this;
130 : }
131 : /// Matrix addition
132 : Matrix<T>& operator+=(const Matrix<T>& m) {
133 : plumed_dbg_assert( m.rw==rw && m.cl==cl );
134 : data+=m.data;
135 : return *this;
136 : }
137 : /// Subtract v from all elements of the Matrix
138 : Matrix<T> operator-=(const T& v) {
139 : for(unsigned i=0; i<sz; ++i) { data-=v; }
140 : return *this;
141 : }
142 : /// Matrix subtraction
143 : Matrix<T>& operator-=(const Matrix<T>& m) {
144 : plumed_dbg_assert( m.rw==rw && m.cl==cl );
145 : data-=m.data;
146 : return *this;
147 : }
148 : /// Test if the matrix is symmetric or not
149 576281 : unsigned isSymmetric() const {
150 576281 : if (rw!=cl) { return 0; }
151 576281 : unsigned sym=1;
152 576281 : 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; }
153 576681 : return sym;
154 : }
155 : };
156 :
157 : /// Multiply matrix by scalar
158 2 : template <typename T> Matrix<T> operator*(T& v, const Matrix<T>& m ) {
159 2 : Matrix<T> new_m(m);
160 2 : new_m*=v;
161 2 : return new_m;
162 : }
163 :
164 13550 : template <typename T> void mult( const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& C ) {
165 13550 : plumed_assert(A.cl==B.rw);
166 13550 : if( A.rw !=C.rw || B.cl !=C.cl ) { C.resize( A.rw, B.cl ); } C=static_cast<T>( 0 );
167 13550 : 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);
168 13550 : }
169 :
170 55 : template <typename T> void mult( const Matrix<T>& A, const std::vector<T>& B, std::vector<T>& C) {
171 55 : plumed_assert( A.cl==B.size() );
172 54 : if( C.size()!=A.rw ) { C.resize(A.rw); }
173 54 : for(unsigned i=0; i<A.rw; ++i) { C[i]= static_cast<T>( 0 ); }
174 55 : for(unsigned i=0; i<A.rw; ++i) for(unsigned k=0; k<A.cl; ++k) C[i]+=A(i,k)*B[k] ;
175 56 : }
176 :
177 : template <typename T> void mult( const std::vector<T>& A, const Matrix<T>& B, std::vector<T>& C) {
178 : plumed_assert( B.rw==A.size() );
179 : if( C.size()!=B.cl ) {C.resize( B.cl );}
180 : for(unsigned i=0; i<B.cl; ++i) { C[i]=static_cast<T>( 0 ); }
181 : for(unsigned i=0; i<B.cl; ++i) for(unsigned k=0; k<B.rw; ++k) C[i]+=A[k]*B(k,i);
182 : }
183 :
184 400 : template <typename T> void transpose( const Matrix<T>& A, Matrix<T>& AT ) {
185 400 : if( A.rw!=AT.cl || A.cl!=AT.rw ) AT.resize( A.cl, A.rw );
186 400 : for(unsigned i=0; i<A.cl; ++i) for(unsigned j=0; j<A.rw; ++j) AT(i,j)=A(j,i);
187 400 : }
188 :
189 : template <typename T> Log& operator<<(Log& ostr, const Matrix<T>& mat) {
190 : for(unsigned i=0; i<mat.sz; ++i) ostr<<mat.data[i]<<" ";
191 : return ostr;
192 : }
193 :
194 : template <typename T> void matrixOut( Log& ostr, const Matrix<T>& mat) {
195 : for(unsigned i=0; i<mat.rw; ++i) {
196 : for(unsigned j=0; j<mat.cl; ++j) { ostr<<mat(i,j)<<" "; }
197 : ostr<<"\n";
198 : }
199 : return;
200 : }
201 :
202 566815 : template <typename T> int diagMat( const Matrix<T>& A, std::vector<double>& eigenvals, Matrix<double>& eigenvecs ) {
203 :
204 : // Check matrix is square and symmetric
205 566815 : plumed_assert( A.rw==A.cl ); plumed_assert( A.isSymmetric()==1 );
206 566816 : double *da=new double[A.sz]; unsigned k=0; double *evals=new double[ A.cl ];
207 : // Transfer the matrix to the local array
208 566905 : for (unsigned i=0; i<A.cl; ++i) for (unsigned j=0; j<A.rw; ++j) da[k++]=static_cast<double>( A(j,i) );
209 :
210 566989 : int n=A.cl; int lwork=-1, liwork=-1, m, info, one=1;
211 566989 : double *work=new double[A.cl]; int *iwork=new int[A.cl];
212 566982 : double vl, vu, abstol=0.0;
213 566982 : int* isup=new int[2*A.cl]; double *evecs=new double[A.sz];
214 :
215 567056 : plumed_lapack_dsyevr("V", "I", "U", &n, da, &n, &vl, &vu, &one, &n,
216 : &abstol, &m, evals, evecs, &n,
217 : isup, work, &lwork, iwork, &liwork, &info);
218 566672 : if (info!=0) return info;
219 :
220 : // Retrieve correct sizes for work and iwork then reallocate
221 566653 : liwork=iwork[0]; delete [] iwork; iwork=new int[liwork];
222 567064 : lwork=static_cast<int>( work[0] ); delete [] work; work=new double[lwork];
223 :
224 567061 : plumed_lapack_dsyevr("V", "I", "U", &n, da, &n, &vl, &vu, &one, &n,
225 : &abstol, &m, evals, evecs, &n,
226 : isup, work, &lwork, iwork, &liwork, &info);
227 566864 : if (info!=0) return info;
228 :
229 566864 : if( eigenvals.size()!=A.cl ) { eigenvals.resize( A.cl ); }
230 567072 : if( eigenvecs.rw!=A.rw || eigenvecs.cl!=A.cl ) { eigenvecs.resize( A.rw, A.cl ); }
231 567279 : k=0;
232 2806147 : for(unsigned i=0; i<A.cl; ++i) {
233 2239094 : eigenvals[i]=evals[i];
234 : // N.B. For ease of producing projectors we store the eigenvectors
235 : // ROW-WISE in the eigenvectors matrix. The first index is the
236 : // eigenvector number and the second the component
237 2240104 : for(unsigned j=0; j<A.rw; ++j) { eigenvecs(i,j)=evecs[k++]; }
238 : }
239 :
240 : // This changes eigenvectors so that the first non-null element
241 : // of each of them is positive
242 : // We can do it because the phase is arbitrary, and helps making
243 : // the result reproducible
244 2805702 : for(int i=0; i<n; ++i) {
245 : int j;
246 2238758 : for(j=0; j<n; j++) if(eigenvecs(i,j)*eigenvecs(i,j)>1e-14) break;
247 2238876 : if(j<n) if(eigenvecs(i,j)<0.0) for(j=0; j<n; j++) eigenvecs(i,j)*=-1;
248 : }
249 :
250 : // Deallocate all the memory used by the various arrays
251 566944 : delete[] da; delete [] work; delete [] evals; delete[] evecs; delete [] iwork; delete [] isup;
252 566938 : return 0;
253 : }
254 :
255 : template <typename T> int pseudoInvert( const Matrix<T>& A, Matrix<double>& pseudoinverse ) {
256 : double *da=new double[A.sz]; unsigned k=0;
257 : // Transfer the matrix to the local array
258 : for (unsigned i=0; i<A.cl; ++i) for (unsigned j=0; j<A.rw; ++j) da[k++]=static_cast<double>( A(j,i) );
259 :
260 : int nsv, info, nrows=A.rw, ncols=A.cl;
261 : if(A.rw>A.cl) {nsv=A.cl;} else {nsv=A.rw;}
262 :
263 : // Create some containers for stuff from single value decomposition
264 : double *S=new double[nsv]; double *U=new double[nrows*nrows];
265 : double *VT=new double[ncols*ncols]; int *iwork=new int[8*nsv];
266 :
267 : // This optimizes the size of the work array used in lapack singular value decomposition
268 : int lwork=-1; double* work=new double[1];
269 : plumed_lapack_dgesdd( "A", &nrows, &ncols, da, &nrows, S, U, &nrows, VT, &ncols, work, &lwork, iwork, &info );
270 : if(info!=0) return info;
271 :
272 : // Retrieve correct sizes for work and rellocate
273 : lwork=(int) work[0]; delete [] work; work=new double[lwork];
274 :
275 : // This does the singular value decomposition
276 : plumed_lapack_dgesdd( "A", &nrows, &ncols, da, &nrows, S, U, &nrows, VT, &ncols, work, &lwork, iwork, &info );
277 : if(info!=0) return info;
278 :
279 : // Compute the tolerance on the singular values ( machine epsilon * number of singular values * maximum singular value )
280 : double tol; tol=S[0]; for(int i=1; i<nsv; ++i) { if( S[i]>tol ) { tol=S[i]; } } tol*=nsv*epsilon;
281 :
282 : // Get the inverses of the singlular values
283 : Matrix<double> Si( ncols, nrows ); Si=0.0;
284 : for(int i=0; i<nsv; ++i) { if( S[i]>tol ) { Si(i,i)=1./S[i]; } else { Si(i,i)=0.0; } }
285 :
286 : // Now extract matrices for pseudoinverse
287 : Matrix<double> V( ncols, ncols ), UT( nrows, nrows ), tmp( ncols, nrows );
288 : k=0; for(int i=0; i<nrows; ++i) { for(int j=0; j<nrows; ++j) { UT(i,j)=U[k++]; } }
289 : k=0; for(int i=0; i<ncols; ++i) { for(int j=0; j<ncols; ++j) { V(i,j)=VT[k++]; } }
290 :
291 : // And do matrix algebra to construct the pseudoinverse
292 : if( pseudoinverse.rw!=ncols || pseudoinverse.cl!=nrows ) pseudoinverse.resize( ncols, nrows );
293 : mult( V, Si, tmp ); mult( tmp, UT, pseudoinverse );
294 :
295 : // Deallocate all the memory
296 : delete [] S; delete [] U; delete [] VT; delete [] work; delete [] iwork; delete [] da;
297 : return 0;
298 : }
299 :
300 9238 : template <typename T> int Invert( const Matrix<T>& A, Matrix<double>& inverse ) {
301 :
302 9238 : if( A.isSymmetric()==1 ) {
303 : // GAT -- I only ever use symmetric matrices so I can invert them like this.
304 : // I choose to do this as I have had problems with the more general way of doing this that
305 : // is implemented below.
306 18364 : std::vector<double> eval(A.rw); Matrix<double> evec(A.rw,A.cl), tevec(A.rw,A.cl);
307 9182 : int err; err=diagMat( A, eval, evec );
308 9182 : if(err!=0) return err;
309 9182 : for (unsigned i=0; i<A.rw; ++i) for (unsigned j=0; j<A.cl; ++j) tevec(i,j)=evec(j,i)/eval[j];
310 18364 : mult(tevec,evec,inverse);
311 : } else {
312 54 : double *da=new double[A.sz]; int *ipiv=new int[A.cl];
313 53 : unsigned k=0; int n=A.rw, info;
314 53 : for(unsigned i=0; i<A.cl; ++i) for(unsigned j=0; j<A.rw; ++j) da[k++]=static_cast<double>( A(j,i) );
315 :
316 57 : plumed_lapack_dgetrf(&n,&n,da,&n,ipiv,&info);
317 55 : if(info!=0) return info;
318 :
319 55 : int lwork=-1; double* work=new double[A.cl];
320 56 : plumed_lapack_dgetri(&n,da,&n,ipiv,work,&lwork,&info);
321 54 : if(info!=0) return info;
322 :
323 54 : lwork=static_cast<int>( work[0] ); delete [] work; work=new double[lwork];
324 56 : plumed_lapack_dgetri(&n,da,&n,ipiv,work,&lwork,&info);
325 54 : if(info!=0) return info;
326 :
327 54 : if( inverse.cl!=A.cl || inverse.rw!=A.rw ) { inverse.resize(A.rw,A.cl); }
328 53 : k=0; for(unsigned i=0; i<A.rw; ++i) for(unsigned j=0; j<A.cl; ++j) inverse(j,i)=da[k++];
329 :
330 57 : delete [] da; delete[] work; delete[] ipiv;
331 : }
332 :
333 9239 : return 0;
334 : }
335 :
336 446 : template <typename T> void cholesky( const Matrix<T>& A, Matrix<T>& B ) {
337 :
338 446 : plumed_assert( A.rw==A.cl && A.isSymmetric() );
339 446 : Matrix<T> L(A.rw,A.cl); L=0.;
340 892 : std::vector<T> D(A.rw,0.);
341 1047 : for(unsigned i=0; i<A.rw; ++i) {
342 601 : L(i,i)=static_cast<T>( 1 );
343 756 : for (unsigned j=0; j<i; ++j) {
344 155 : L(i,j)=A(i,j);
345 155 : for (unsigned k=0; k<j; ++k) L(i,j)-=L(i,k)*L(j,k)*D[k];
346 155 : if (D[j]!=0.) L(i,j)/=D[j]; else L(i,j)=static_cast<T>( 0 );
347 : }
348 601 : D[i]=A(i,i);
349 601 : for (unsigned k=0; k<i; ++k) D[i]-=L(i,k)*L(i,k)*D[k];
350 : }
351 :
352 446 : for(unsigned i=0; i<A.rw; ++i) D[i]=(D[i]>0.?sqrt(D[i]):0.);
353 446 : if( B.rw!=A.rw || B.cl!=A.cl ) { B.resize( A.rw, A.cl); }
354 892 : 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];
355 446 : }
356 :
357 : template <typename T> void chol_elsolve( const Matrix<T>& M, const std::vector<T>& b, std::vector<T>& y ) {
358 :
359 : plumed_assert( M.rw==M.cl && M(0,1)==0.0 && b.size()==M.rw );
360 : if( y.size()!=M.rw ) { y.resize( M.rw ); }
361 : for(unsigned i=0; i<M.rw; ++i) {
362 : y[i]=b[i];
363 : for(unsigned j=0; j<i; ++j) y[i]-=M(i,j)*y[j];
364 : y[i]*=1.0/M(i,i);
365 : }
366 : }
367 :
368 0 : template <typename T> int logdet( const Matrix<T>& M, double& ldet ) {
369 : // Check matrix is square and symmetric
370 0 : plumed_assert( M.rw==M.cl || M.isSymmetric() );
371 :
372 0 : double *da=new double[M.sz]; unsigned k=0; double *evals=new double[M.cl];
373 : // Transfer the matrix to the local array
374 0 : for (unsigned i=0; i<M.rw; ++i) for (unsigned j=0; j<M.cl; ++j) da[k++]=static_cast<double>( M(j,i) );
375 :
376 0 : int n=M.cl; int lwork=-1, liwork=-1, info, m, one=1;
377 0 : double *work=new double[M.rw]; int *iwork=new int[M.rw];
378 0 : double vl, vu, abstol=0.0;
379 0 : int* isup=new int[2*M.rw]; double *evecs=new double[M.sz];
380 0 : plumed_lapack_dsyevr("N", "I", "U", &n, da, &n, &vl, &vu, &one, &n,
381 : &abstol, &m, evals, evecs, &n,
382 : isup, work, &lwork, iwork, &liwork, &info);
383 0 : if (info!=0) return info;
384 :
385 : // Retrieve correct sizes for work and iwork then reallocate
386 0 : lwork=static_cast<int>( work[0] ); delete [] work; work=new double[lwork];
387 0 : liwork=iwork[0]; delete [] iwork; iwork=new int[liwork];
388 :
389 0 : plumed_lapack_dsyevr("N", "I", "U", &n, da, &n, &vl, &vu, &one, &n,
390 : &abstol, &m, evals, evecs, &n,
391 : isup, work, &lwork, iwork, &liwork, &info);
392 0 : if (info!=0) return info;
393 :
394 : // Transfer the eigenvalues and eigenvectors to the output
395 0 : ldet=0; for(unsigned i=0; i<M.cl; i++) { ldet+=log(evals[i]); }
396 :
397 : // Deallocate all the memory used by the various arrays
398 0 : delete[] da; delete [] work; delete [] evals; delete[] evecs; delete [] iwork; delete [] isup;
399 :
400 0 : return 0;
401 : }
402 :
403 :
404 :
405 : }
406 : #endif
|
