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 3256 : void ActionWithMatrix::registerKeywords( Keywords& keys ) {
28 3256 : ActionWithVector::registerKeywords( keys );
29 3256 : keys.use("ARG");
30 3256 : }
31 :
32 1482 : ActionWithMatrix::ActionWithMatrix(const ActionOptions&ao):
33 : Action(ao),
34 : ActionWithVector(ao),
35 1482 : next_action_in_chain(NULL),
36 1482 : matrix_to_do_before(NULL),
37 1482 : matrix_to_do_after(NULL),
38 1482 : clearOnEachCycle(true) {
39 1482 : }
40 :
41 1482 : ActionWithMatrix::~ActionWithMatrix() {
42 1482 : if( matrix_to_do_before ) {
43 830 : matrix_to_do_before->matrix_to_do_after=NULL;
44 830 : matrix_to_do_before->next_action_in_chain=NULL;
45 : }
46 1482 : }
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 56158 : void ActionWithMatrix::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
64 56158 : ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
65 :
66 122188 : for(int i=0; i<getNumberOfComponents(); ++i) {
67 66030 : Value* myval=getPntrToComponent(i);
68 66030 : if( myval->getRank()!=2 || myval->hasDerivatives() ) {
69 30514 : continue;
70 : }
71 35516 : myval->setPositionInMatrixStash(nmat);
72 35516 : nmat++;
73 35516 : if( !myval->valueIsStored() ) {
74 27718 : continue;
75 : }
76 7798 : if( myval->getShape()[1]>maxcol ) {
77 7137 : maxcol=myval->getShape()[1];
78 : }
79 : myval->setMatrixBookeepingStart(nbookeeping);
80 7798 : 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 56158 : clearOnEachCycle = false;
84 86672 : for(int i=0; i<getNumberOfComponents(); ++i) {
85 61216 : const Value* myval=getConstPntrToComponent(i);
86 61216 : if( myval->getRank()==2 && !myval->hasDerivatives() ) {
87 30702 : clearOnEachCycle = true;
88 30702 : break;
89 : }
90 : }
91 : // Turn off clearing of derivatives if we have only the values of adjacency matrices
92 56158 : if( doNotCalculateDerivatives() && isAdjacencyMatrix() ) {
93 852 : clearOnEachCycle = false;
94 : }
95 56158 : }
96 :
97 3178 : void ActionWithMatrix::finishChainBuild( ActionWithVector* act ) {
98 3178 : ActionWithMatrix* am=dynamic_cast<ActionWithMatrix*>(act);
99 3178 : 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 2267 : if( next_action_in_chain ) {
104 1380 : 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 : return ;
110 : }
111 830 : plumed_massert( !matrix_to_do_after, "cannot add " + act->getLabel() + " in " + getLabel() + " as have already added " + matrix_to_do_after->getLabel() );
112 830 : matrix_to_do_after=am;
113 830 : 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 30978 : void ActionWithMatrix::getTotalMatrixBookeeping( unsigned& nbookeeping ) {
125 67937 : for(int i=0; i<getNumberOfComponents(); ++i) {
126 36959 : Value* myval=getPntrToComponent(i);
127 36959 : if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() ) {
128 32174 : continue;
129 : }
130 4785 : myval->reshapeMatrixStore( getNumberOfColumns() );
131 4785 : nbookeeping += myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
132 : }
133 30978 : if( next_action_in_chain ) {
134 13479 : next_action_in_chain->getTotalMatrixBookeeping( nbookeeping );
135 : }
136 30978 : }
137 :
138 31020 : void ActionWithMatrix::calculate() {
139 31020 : if( actionInChain() ) {
140 13521 : return ;
141 : }
142 : // Update all the neighbour lists
143 17499 : updateAllNeighbourLists();
144 : // Setup the matrix indices
145 17499 : unsigned nbookeeping=0;
146 17499 : getTotalMatrixBookeeping( nbookeeping );
147 17499 : if( matrix_bookeeping.size()!=nbookeeping ) {
148 380 : matrix_bookeeping.resize( nbookeeping );
149 : }
150 : std::fill( matrix_bookeeping.begin(), matrix_bookeeping.end(), 0 );
151 : // And run all the tasks
152 17499 : runAllTasks();
153 : }
154 :
155 30978 : void ActionWithMatrix::updateAllNeighbourLists() {
156 30978 : updateNeighbourList();
157 30978 : if( next_action_in_chain ) {
158 13479 : next_action_in_chain->updateAllNeighbourLists();
159 : }
160 30978 : }
161 :
162 2404351 : void ActionWithMatrix::performTask( const unsigned& task_index, MultiValue& myvals ) const {
163 : std::vector<unsigned> & indices( myvals.getIndices() );
164 2404351 : if( matrix_to_do_before ) {
165 : plumed_dbg_assert( myvals.inVectorCall() );
166 1314856 : runEndOfRowJobs( task_index, indices, myvals );
167 1314856 : return;
168 : }
169 1089495 : setupForTask( task_index, indices, myvals );
170 :
171 : // Now loop over the row of the matrix
172 : unsigned ntwo_atoms = myvals.getSplitIndex();
173 44938581 : for(unsigned i=1; i<ntwo_atoms; ++i) {
174 : // This does everything in the stream that is done with single matrix elements
175 43849086 : runTask( getLabel(), task_index, indices[i], myvals );
176 : // Now clear only elements that are not accumulated over whole row
177 43849086 : 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 1089495 : runEndOfRowJobs( task_index, indices, myvals );
181 : }
182 :
183 132572654 : void ActionWithMatrix::runTask( const std::string& controller, const unsigned& current, const unsigned colno, MultiValue& myvals ) const {
184 : double outval=0;
185 132572654 : myvals.setTaskIndex(current);
186 132572654 : myvals.setSecondTaskIndex( colno );
187 132572654 : if( isActive() ) {
188 124764100 : performTask( controller, current, colno, myvals );
189 : }
190 132572654 : bool hasval = !isAdjacencyMatrix();
191 344330772 : for(int i=0; i<getNumberOfComponents(); ++i) {
192 292594234 : if( fabs(myvals.get( getConstPntrToComponent(i)->getPositionInStream()) )>0 ) {
193 : hasval=true;
194 : break;
195 : }
196 : }
197 :
198 132572654 : if( hasval ) {
199 521241818 : for(int i=0; i<getNumberOfComponents(); ++i) {
200 402540890 : const Value* myval=getConstPntrToComponent(i);
201 402540890 : if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() ) {
202 384052058 : continue;
203 : }
204 18488832 : unsigned matindex = myval->getPositionInMatrixStash(), matbook_start = myval->getMatrixBookeepingStart(), col_stash_index = colno;
205 18488832 : if( colno>=myval->getShape()[0] ) {
206 11949150 : col_stash_index = colno - myval->getShape()[0];
207 : }
208 18488832 : unsigned rowstart = matbook_start+current*(1+myval->getNumberOfColumns());
209 18488832 : 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 21864543 : for(unsigned j=0; j<myvals.getNumberActive(sind); ++j) {
216 20313477 : unsigned kindex = myvals.getActiveIndex(sind,j);
217 20313477 : myvals.addMatrixForce( matindex, kindex, fforce*myvals.getDerivative(sind,kindex ) );
218 : }
219 : }
220 18488832 : double finalval = myvals.get( myval->getPositionInStream() );
221 18488832 : if( fabs(finalval)>0 ) {
222 12751774 : myvals.stashMatrixElement( matindex, rowstart, col_stash_index, finalval );
223 : }
224 : }
225 : }
226 132572654 : if( matrix_to_do_after ) {
227 88723568 : matrix_to_do_after->runTask( controller, current, colno, myvals );
228 : }
229 132572654 : }
230 :
231 20037 : void ActionWithMatrix::gatherThreads( const unsigned& nt, const unsigned& bufsize, const std::vector<double>& omp_buffer, std::vector<double>& buffer, MultiValue& myvals ) {
232 20037 : ActionWithVector::gatherThreads( nt, bufsize, omp_buffer, buffer, myvals );
233 64708733 : for(unsigned i=0; i<matrix_bookeeping.size(); ++i) {
234 64688696 : matrix_bookeeping[i] += myvals.getMatrixBookeeping()[i];
235 : }
236 20037 : }
237 :
238 17499 : void ActionWithMatrix::gatherProcesses( std::vector<double>& buffer ) {
239 17499 : ActionWithVector::gatherProcesses( buffer );
240 17499 : if( matrix_bookeeping.size()>0 && !runInSerial() ) {
241 4251 : comm.Sum( matrix_bookeeping );
242 : }
243 17499 : unsigned nval=0;
244 17499 : transferNonZeroMatrixElementsToValues( nval, matrix_bookeeping );
245 17499 : }
246 :
247 30978 : void ActionWithMatrix::transferNonZeroMatrixElementsToValues( unsigned& nval, const std::vector<unsigned>& matbook ) {
248 67937 : for(int i=0; i<getNumberOfComponents(); ++i) {
249 36959 : Value* myval=getPntrToComponent(i);
250 36959 : if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() || getNumberOfColumns()>=myval->getShape()[1] ) {
251 36860 : 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 30978 : if( next_action_in_chain ) {
260 13479 : next_action_in_chain->transferNonZeroMatrixElementsToValues( nval, matbook );
261 : }
262 30978 : }
263 :
264 731298 : void ActionWithMatrix::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
265 : const unsigned& bufstart, std::vector<double>& buffer ) const {
266 731298 : if( getConstPntrToComponent(valindex)->getRank()==1 ) {
267 335650 : ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer );
268 : return;
269 : }
270 395648 : const Value* myval=getConstPntrToComponent(valindex);
271 : unsigned ncols = myval->getNumberOfColumns(), matind = myval->getPositionInMatrixStash();
272 395648 : unsigned matbook_start = myval->getMatrixBookeepingStart(), vindex = bufstart + code*myval->getNumberOfColumns();
273 : const std::vector<unsigned> & matbook( myvals.getMatrixBookeeping() );
274 395648 : unsigned nelements = matbook[matbook_start+code*(1+ncols)];
275 395648 : if( ncols>=myval->getShape()[1] ) {
276 : // In this case we store the full matrix
277 10789164 : for(unsigned j=0; j<nelements; ++j) {
278 10428057 : 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 10428057 : 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 896043 : bool ActionWithMatrix::checkForTaskForce( const unsigned& itask, const Value* myval ) const {
293 896043 : if( myval->getRank()<2 ) {
294 416154 : 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 212846 : void ActionWithMatrix::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
306 212846 : if( myval->getRank()==1 ) {
307 144890 : 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 132572654 : void ActionWithMatrix::clearMatrixElements( MultiValue& myvals ) const {
317 132572654 : if( isActive() && clearOnEachCycle ) {
318 300709119 : for(int i=0; i<getNumberOfComponents(); ++i) {
319 223404104 : const Value* myval=getConstPntrToComponent(i);
320 223404104 : if( myval->getRank()==2 && !myval->hasDerivatives() ) {
321 223404104 : myvals.clearDerivatives( myval->getPositionInStream() );
322 : }
323 : }
324 : }
325 132572654 : if( matrix_to_do_after ) {
326 88723568 : matrix_to_do_after->clearMatrixElements( myvals );
327 : }
328 132572654 : }
329 :
330 : }
|
