Skip to content

Commit

Permalink
Fixed incorrect array size in diagMatSym
Browse files Browse the repository at this point in the history
According to dsyevr documentation:
- the matrix being diagonalized is overwritten.
  this is not happning in practice, but I added a temporary matrix
- the vector of eigenvalues should hold N elements even if we
  only compute M<N elements.

The latter point is important when we diagonalize a 4x4 matrix and
only need the lowest eigenvalue (as in RMSD calculations).
Notice that in practice the next element of the array
is only accessed if there are two numerically identical eigenvalues,
which happens when atoms are aligned on a straight line
(unlikely)
  • Loading branch information
GiovanniBussi committed Mar 4, 2024
1 parent e71a88a commit f1da0a9
Showing 1 changed file with 7 additions and 2 deletions.
9 changes: 7 additions & 2 deletions src/tools/Tensor.h
Expand Up @@ -535,6 +535,10 @@ void diagMatSym(const TensorGeneric<n,n>&mat,VectorGeneric<m>&evals,TensorGeneri
std::array<int,10*n> iwork;
std::array<double,(6+bs)*n> work;
std::array<int,2*m> isup;
// documentation says that "matrix is destroyed" !!!
auto mat_copy=mat;
// documentation says this is size n (even if m<n)
std::array<double,n> evals_tmp;
int nn=n; // dimension of matrix
double vl=0.0, vu=1.0; // ranges - not used
int one=1,mm=m; // minimum and maximum index
Expand All @@ -543,13 +547,14 @@ void diagMatSym(const TensorGeneric<n,n>&mat,VectorGeneric<m>&evals,TensorGeneri
int info=0; // result
int liwork=iwork.size();
int lwork=work.size();
TensorGenericAux::local_dsyevr("V", (n==m?"A":"I"), "U", &nn, const_cast<double*>(&mat[0][0]), &nn, &vl, &vu, &one, &mm,
&abstol, &mout, &evals[0], &evec[0][0], &nn,
TensorGenericAux::local_dsyevr("V", (n==m?"A":"I"), "U", &nn, const_cast<double*>(&mat_copy[0][0]), &nn, &vl, &vu, &one, &mm,
&abstol, &mout, &evals_tmp[0], &evec[0][0], &nn,
isup.data(), work.data(), &lwork, iwork.data(), &liwork, &info);
if(info!=0) plumed_error()<<"Error diagonalizing matrix\n"
<<"Matrix:\n"<<mat<<"\n"
<<"Info: "<<info<<"\n";
plumed_assert(mout==m);
for(unsigned i=0; i<m; i++) evals[i]=evals_tmp[i];
// This changes eigenvectors so that the first non-null element
// of each of them is positive
// We can do it because the phase is arbitrary, and helps making
Expand Down

0 comments on commit f1da0a9

Please sign in to comment.