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 : #include "ActionWithMatrix.h"
23 : #include "tools/Communicator.h"
24 :
25 : namespace PLMD {
26 :
27 3240 : void ActionWithMatrix::registerKeywords( Keywords& keys ) {
28 3240 : ActionWithVector::registerKeywords( keys );
29 3240 : keys.use("ARG");
30 3240 : }
31 :
32 1474 : ActionWithMatrix::ActionWithMatrix(const ActionOptions&ao):
33 : Action(ao),
34 : ActionWithVector(ao),
35 1474 : next_action_in_chain(NULL),
36 1474 : matrix_to_do_before(NULL),
37 1474 : matrix_to_do_after(NULL),
38 1474 : clearOnEachCycle(true) {
39 1474 : }
40 :
41 1474 : ActionWithMatrix::~ActionWithMatrix() {
42 1474 : if( matrix_to_do_before ) {
43 828 : matrix_to_do_before->matrix_to_do_after=NULL;
44 828 : matrix_to_do_before->next_action_in_chain=NULL;
45 : }
46 1474 : }
47 :
48 16 : void ActionWithMatrix::getAllActionLabelsInMatrixChain( std::vector<std::string>& mylabels ) const {
49 : bool found=false;
50 24 : for(unsigned i=0; i<mylabels.size(); ++i) {
51 8 : if( getLabel()==mylabels[i] ) {
52 : found=true;
53 : }
54 : }
55 16 : if( !found ) {
56 16 : mylabels.push_back( getLabel() );
57 : }
58 16 : if( matrix_to_do_after ) {
59 8 : matrix_to_do_after->getAllActionLabelsInMatrixChain( mylabels );
60 : }
61 16 : }
62 :
63 55936 : void ActionWithMatrix::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
64 55936 : ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
65 :
66 121744 : for(int i=0; i<getNumberOfComponents(); ++i) {
67 65808 : Value* myval=getPntrToComponent(i);
68 65808 : if( myval->getRank()!=2 || myval->hasDerivatives() ) {
69 30419 : continue;
70 : }
71 35389 : myval->setPositionInMatrixStash(nmat);
72 35389 : nmat++;
73 35389 : if( !myval->valueIsStored() ) {
74 27595 : continue;
75 : }
76 7794 : if( myval->getShape()[1]>maxcol ) {
77 7133 : maxcol=myval->getShape()[1];
78 : }
79 : myval->setMatrixBookeepingStart(nbookeeping);
80 7794 : nbookeeping += myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
81 : }
82 : // Turn off clearning of derivatives after each matrix run if there are no matrices in the output of this action
83 55936 : clearOnEachCycle = false;
84 86355 : for(int i=0; i<getNumberOfComponents(); ++i) {
85 60994 : const Value* myval=getConstPntrToComponent(i);
86 60994 : if( myval->getRank()==2 && !myval->hasDerivatives() ) {
87 30575 : clearOnEachCycle = true;
88 30575 : break;
89 : }
90 : }
91 : // Turn off clearing of derivatives if we have only the values of adjacency matrices
92 55936 : if( doNotCalculateDerivatives() && isAdjacencyMatrix() ) {
93 852 : clearOnEachCycle = false;
94 : }
95 55936 : }
96 :
97 3317 : void ActionWithMatrix::finishChainBuild( ActionWithVector* act ) {
98 3317 : ActionWithMatrix* am=dynamic_cast<ActionWithMatrix*>(act);
99 3317 : if( !am || act==this ) {
100 : return;
101 : }
102 : // Build the list that contains everything we are going to loop over in getTotalMatrixBookeepgin and updateAllNeighbourLists
103 2269 : if( next_action_in_chain ) {
104 1382 : next_action_in_chain->finishChainBuild( act );
105 : } else {
106 887 : next_action_in_chain=am;
107 : // Build the list of things we are going to loop over in runTask
108 887 : if( am->isAdjacencyMatrix() || act->getName()=="VSTACK" ) {
109 59 : return ;
110 : }
111 828 : plumed_massert( !matrix_to_do_after, "cannot add " + act->getLabel() + " in " + getLabel() + " as have already added " + matrix_to_do_after->getLabel() );
112 828 : matrix_to_do_after=am;
113 828 : am->matrix_to_do_before=this;
114 : }
115 : }
116 :
117 1263 : const ActionWithMatrix* ActionWithMatrix::getFirstMatrixInChain() const {
118 1263 : if( !actionInChain() ) {
119 : return this;
120 : }
121 410 : return matrix_to_do_before->getFirstMatrixInChain();
122 : }
123 :
124 30970 : void ActionWithMatrix::getTotalMatrixBookeeping( unsigned& nbookeeping ) {
125 67921 : for(int i=0; i<getNumberOfComponents(); ++i) {
126 36951 : Value* myval=getPntrToComponent(i);
127 36951 : if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() ) {
128 32170 : continue;
129 : }
130 4781 : myval->reshapeMatrixStore( getNumberOfColumns() );
131 4781 : nbookeeping += myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
132 : }
133 30970 : if( next_action_in_chain ) {
134 13487 : next_action_in_chain->getTotalMatrixBookeeping( nbookeeping );
135 : }
136 30970 : }
137 :
138 31012 : void ActionWithMatrix::calculate() {
139 31012 : if( actionInChain() ) {
140 13529 : return ;
141 : }
142 : // Update all the neighbour lists
143 17483 : updateAllNeighbourLists();
144 : // Setup the matrix indices
145 17483 : unsigned nbookeeping=0;
146 17483 : getTotalMatrixBookeeping( nbookeeping );
147 17483 : if( matrix_bookeeping.size()!=nbookeeping ) {
148 376 : matrix_bookeeping.resize( nbookeeping );
149 : }
150 : std::fill( matrix_bookeeping.begin(), matrix_bookeeping.end(), 0 );
151 : // And run all the tasks
152 17483 : runAllTasks();
153 : }
154 :
155 30970 : void ActionWithMatrix::updateAllNeighbourLists() {
156 30970 : updateNeighbourList();
157 30970 : if( next_action_in_chain ) {
158 13487 : next_action_in_chain->updateAllNeighbourLists();
159 : }
160 30970 : }
161 :
162 1379848 : void ActionWithMatrix::performTask( const unsigned& task_index, MultiValue& myvals ) const {
163 : std::vector<unsigned> & indices( myvals.getIndices() );
164 1379848 : if( matrix_to_do_before ) {
165 : plumed_dbg_assert( myvals.inVectorCall() );
166 714734 : runEndOfRowJobs( task_index, indices, myvals );
167 714734 : return;
168 : }
169 665114 : setupForTask( task_index, indices, myvals );
170 :
171 : // Now loop over the row of the matrix
172 : unsigned ntwo_atoms = myvals.getSplitIndex();
173 29791135 : for(unsigned i=1; i<ntwo_atoms; ++i) {
174 : // This does everything in the stream that is done with single matrix elements
175 29126021 : runTask( getLabel(), task_index, indices[i], myvals );
176 : // Now clear only elements that are not accumulated over whole row
177 29126021 : clearMatrixElements( myvals );
178 : }
179 : // This updates the jobs that need to be completed when we get to the end of a row of the matrix
180 665114 : runEndOfRowJobs( task_index, indices, myvals );
181 : }
182 :
183 88131057 : void ActionWithMatrix::runTask( const std::string& controller, const unsigned& current, const unsigned colno, MultiValue& myvals ) const {
184 : double outval=0;
185 88131057 : myvals.setTaskIndex(current);
186 88131057 : myvals.setSecondTaskIndex( colno );
187 88131057 : if( isActive() ) {
188 87054595 : performTask( controller, current, colno, myvals );
189 : }
190 88131057 : bool hasval = !isAdjacencyMatrix();
191 172328384 : for(int i=0; i<getNumberOfComponents(); ++i) {
192 145389193 : if( fabs(myvals.get( getConstPntrToComponent(i)->getPositionInStream()) )>0 ) {
193 : hasval=true;
194 : break;
195 : }
196 : }
197 :
198 88131057 : if( hasval ) {
199 284623420 : for(int i=0; i<getNumberOfComponents(); ++i) {
200 201333323 : const Value* myval=getConstPntrToComponent(i);
201 201333323 : if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() ) {
202 187163511 : continue;
203 : }
204 14169812 : unsigned matindex = myval->getPositionInMatrixStash(), matbook_start = myval->getMatrixBookeepingStart(), col_stash_index = colno;
205 14169812 : if( colno>=myval->getShape()[0] ) {
206 7630130 : col_stash_index = colno - myval->getShape()[0];
207 : }
208 14169812 : unsigned rowstart = matbook_start+current*(1+myval->getNumberOfColumns());
209 14169812 : if( myval->forcesWereAdded() ) {
210 1551066 : unsigned sind = myval->getPositionInStream(), find = myvals.getMatrixBookeeping()[rowstart];
211 1551066 : double fforce = myval->getForce( myvals.getTaskIndex()*myval->getNumberOfColumns() + find );
212 1551066 : if( getNumberOfColumns()>=myval->getShape()[1] ) {
213 1414660 : fforce = myval->getForce( myvals.getTaskIndex()*myval->getShape()[1] + col_stash_index );
214 : }
215 20795223 : for(unsigned j=0; j<myvals.getNumberActive(sind); ++j) {
216 19244157 : unsigned kindex = myvals.getActiveIndex(sind,j);
217 19244157 : myvals.addMatrixForce( matindex, kindex, fforce*myvals.getDerivative(sind,kindex ) );
218 : }
219 : }
220 14169812 : double finalval = myvals.get( myval->getPositionInStream() );
221 14169812 : if( fabs(finalval)>0 ) {
222 8598870 : myvals.stashMatrixElement( matindex, rowstart, col_stash_index, finalval );
223 : }
224 : }
225 : }
226 88131057 : if( matrix_to_do_after ) {
227 59005036 : matrix_to_do_after->runTask( controller, current, colno, myvals );
228 : }
229 88131057 : }
230 :
231 19984 : void ActionWithMatrix::gatherThreads( const unsigned& nt, const unsigned& bufsize, const std::vector<double>& omp_buffer, std::vector<double>& buffer, MultiValue& myvals ) {
232 19984 : ActionWithVector::gatherThreads( nt, bufsize, omp_buffer, buffer, myvals );
233 64708664 : for(unsigned i=0; i<matrix_bookeeping.size(); ++i) {
234 64688680 : matrix_bookeeping[i] += myvals.getMatrixBookeeping()[i];
235 : }
236 19984 : }
237 :
238 17483 : void ActionWithMatrix::gatherProcesses( std::vector<double>& buffer ) {
239 17483 : ActionWithVector::gatherProcesses( buffer );
240 17483 : if( matrix_bookeeping.size()>0 && !runInSerial() ) {
241 4247 : comm.Sum( matrix_bookeeping );
242 : }
243 17483 : unsigned nval=0;
244 17483 : transferNonZeroMatrixElementsToValues( nval, matrix_bookeeping );
245 17483 : }
246 :
247 30970 : void ActionWithMatrix::transferNonZeroMatrixElementsToValues( unsigned& nval, const std::vector<unsigned>& matbook ) {
248 67921 : for(int i=0; i<getNumberOfComponents(); ++i) {
249 36951 : Value* myval=getPntrToComponent(i);
250 36951 : if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() || getNumberOfColumns()>=myval->getShape()[1] ) {
251 36852 : continue;
252 : }
253 99 : unsigned nelements = myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
254 5060331 : for(unsigned j=0; j<nelements; ++j) {
255 5060232 : myval->setMatrixBookeepingElement( j, matbook[nval+j] );
256 : }
257 99 : nval += nelements;
258 : }
259 30970 : if( next_action_in_chain ) {
260 13487 : next_action_in_chain->transferNonZeroMatrixElementsToValues( nval, matbook );
261 : }
262 30970 : }
263 :
264 389280 : void ActionWithMatrix::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
265 : const unsigned& bufstart, std::vector<double>& buffer ) const {
266 389280 : if( getConstPntrToComponent(valindex)->getRank()==1 ) {
267 159752 : ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer );
268 : return;
269 : }
270 229528 : const Value* myval=getConstPntrToComponent(valindex);
271 : unsigned ncols = myval->getNumberOfColumns(), matind = myval->getPositionInMatrixStash();
272 229528 : unsigned matbook_start = myval->getMatrixBookeepingStart(), vindex = bufstart + code*myval->getNumberOfColumns();
273 : const std::vector<unsigned> & matbook( myvals.getMatrixBookeeping() );
274 229528 : unsigned nelements = matbook[matbook_start+code*(1+ncols)];
275 229528 : if( ncols>=myval->getShape()[1] ) {
276 : // In this case we store the full matrix
277 6470140 : for(unsigned j=0; j<nelements; ++j) {
278 6275153 : unsigned jind = matbook[matbook_start+code*(1+ncols)+1+j];
279 : plumed_dbg_massert( vindex+j<buffer.size(), "failing in " + getLabel() + " on value " + myval->getName() );
280 6275153 : buffer[vindex + jind] += myvals.getStashedMatrixElement( matind, jind );
281 : }
282 : } else {
283 : // This is for storing sparse matrices when we can
284 965403 : for(unsigned j=0; j<nelements; ++j) {
285 930862 : unsigned jind = matbook[matbook_start+code*(1+ncols)+1+j];
286 : plumed_dbg_massert( vindex+j<buffer.size(), "failing in " + getLabel() + " on value " + myval->getName() );
287 930862 : buffer[vindex + j] += myvals.getStashedMatrixElement( matind, jind );
288 : }
289 : }
290 : }
291 :
292 721141 : bool ActionWithMatrix::checkForTaskForce( const unsigned& itask, const Value* myval ) const {
293 721141 : if( myval->getRank()<2 ) {
294 241252 : return ActionWithVector::checkForTaskForce( itask, myval );
295 : }
296 479889 : unsigned nelements = myval->getRowLength(itask), startr = itask*myval->getNumberOfColumns();
297 10439083 : for(unsigned j=0; j<nelements; ++j ) {
298 10051254 : if( fabs( myval->getForce( startr + j ) )>epsilon ) {
299 : return true;
300 : }
301 : }
302 : return false;
303 : }
304 :
305 210505 : void ActionWithMatrix::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
306 210505 : if( myval->getRank()==1 ) {
307 142549 : ActionWithVector::gatherForcesOnStoredValue( myval, itask, myvals, forces );
308 : return;
309 : }
310 : unsigned matind = myval->getPositionInMatrixStash();
311 842701414 : for(unsigned j=0; j<forces.size(); ++j) {
312 842633458 : forces[j] += myvals.getStashedMatrixForce( matind, j );
313 : }
314 : }
315 :
316 88131057 : void ActionWithMatrix::clearMatrixElements( MultiValue& myvals ) const {
317 88131057 : if( isActive() && clearOnEachCycle ) {
318 153880620 : for(int i=0; i<getNumberOfComponents(); ++i) {
319 103870566 : const Value* myval=getConstPntrToComponent(i);
320 103870566 : if( myval->getRank()==2 && !myval->hasDerivatives() ) {
321 103870566 : myvals.clearDerivatives( myval->getPositionInStream() );
322 : }
323 : }
324 : }
325 88131057 : if( matrix_to_do_after ) {
326 59005036 : matrix_to_do_after->clearMatrixElements( myvals );
327 : }
328 88131057 : }
329 :
330 : }
|
