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