Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2017 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 : #include "MatrixOperationBase.h"
23 : #include "core/ActionRegister.h"
24 :
25 : #include <algorithm>
26 :
27 : //+PLUMEDOC ANALYSIS DIAGONALIZE
28 : /*
29 : Calculate the eigenvalues and eigenvectors of a square matrix
30 :
31 : This action allows you to use the [dsyevr](https://www.netlib.org/lapack/explore-html/d1/d56/group__heevr_gaa334ac0c11113576db0fc37b7565e8b5.html#gaa334ac0c11113576db0fc37b7565e8b5)
32 : function from the [LAPACK](https://www.netlib.org/lapack/explore-html/) library to calculate the [eigenvalues and eigenvectors](https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix)
33 : of a real symmetric matrix. For example, the following input can be used to calculate and print all four eigenvalues of the input [DISTANCE_MATRIX](DISTANCE_MATRIX.md).
34 :
35 : ```plumed
36 : d: DISTANCE_MATRIX GROUP=1-4
37 : diag: DIAGONALIZE ARG=d
38 : PRINT ARG=diag.vals-1,diag.vals-2,diag.vals-3,diag.vals-4 FILE=colvar
39 : ```
40 :
41 : If you wish to calculate only a subset of the eigenvalues and eigenvectors you would use an input like the one shown below. This input calculates and outputs the largest eigenvalue
42 : and its corresponding eigenvector.
43 :
44 : ```plumed
45 : d: DISTANCE_MATRIX GROUP=1-4
46 : diag: DIAGONALIZE ARG=d VECTORS=1
47 : PRINT ARG=diag.vals-1,diag.vecs-1 FILE=colvar
48 : ```
49 :
50 : */
51 : //+ENDPLUMEDOC
52 :
53 : namespace PLMD {
54 : namespace matrixtools {
55 :
56 : class DiagonalizeMatrix : public MatrixOperationBase {
57 : private:
58 : std::vector<unsigned> desired_vectors;
59 : Matrix<double> mymatrix;
60 : std::vector<double> eigvals;
61 : Matrix<double> eigvecs;
62 : /// For each requested eigenpair that has forces: the index of the eigenpair in eigvecs,
63 : /// the force on the eigenvalue and a vector used to compute the forces from the eigenvector
64 2436 : struct ForceTerm {
65 : unsigned ncol;
66 : double evalForce;
67 : bool hasEvecForce;
68 : std::vector<double> w;
69 : };
70 : std::vector<ForceTerm> forceTerms;
71 : public:
72 : static void registerKeywords( Keywords& keys );
73 : /// Constructor
74 : explicit DiagonalizeMatrix(const ActionOptions&);
75 : /// This is required to set the number of derivatives for the eigenvalues
76 328 : unsigned getNumberOfDerivatives() override {
77 328 : return getPntrToArgument(0)->getNumberOfValues();
78 : }
79 : ///
80 : void prepare() override ;
81 : ///
82 : void calculate() override ;
83 : /// Prepare the quantities needed to compute the forces on the matrix elements
84 : void apply() override ;
85 : ///
86 : double getForceOnMatrixElement( const unsigned& jrow, const unsigned& krow ) const override;
87 : };
88 :
89 : PLUMED_REGISTER_ACTION(DiagonalizeMatrix,"DIAGONALIZE")
90 :
91 58 : void DiagonalizeMatrix::registerKeywords( Keywords& keys ) {
92 58 : MatrixOperationBase::registerKeywords( keys );
93 58 : keys.add("compulsory","VECTORS","all","the eigenvalues and vectors that you would like to calculate. 1=largest, 2=second largest and so on");
94 116 : keys.addOutputComponent("vals","default","scalar","the eigevalues of the input matrix");
95 116 : keys.addOutputComponent("vecs","default","vector","the eigenvectors of the input matrix");
96 58 : }
97 :
98 30 : DiagonalizeMatrix::DiagonalizeMatrix(const ActionOptions& ao):
99 : Action(ao),
100 30 : MatrixOperationBase(ao) {
101 30 : if( getPntrToArgument(0)->getShape()[0]!=getPntrToArgument(0)->getShape()[1] ) {
102 0 : error("input matrix should be square");
103 : }
104 :
105 30 : std::vector<std::string> eigv;
106 60 : parseVector("VECTORS",eigv);
107 30 : if( eigv.size()>1 ) {
108 10 : Tools::interpretRanges(eigv);
109 10 : desired_vectors.resize( eigv.size() );
110 30 : for(unsigned i=0; i<eigv.size(); ++i) {
111 20 : Tools::convert( eigv[i], desired_vectors[i] );
112 : }
113 : } else {
114 20 : if( eigv.size()==0 ) {
115 0 : error("missing input to VECTORS keyword");
116 : }
117 : unsigned ivec;
118 20 : if( eigv[0]=="all" ) {
119 9 : desired_vectors.resize( getPntrToArgument(0)->getShape()[0] );
120 30 : for(unsigned i=0; i<desired_vectors.size(); ++i) {
121 21 : desired_vectors[i] = i + 1;
122 : }
123 : } else {
124 11 : Tools::convert( eigv[0], ivec );
125 11 : desired_vectors.resize(1);
126 11 : desired_vectors[0]=ivec;
127 : }
128 : }
129 :
130 : std::string num;
131 30 : std::vector<std::size_t> eval_shape(0);
132 30 : std::vector<std::size_t> evec_shape(1);
133 30 : evec_shape[0] = getPntrToArgument(0)->getShape()[0];
134 82 : for(unsigned i=0; i<desired_vectors.size(); ++i) {
135 52 : Tools::convert( desired_vectors[i], num );
136 52 : addComponent( "vals-" + num, eval_shape );
137 52 : componentIsNotPeriodic( "vals-" + num );
138 52 : addComponent( "vecs-" + num, evec_shape );
139 104 : componentIsNotPeriodic( "vecs-" + num );
140 : // Make sure eigenvalues are always stored
141 : }
142 :
143 30 : std::vector<unsigned> eigvecs_shape(2);
144 30 : eigvecs_shape[0]=eigvecs_shape[1]=getPntrToArgument(0)->getShape()[0];
145 30 : mymatrix.resize( eigvecs_shape[0], eigvecs_shape[1] );
146 30 : eigvals.resize( eigvecs_shape[0] );
147 30 : eigvecs.resize( eigvecs_shape[0], eigvecs_shape[1] );
148 60 : }
149 :
150 816 : void DiagonalizeMatrix::prepare() {
151 816 : std::vector<std::size_t> shape(1);
152 816 : shape[0]=getPntrToArgument(0)->getShape()[0];
153 2874 : for(unsigned i=0; i<desired_vectors.size(); ++i) {
154 2058 : if( getPntrToComponent( 2*i+1 )->getShape()[0]!=shape[0] ) {
155 9 : getPntrToComponent( 2*i+1 )->setShape( shape );
156 : }
157 : }
158 :
159 816 : }
160 :
161 815 : void DiagonalizeMatrix::calculate() {
162 815 : if( getPntrToArgument(0)->getShape()[0]==0 ) {
163 : return ;
164 : }
165 : // Resize stuff that might need resizing
166 815 : unsigned nvals=getPntrToArgument(0)->getShape()[0];
167 815 : if( eigvals.size()!=nvals ) {
168 5 : mymatrix.resize( nvals, nvals );
169 5 : eigvals.resize( nvals );
170 5 : eigvecs.resize( nvals, nvals );
171 : }
172 :
173 : // Retrieve the matrix from input
174 815 : retrieveFullMatrix( mymatrix );
175 : // Now diagonalize the matrix. The forces on an eigenvector depend on all the eigenvectors,
176 : // so they are all computed when derivatives are needed. Otherwise only the eigenpairs from
177 : // the smallest requested one to the largest eigenvalue are computed, which is faster
178 : // when only a few of the largest eigenvalues are requested (e.g. VECTORS=1)
179 815 : if( doNotCalculateDerivatives() ) {
180 19 : const unsigned maxvec=*std::max_element( desired_vectors.begin(), desired_vectors.end() );
181 19 : diagMatRange( mymatrix, nvals-maxvec, nvals-1, eigvals, eigvecs );
182 : } else {
183 796 : diagMat( mymatrix, eigvals, eigvecs );
184 : }
185 : // And set the eigenvalues and eigenvectors
186 2871 : for(unsigned i=0; i<desired_vectors.size(); ++i) {
187 2056 : getPntrToComponent(2*i)->set( eigvals[ mymatrix.ncols()-desired_vectors[i]] );
188 2056 : Value* evec_out = getPntrToComponent(2*i+1);
189 2056 : unsigned vreq = mymatrix.ncols()-desired_vectors[i];
190 12979 : for(unsigned j=0; j<mymatrix.ncols(); ++j) {
191 10923 : evec_out->set( j, eigvecs( vreq, j ) );
192 : }
193 : }
194 : }
195 :
196 799 : void DiagonalizeMatrix::apply() {
197 799 : if( doNotCalculateDerivatives() ) {
198 : return;
199 : }
200 : // The derivative of eigenvector c with respect to element (j,k) of the matrix is
201 : // d v_c[n] / d A(j,k) = sum_{m!=c} v_m[n] v_m[j] v_c[k] / (lambda_c-lambda_m)
202 : // so the force on (j,k) from the forces f on v_c is v_c[k]*w[j], with
203 : // w[j] = sum_{m!=c} v_m[j] (v_m.f) / (lambda_c-lambda_m)
204 : // computing w once here costs O(N^2), instead of O(N^2) for each matrix element
205 796 : const unsigned nvals=eigvals.size();
206 796 : forceTerms.clear();
207 2815 : for(unsigned i=0; i<desired_vectors.size(); ++i) {
208 2019 : const bool evalForce=getConstPntrToComponent(2*i)->forcesWereAdded();
209 2019 : const bool evecForce=getConstPntrToComponent(2*i+1)->forcesWereAdded();
210 2019 : if( !evalForce && !evecForce ) {
211 1208 : continue;
212 : }
213 : ForceTerm ft;
214 811 : ft.ncol=nvals-desired_vectors[i];
215 811 : ft.evalForce=( evalForce ? getConstPntrToComponent(2*i)->getForce(0) : 0.0 );
216 811 : ft.hasEvecForce=evecForce;
217 811 : if( evecForce ) {
218 136 : const Value* evec=getConstPntrToComponent(2*i+1);
219 136 : ft.w.assign(nvals,0.0);
220 1921 : for(unsigned m=0; m<nvals; ++m) {
221 1785 : if( m==ft.ncol ) {
222 136 : continue;
223 : }
224 : double vf=0;
225 24021 : for(unsigned n=0; n<nvals; ++n) {
226 22372 : vf += eigvecs(m,n)*evec->getForce(n);
227 : }
228 1649 : vf /= (eigvals[ft.ncol]-eigvals[m]);
229 24021 : for(unsigned j=0; j<nvals; ++j) {
230 22372 : ft.w[j] += vf*eigvecs(m,j);
231 : }
232 : }
233 : }
234 : forceTerms.push_back(std::move(ft));
235 : }
236 796 : MatrixOperationBase::apply();
237 : }
238 :
239 34161 : double DiagonalizeMatrix::getForceOnMatrixElement( const unsigned& jrow, const unsigned& kcol ) const {
240 : double ff = 0;
241 68562 : for(const auto & ft : forceTerms) {
242 : // Forces on eigenvalues
243 34401 : ff += ft.evalForce*eigvecs(ft.ncol,jrow)*eigvecs(ft.ncol,kcol);
244 : // Forces on eigenvectors
245 34401 : if( ft.hasEvecForce ) {
246 23086 : ff += ft.w[jrow]*eigvecs(ft.ncol,kcol);
247 : }
248 : }
249 34161 : return ff;
250 : }
251 :
252 : }
253 : }
|
