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 : template <typename T> T norm( const std::vector<T>& A ) {
48 : T val{};
49 : for(unsigned i=0; i<A.size(); ++i) {
50 : val+=A[i]*A[i];
51 : }
52 : 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 9736587 : 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<double>&, Matrix<double>& );
75 : /// Calculate the Moore-Penrose Pseudoinverse of a matrix
76 : template <typename U> friend int pseudoInvert( const Matrix<U>&, Matrix<double>& );
77 : /// Calculate the logarithm of the determinant of a symmetric matrix - returns zero if succesfull
78 : template <typename U> friend int logdet( const Matrix<U>&, double& );
79 : /// Invert a matrix (works for both symmetric and asymmetric matrices) - returns zero if sucesfull
80 : template <typename U> friend int Invert( const Matrix<U>&, Matrix<double>& );
81 : /// Do a cholesky decomposition of a matrix
82 : template <typename U> friend void cholesky( const Matrix<U>&, Matrix<U>& );
83 : /// Solve a system of equations using the cholesky decomposition
84 : template <typename U> friend void chol_elsolve( const Matrix<U>&, const std::vector<U>&, std::vector<U>& );
85 : private:
86 : /// Number of elements in matrix (nrows*ncols)
87 : unsigned sz;
88 : /// Number of rows in matrix
89 : unsigned rw;
90 : /// Number of columns in matrix
91 : unsigned cl;
92 : /// The data in the matrix
93 : std::vector<T> data;
94 : public:
95 10771705 : explicit Matrix(const unsigned nr=0, const unsigned nc=0 ) : sz(nr*nc), rw(nr), cl(nc), data(nr*nc) {}
96 259 : Matrix(const Matrix<T>& t) : sz(t.sz), rw(t.rw), cl(t.cl), data(t.data) {}
97 : /// Resize the matrix
98 : void resize( const unsigned nr, const unsigned nc ) {
99 390 : rw=nr;
100 392 : cl=nc;
101 392 : sz=nr*nc;
102 392 : data.resize(sz);
103 177 : }
104 : /// Return the number of rows
105 : inline unsigned nrows() const {
106 37740101 : return rw;
107 : }
108 : /// Return the number of columns
109 : inline unsigned ncols() const {
110 116552290 : return cl;
111 : }
112 : /// Return the contents of the matrix as a vector of length rw*cl
113 : inline std::vector<T>& getVector() {
114 26 : return data;
115 : }
116 : /// Set the matrix from a vector input
117 4 : inline void setFromVector( const std::vector<T>& vecin ) {
118 4 : plumed_assert( vecin.size()==sz );
119 904 : for(unsigned i=0; i<sz; ++i) {
120 900 : data[i]=vecin[i];
121 : }
122 4 : }
123 : /// Return element i,j of the matrix
124 : inline const T& operator () (const unsigned& i, const unsigned& j) const {
125 2267477073 : return data[j+i*cl];
126 : }
127 : /// Return a referenre to element i,j of the matrix
128 : inline T& operator () (const unsigned& i, const unsigned& j) {
129 230283586 : return data[j+i*cl];
130 : }
131 : /// Set all elements of the matrix equal to the value of v
132 : Matrix<T>& operator=(const T& v) {
133 6748667 : for(unsigned i=0; i<sz; ++i) {
134 6733743 : data[i]=v;
135 : }
136 : return *this;
137 : }
138 : /// Set the Matrix equal to another Matrix
139 : Matrix<T>& operator=(const Matrix<T>& m) {
140 5208 : sz=m.sz;
141 5208 : rw=m.rw;
142 5208 : cl=m.cl;
143 5208 : data=m.data;
144 5208 : return *this;
145 : }
146 : /// Set the Matrix equal to the value of a standard vector - used for readin
147 : Matrix<T>& operator=(const std::vector<T>& v) {
148 : plumed_dbg_assert( v.size()==sz );
149 : for(unsigned i=0; i<sz; ++i) {
150 : data[i]=v[i];
151 : }
152 : return *this;
153 : }
154 : /// Add v to all elements of the Matrix
155 : Matrix<T> operator+=(const T& v) {
156 : for(unsigned i=0; i<sz; ++i) {
157 : data[i]+=v;
158 : }
159 : return *this;
160 : }
161 : /// Multiply all elements by v
162 112 : Matrix<T> operator*=(const T& v) {
163 440804 : for(unsigned i=0; i<sz; ++i) {
164 440692 : data[i]*=v;
165 : }
166 112 : return *this;
167 : }
168 : /// Matrix addition
169 : Matrix<T>& operator+=(const Matrix<T>& m) {
170 : plumed_dbg_assert( m.rw==rw && m.cl==cl );
171 : data+=m.data;
172 : return *this;
173 : }
174 : /// Subtract v from all elements of the Matrix
175 : Matrix<T> operator-=(const T& v) {
176 : for(unsigned i=0; i<sz; ++i) {
177 : data-=v;
178 : }
179 : return *this;
180 : }
181 : /// Matrix subtraction
182 : Matrix<T>& operator-=(const Matrix<T>& m) {
183 : plumed_dbg_assert( m.rw==rw && m.cl==cl );
184 : data-=m.data;
185 : return *this;
186 : }
187 : /// Test if the matrix is symmetric or not
188 23980 : unsigned isSymmetric() const {
189 23980 : if (rw!=cl) {
190 : return 0;
191 : }
192 : unsigned sym=1;
193 48894 : for(unsigned i=1; i<rw; ++i)
194 244923 : for(unsigned j=0; j<i; ++j)
195 220066 : if( std::fabs(data[i+j*cl]-data[j+i*cl])>1.e-10 ) {
196 : sym=0;
197 : break;
198 : }
199 : return sym;
200 : }
201 : };
202 :
203 : /// Multiply matrix by scalar
204 12 : template <typename T> Matrix<T> operator*(T& v, const Matrix<T>& m ) {
205 : Matrix<T> new_m(m);
206 12 : new_m*=v;
207 12 : return new_m;
208 : }
209 :
210 13725 : template <typename T> void mult( const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& C ) {
211 13725 : plumed_assert(A.cl==B.rw);
212 13725 : if( A.rw !=C.rw || B.cl !=C.cl ) {
213 0 : C.resize( A.rw, B.cl );
214 : }
215 : C=static_cast<T>( 0 );
216 55954 : for(unsigned i=0; i<A.rw; ++i)
217 4110794 : for(unsigned j=0; j<B.cl; ++j)
218 2011178606 : for (unsigned k=0; k<A.cl; ++k) {
219 2007110041 : C(i,j)+=A(i,k)*B(k,j);
220 : }
221 13725 : }
222 :
223 2141 : template <typename T> void mult( const Matrix<T>& A, const std::vector<T>& B, std::vector<T>& C) {
224 2141 : plumed_assert( A.cl==B.size() );
225 2141 : if( C.size()!=A.rw ) {
226 94 : C.resize(A.rw);
227 : }
228 37323 : for(unsigned i=0; i<A.rw; ++i) {
229 35182 : C[i]= static_cast<T>( 0 );
230 : }
231 37323 : for(unsigned i=0; i<A.rw; ++i)
232 679318 : for(unsigned k=0; k<A.cl; ++k) {
233 644136 : C[i]+=A(i,k)*B[k] ;
234 : }
235 2141 : }
236 :
237 : template <typename T> void mult( const std::vector<T>& A, const Matrix<T>& B, std::vector<T>& C) {
238 : plumed_assert( B.rw==A.size() );
239 : if( C.size()!=B.cl ) {
240 : C.resize( B.cl );
241 : }
242 : for(unsigned i=0; i<B.cl; ++i) {
243 : C[i]=static_cast<T>( 0 );
244 : }
245 : for(unsigned i=0; i<B.cl; ++i)
246 : for(unsigned k=0; k<B.rw; ++k) {
247 : C[i]+=A[k]*B(k,i);
248 : }
249 : }
250 :
251 : template <typename T> void transpose( const Matrix<T>& A, Matrix<T>& AT ) {
252 : if( A.rw!=AT.cl || A.cl!=AT.rw ) {
253 : AT.resize( A.cl, A.rw );
254 : }
255 : for(unsigned i=0; i<A.cl; ++i)
256 : for(unsigned j=0; j<A.rw; ++j) {
257 : AT(i,j)=A(j,i);
258 : }
259 : }
260 :
261 : template <typename T> Log& operator<<(Log& ostr, const Matrix<T>& mat) {
262 : for(unsigned i=0; i<mat.sz; ++i) {
263 : ostr<<mat.data[i]<<" ";
264 : }
265 : return ostr;
266 : }
267 :
268 : template <typename T> void matrixOut( Log& ostr, const Matrix<T>& mat) {
269 : for(unsigned i=0; i<mat.rw; ++i) {
270 : for(unsigned j=0; j<mat.cl; ++j) {
271 : ostr<<mat(i,j)<<" ";
272 : }
273 : ostr<<"\n";
274 : }
275 : return;
276 : }
277 :
278 14151 : template <typename T> int diagMat( const Matrix<T>& A, std::vector<double>& eigenvals, Matrix<double>& eigenvecs ) {
279 :
280 : // Check matrix is square and symmetric
281 14151 : plumed_assert( A.rw==A.cl );
282 14151 : plumed_assert( A.isSymmetric()==1 );
283 14151 : std::vector<double> da(A.sz);
284 : unsigned k=0;
285 14151 : std::vector<double> evals(A.cl);
286 : // Transfer the matrix to the local array
287 43840 : for (unsigned i=0; i<A.cl; ++i)
288 480492 : for (unsigned j=0; j<A.rw; ++j) {
289 450803 : da[k++]=static_cast<double>( A(j,i) );
290 : }
291 :
292 14151 : int n=A.cl;
293 14151 : int lwork=-1, liwork=-1, m, info, one=1;
294 14151 : std::vector<double> work(A.cl);
295 14151 : std::vector<int> iwork(A.cl);
296 14151 : double vl, vu, abstol=0.0;
297 14151 : std::vector<int> isup(2*A.cl);
298 14151 : std::vector<double> evecs(A.sz);
299 :
300 14151 : plumed_lapack_dsyevr("V", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
301 : &abstol, &m, evals.data(), evecs.data(), &n,
302 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
303 14151 : if (info!=0) {
304 : return info;
305 : }
306 :
307 : // Retrieve correct sizes for work and iwork then reallocate
308 14151 : liwork=iwork[0];
309 14151 : iwork.resize(liwork);
310 14151 : lwork=static_cast<int>( work[0] );
311 14151 : work.resize(lwork);
312 :
313 14151 : plumed_lapack_dsyevr("V", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
314 : &abstol, &m, evals.data(), evecs.data(), &n,
315 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
316 14151 : if (info!=0) {
317 : return info;
318 : }
319 :
320 14151 : if( eigenvals.size()!=A.cl ) {
321 0 : eigenvals.resize( A.cl );
322 : }
323 14151 : if( eigenvecs.rw!=A.rw || eigenvecs.cl!=A.cl ) {
324 0 : eigenvecs.resize( A.rw, A.cl );
325 : }
326 : k=0;
327 43840 : for(unsigned i=0; i<A.cl; ++i) {
328 29689 : eigenvals[i]=evals[i];
329 : // N.B. For ease of producing projectors we store the eigenvectors
330 : // ROW-WISE in the eigenvectors matrix. The first index is the
331 : // eigenvector number and the second the component
332 480492 : for(unsigned j=0; j<A.rw; ++j) {
333 450803 : eigenvecs(i,j)=evecs[k++];
334 : }
335 : }
336 :
337 : // This changes eigenvectors so that the first non-null element
338 : // of each of them is positive
339 : // We can do it because the phase is arbitrary, and helps making
340 : // the result reproducible
341 43840 : for(int i=0; i<n; ++i) {
342 : int j;
343 29689 : for(j=0; j<n; j++)
344 29689 : if(eigenvecs(i,j)*eigenvecs(i,j)>1e-14) {
345 : break;
346 : }
347 29689 : if(j<n)
348 29689 : if(eigenvecs(i,j)<0.0)
349 208257 : for(j=0; j<n; j++) {
350 203093 : eigenvecs(i,j)*=-1;
351 : }
352 : }
353 : return 0;
354 : }
355 :
356 1 : template <typename T> int pseudoInvert( const Matrix<T>& A, Matrix<double>& pseudoinverse ) {
357 1 : std::vector<double> da(A.sz);
358 : unsigned k=0;
359 : // Transfer the matrix to the local array
360 501 : for (unsigned i=0; i<A.cl; ++i)
361 250500 : for (unsigned j=0; j<A.rw; ++j) {
362 250000 : da[k++]=static_cast<double>( A(j,i) );
363 : }
364 :
365 1 : int nsv, info, nrows=A.rw, ncols=A.cl;
366 1 : if(A.rw>A.cl) {
367 : nsv=A.cl;
368 : } else {
369 : nsv=A.rw;
370 : }
371 :
372 : // Create some containers for stuff from single value decomposition
373 1 : std::vector<double> S(nsv);
374 1 : std::vector<double> U(nrows*nrows);
375 1 : std::vector<double> VT(ncols*ncols);
376 1 : std::vector<int> iwork(8*nsv);
377 :
378 : // This optimizes the size of the work array used in lapack singular value decomposition
379 1 : int lwork=-1;
380 1 : std::vector<double> work(1);
381 1 : plumed_lapack_dgesdd( "A", &nrows, &ncols, da.data(), &nrows, S.data(), U.data(), &nrows, VT.data(), &ncols, work.data(), &lwork, iwork.data(), &info );
382 1 : if(info!=0) {
383 : return info;
384 : }
385 :
386 : // Retrieve correct sizes for work and rellocate
387 1 : lwork=(int) work[0];
388 1 : work.resize(lwork);
389 :
390 : // This does the singular value decomposition
391 1 : plumed_lapack_dgesdd( "A", &nrows, &ncols, da.data(), &nrows, S.data(), U.data(), &nrows, VT.data(), &ncols, work.data(), &lwork, iwork.data(), &info );
392 1 : if(info!=0) {
393 : return info;
394 : }
395 :
396 : // Compute the tolerance on the singular values ( machine epsilon * number of singular values * maximum singular value )
397 : double tol;
398 1 : tol=S[0];
399 500 : for(int i=1; i<nsv; ++i) {
400 499 : if( S[i]>tol ) {
401 : tol=S[i];
402 : }
403 : }
404 1 : tol*=nsv*epsilon;
405 :
406 : // Get the inverses of the singlular values
407 1 : Matrix<double> Si( ncols, nrows );
408 : Si=0.0;
409 501 : for(int i=0; i<nsv; ++i) {
410 500 : if( S[i]>tol ) {
411 499 : Si(i,i)=1./S[i];
412 : } else {
413 1 : Si(i,i)=0.0;
414 : }
415 : }
416 :
417 : // Now extract matrices for pseudoinverse
418 1 : Matrix<double> V( ncols, ncols ), UT( nrows, nrows ), tmp( ncols, nrows );
419 : k=0;
420 501 : for(int i=0; i<nrows; ++i) {
421 250500 : for(int j=0; j<nrows; ++j) {
422 250000 : UT(i,j)=U[k++];
423 : }
424 : }
425 : k=0;
426 501 : for(int i=0; i<ncols; ++i) {
427 250500 : for(int j=0; j<ncols; ++j) {
428 250000 : V(i,j)=VT[k++];
429 : }
430 : }
431 :
432 : // And do matrix algebra to construct the pseudoinverse
433 1 : if( pseudoinverse.rw!=ncols || pseudoinverse.cl!=nrows ) {
434 0 : pseudoinverse.resize( ncols, nrows );
435 : }
436 1 : mult( V, Si, tmp );
437 1 : mult( tmp, UT, pseudoinverse );
438 :
439 : return 0;
440 : }
441 :
442 9383 : template <typename T> int Invert( const Matrix<T>& A, Matrix<double>& inverse ) {
443 :
444 9383 : if( A.isSymmetric()==1 ) {
445 : // GAT -- I only ever use symmetric matrices so I can invert them like this.
446 : // I choose to do this as I have had problems with the more general way of doing this that
447 : // is implemented below.
448 9326 : std::vector<double> eval(A.rw);
449 9326 : Matrix<double> evec(A.rw,A.cl), tevec(A.rw,A.cl);
450 : int err;
451 9326 : err=diagMat( A, eval, evec );
452 9326 : if(err!=0) {
453 : return err;
454 : }
455 27816 : for (unsigned i=0; i<A.rw; ++i)
456 55574 : for (unsigned j=0; j<A.cl; ++j) {
457 37084 : tevec(i,j)=evec(j,i)/eval[j];
458 : }
459 9326 : mult(tevec,evec,inverse);
460 : } else {
461 57 : std::vector<double> da(A.sz);
462 57 : std::vector<int> ipiv(A.cl);
463 : unsigned k=0;
464 57 : int n=A.rw, info;
465 171 : for(unsigned i=0; i<A.cl; ++i)
466 342 : for(unsigned j=0; j<A.rw; ++j) {
467 228 : da[k++]=static_cast<double>( A(j,i) );
468 : }
469 :
470 57 : plumed_lapack_dgetrf(&n,&n,da.data(),&n,ipiv.data(),&info);
471 57 : if(info!=0) {
472 : return info;
473 : }
474 :
475 57 : int lwork=-1;
476 57 : std::vector<double> work(A.cl);
477 57 : plumed_lapack_dgetri(&n,da.data(),&n,ipiv.data(),work.data(),&lwork,&info);
478 57 : if(info!=0) {
479 : return info;
480 : }
481 :
482 57 : lwork=static_cast<int>( work[0] );
483 57 : work.resize(lwork);
484 57 : plumed_lapack_dgetri(&n,da.data(),&n,ipiv.data(),work.data(),&lwork,&info);
485 57 : if(info!=0) {
486 : return info;
487 : }
488 :
489 57 : if( inverse.cl!=A.cl || inverse.rw!=A.rw ) {
490 0 : inverse.resize(A.rw,A.cl);
491 : }
492 : k=0;
493 171 : for(unsigned i=0; i<A.rw; ++i)
494 342 : for(unsigned j=0; j<A.cl; ++j) {
495 228 : inverse(j,i)=da[k++];
496 : }
497 : }
498 :
499 : return 0;
500 : }
501 :
502 446 : template <typename T> void cholesky( const Matrix<T>& A, Matrix<T>& B ) {
503 :
504 446 : plumed_assert( A.rw==A.cl && A.isSymmetric() );
505 : Matrix<T> L(A.rw,A.cl);
506 : L=0.;
507 446 : std::vector<T> D(A.rw,0.);
508 1047 : for(unsigned i=0; i<A.rw; ++i) {
509 601 : L(i,i)=static_cast<T>( 1 );
510 756 : for (unsigned j=0; j<i; ++j) {
511 155 : L(i,j)=A(i,j);
512 155 : for (unsigned k=0; k<j; ++k) {
513 0 : L(i,j)-=L(i,k)*L(j,k)*D[k];
514 : }
515 155 : if (D[j]!=0.) {
516 155 : L(i,j)/=D[j];
517 : } else {
518 0 : L(i,j)=static_cast<T>( 0 );
519 : }
520 : }
521 601 : D[i]=A(i,i);
522 756 : for (unsigned k=0; k<i; ++k) {
523 155 : D[i]-=L(i,k)*L(i,k)*D[k];
524 : }
525 : }
526 :
527 1047 : for(unsigned i=0; i<A.rw; ++i) {
528 601 : D[i]=(D[i]>0.?std::sqrt(D[i]):0.);
529 : }
530 446 : if( B.rw!=A.rw || B.cl!=A.cl ) {
531 0 : B.resize( A.rw, A.cl);
532 : }
533 : B=0.;
534 1047 : for(unsigned i=0; i<A.rw; ++i)
535 1357 : for(unsigned j=0; j<=i; ++j) {
536 756 : B(i,j)+=L(i,j)*D[j];
537 : }
538 446 : }
539 :
540 : template <typename T> void chol_elsolve( const Matrix<T>& M, const std::vector<T>& b, std::vector<T>& y ) {
541 :
542 : plumed_assert( M.rw==M.cl && M(0,1)==0.0 && b.size()==M.rw );
543 : if( y.size()!=M.rw ) {
544 : y.resize( M.rw );
545 : }
546 : for(unsigned i=0; i<M.rw; ++i) {
547 : y[i]=b[i];
548 : for(unsigned j=0; j<i; ++j) {
549 : y[i]-=M(i,j)*y[j];
550 : }
551 : y[i]*=1.0/M(i,i);
552 : }
553 : }
554 :
555 46 : template <typename T> int logdet( const Matrix<T>& M, double& ldet ) {
556 : // Check matrix is square and symmetric
557 46 : plumed_assert( M.rw==M.cl || M.isSymmetric() );
558 :
559 46 : std::vector<double> da(M.sz);
560 : unsigned k=0;
561 46 : std::vector<double> evals(M.cl);
562 : // Transfer the matrix to the local array
563 180 : for (unsigned i=0; i<M.rw; ++i)
564 532 : for (unsigned j=0; j<M.cl; ++j) {
565 398 : da[k++]=static_cast<double>( M(j,i) );
566 : }
567 :
568 46 : int n=M.cl;
569 46 : int lwork=-1, liwork=-1, info, m, one=1;
570 46 : std::vector<double> work(M.rw);
571 46 : std::vector<int> iwork(M.rw);
572 46 : double vl, vu, abstol=0.0;
573 46 : std::vector<int> isup(2*M.rw);
574 46 : std::vector<double> evecs(M.sz);
575 46 : plumed_lapack_dsyevr("N", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
576 : &abstol, &m, evals.data(), evecs.data(), &n,
577 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
578 46 : if (info!=0) {
579 : return info;
580 : }
581 :
582 : // Retrieve correct sizes for work and iwork then reallocate
583 46 : lwork=static_cast<int>( work[0] );
584 46 : work.resize(lwork);
585 46 : liwork=iwork[0];
586 46 : iwork.resize(liwork);
587 :
588 46 : plumed_lapack_dsyevr("N", "I", "U", &n, da.data(), &n, &vl, &vu, &one, &n,
589 : &abstol, &m, evals.data(), evecs.data(), &n,
590 : isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
591 46 : if (info!=0) {
592 : return info;
593 : }
594 :
595 : // Transfer the eigenvalues and eigenvectors to the output
596 46 : ldet=0;
597 180 : for(unsigned i=0; i<M.cl; i++) {
598 134 : ldet+=log(evals[i]);
599 : }
600 :
601 : return 0;
602 : }
603 :
604 :
605 :
606 : }
607 : #endif
|
