Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2020 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 : #ifndef __PLUMED_function_FunctionOfMatrix_h
23 : #define __PLUMED_function_FunctionOfMatrix_h
24 :
25 : #include "core/ActionWithMatrix.h"
26 : #include "FunctionOfVector.h"
27 : #include "Sum.h"
28 : #include "tools/Matrix.h"
29 :
30 : namespace PLMD {
31 : namespace function {
32 :
33 : template <class T>
34 : class FunctionOfMatrix : public ActionWithMatrix {
35 : private:
36 : /// Is this the first step of the calculation
37 : bool firststep;
38 : /// The function that is being computed
39 : T myfunc;
40 : /// The number of derivatives for this action
41 : unsigned nderivatives;
42 : /// A vector that tells us if we have stored the input value
43 : std::vector<bool> stored_arguments;
44 : /// Switch off updating the arguments for this action
45 : std::vector<bool> update_arguments;
46 : /// The list of actiosn in this chain
47 : std::vector<std::string> actionsLabelsInChain;
48 : /// Get the shape of the output matrix
49 : std::vector<unsigned> getValueShapeFromArguments();
50 : public:
51 : static void registerKeywords(Keywords&);
52 : explicit FunctionOfMatrix(const ActionOptions&);
53 : /// Get the label to write in the graph
54 0 : std::string writeInGraph() const override { return myfunc.getGraphInfo( getName() ); }
55 : /// Make sure the derivatives are turned on
56 : void turnOnDerivatives() override;
57 : /// Get the number of derivatives for this action
58 : unsigned getNumberOfDerivatives() override ;
59 : /// Resize the matrices
60 : void prepare() override ;
61 : /// This gets the number of columns
62 : unsigned getNumberOfColumns() const override ;
63 : /// This checks for tasks in the parent class
64 : // void buildTaskListFromArgumentRequests( const unsigned& ntasks, bool& reduce, std::set<AtomNumber>& otasks ) override ;
65 : /// This ensures that we create some bookeeping stuff during the first step
66 : void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) override ;
67 : /// This sets up for the task
68 : void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
69 : /// Calculate the full matrix
70 : void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override ;
71 : /// This updates the indices for the matrix
72 : // void updateCentralMatrixIndex( const unsigned& ind, const std::vector<unsigned>& indices, MultiValue& myvals ) const override ;
73 : void runEndOfRowJobs( const unsigned& ind, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
74 : };
75 :
76 : template <class T>
77 1013 : void FunctionOfMatrix<T>::registerKeywords(Keywords& keys ) {
78 2026 : ActionWithMatrix::registerKeywords(keys); keys.use("ARG"); std::string name = keys.getDisplayName();
79 1013 : std::size_t und=name.find("_MATRIX"); keys.setDisplayName( name.substr(0,und) );
80 2026 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
81 2026 : keys.reserve("compulsory","PERIODIC","if the output of your function is periodic then you should specify the periodicity of the function. If the output is not periodic you must state this using PERIODIC=NO");
82 1013 : T tfunc; tfunc.registerKeywords( keys );
83 2026 : if( keys.getDisplayName()=="SUM" ) {
84 168 : keys.setValueDescription("the sum of all the elements in the input matrix");
85 1858 : } else if( keys.getDisplayName()=="HIGHEST" ) {
86 0 : keys.setValueDescription("the largest element of the input matrix");
87 1858 : } else if( keys.getDisplayName()=="LOWEST" ) {
88 0 : keys.setValueDescription("the smallest element in the input matrix");
89 1858 : } else if( keys.outputComponentExists(".#!value") ) {
90 1672 : keys.setValueDescription("the matrix obtained by doing an element-wise application of " + keys.getOutputComponentDescription(".#!value") + " to the input matrix");
91 : }
92 1899 : }
93 :
94 : template <class T>
95 495 : FunctionOfMatrix<T>::FunctionOfMatrix(const ActionOptions&ao):
96 : Action(ao),
97 : ActionWithMatrix(ao),
98 495 : firststep(true)
99 : {
100 451 : if( myfunc.getArgStart()>0 ) error("this has not beeen implemented -- if you are interested email gareth.tribello@gmail.com");
101 : // Get the shape of the output
102 495 : std::vector<unsigned> shape( getValueShapeFromArguments() );
103 : // Check if the output matrix is symmetric
104 495 : bool symmetric=true; unsigned argstart=myfunc.getArgStart();
105 1508 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
106 1013 : if( getPntrToArgument(i)->getRank()==2 ) {
107 948 : if( !getPntrToArgument(i)->isSymmetric() ) { symmetric=false; }
108 : }
109 : }
110 : // Read the input and do some checks
111 495 : myfunc.read( this );
112 : // Setup to do this in chain if possible
113 495 : if( myfunc.doWithTasks() ) done_in_chain=true;
114 : // Check we are not calculating a sum
115 41 : if( myfunc.zeroRank() ) shape.resize(0);
116 : // Get the names of the components
117 495 : std::vector<std::string> components( keywords.getOutputComponents() );
118 : // Create the values to hold the output
119 42 : std::vector<std::string> str_ind( myfunc.getComponentsPerLabel() );
120 1034 : for(unsigned i=0; i<components.size(); ++i) {
121 84 : if( str_ind.size()>0 ) {
122 168 : std::string compstr = components[i]; if( components[i]==".#!value" ) compstr = "";
123 760 : for(unsigned j=0; j<str_ind.size(); ++j) {
124 : if( myfunc.zeroRank() ) {
125 : addComponentWithDerivatives( compstr + str_ind[j], shape );
126 : } else {
127 1352 : addComponent( compstr + str_ind[j], shape );
128 676 : getPntrToComponent(i*str_ind.size()+j)->setSymmetric( symmetric );
129 : }
130 : }
131 41 : } else if( components[i]==".#!value" && myfunc.zeroRank() ) {
132 41 : addValueWithDerivatives( shape );
133 414 : } else if( components[i]==".#!value" ) {
134 410 : addValue( shape ); getPntrToComponent(0)->setSymmetric( symmetric );
135 4 : } else if( components[i].find_first_of("_")!=std::string::npos ) {
136 0 : if( getNumberOfArguments()-argstart==1 ) { addValue( shape ); getPntrToComponent(0)->setSymmetric( symmetric ); }
137 : else {
138 0 : for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
139 0 : addComponent( getPntrToArgument(j)->getName() + components[i], shape );
140 0 : getPntrToComponent(i*(getNumberOfArguments()-argstart)+j-argstart)->setSymmetric( symmetric );
141 : }
142 : }
143 4 : } else { addComponent( components[i], shape ); getPntrToComponent(i)->setSymmetric( symmetric ); }
144 : }
145 : // Check if this can be sped up
146 370 : if( myfunc.getDerivativeZeroIfValueIsZero() ) {
147 174 : for(int i=0; i<getNumberOfComponents(); ++i) getPntrToComponent(i)->setDerivativeIsZeroWhenValueIsZero();
148 : }
149 : // Set the periodicities of the output components
150 495 : myfunc.setPeriodicityForOutputs( this );
151 : // We can't do this with if we are dividing a stack by some a product v.v^T product as we need to store the vector
152 : // In order to do this type of calculation. There should be a neater fix than this but I can't see it.
153 : bool foundneigh=false; const ActionWithMatrix* chainstart = NULL;
154 1503 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
155 1207 : if( getPntrToArgument(i)->isConstant() && getNumberOfArguments()>1 ) continue ;
156 932 : std::string argname=(getPntrToArgument(i)->getPntrToAction())->getName();
157 932 : if( argname=="NEIGHBORS" ) { foundneigh=true; break; }
158 929 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
159 929 : if( !av ) done_in_chain=false;
160 929 : if( getPntrToArgument(i)->getRank()==0 ) {
161 0 : function::FunctionOfVector<function::Sum>* as = dynamic_cast<function::FunctionOfVector<function::Sum>*>( getPntrToArgument(i)->getPntrToAction() );
162 0 : if(as) done_in_chain=false;
163 929 : } else if( getPntrToArgument(i)->ignoreStoredValue( getLabel() ) ) {
164 : // This option deals with the case when you have two adjacency matrices, A_ij and B_ij, multiplied together. This cannot be done in the chain as the rows
165 : // of the two adjacency matrix are run over separately. The value A_ij is thus not available when B_ij is calculated.
166 853 : ActionWithMatrix* am = dynamic_cast<ActionWithMatrix*>( getPntrToArgument(i)->getPntrToAction() );
167 853 : plumed_assert( am ); const ActionWithMatrix* thischain = am->getFirstMatrixInChain();
168 853 : if( !thischain->isAdjacencyMatrix() && thischain->getName()!="VSTACK" ) continue;
169 657 : if( !chainstart ) chainstart = thischain;
170 317 : else if( thischain!=chainstart ) done_in_chain=false;
171 : }
172 : }
173 : // If we are working with neighbors we trick PLUMED into storing ALL the components of the other arguments
174 : // in this way we can ensure that the function of the neighbours matrix is in a chain starting from the
175 : // Neighbours matrix action.
176 : if( foundneigh ) {
177 9 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
178 6 : ActionWithValue* av=getPntrToArgument(i)->getPntrToAction();
179 6 : if( av->getName()!="NEIGHBORS" ) {
180 8 : for(int i=0; i<av->getNumberOfComponents(); ++i) (av->copyOutput(i))->buildDataStore();
181 : }
182 : }
183 : }
184 : // Now setup the action in the chain if we can
185 495 : nderivatives = buildArgumentStore(myfunc.getArgStart());
186 990 : }
187 :
188 : template <class T>
189 1921 : void FunctionOfMatrix<T>::turnOnDerivatives() {
190 1921 : if( !myfunc.derivativesImplemented() ) error("derivatives have not been implemended for " + getName() );
191 1921 : ActionWithValue::turnOnDerivatives(); myfunc.setup(this);
192 1921 : }
193 :
194 : template <class T>
195 30411 : unsigned FunctionOfMatrix<T>::getNumberOfDerivatives() {
196 30411 : return nderivatives;
197 : }
198 :
199 : template <class T>
200 2229 : void FunctionOfMatrix<T>::prepare() {
201 2229 : unsigned argstart = myfunc.getArgStart(); std::vector<unsigned> shape(2);
202 2229 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
203 2229 : if( getPntrToArgument(i)->getRank()==2 ) {
204 2229 : shape[0] = getPntrToArgument(i)->getShape()[0];
205 2229 : shape[1] = getPntrToArgument(i)->getShape()[1];
206 2229 : break;
207 : }
208 : }
209 6682 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
210 4453 : Value* myval = getPntrToComponent(i);
211 4453 : if( myval->getRank()==2 && (myval->getShape()[0]!=shape[0] || myval->getShape()[1]!=shape[1]) ) {
212 18 : myval->setShape(shape); if( myval->valueIsStored() ) myval->reshapeMatrixStore( shape[1] );
213 : }
214 : }
215 2229 : ActionWithVector::prepare();
216 2229 : }
217 :
218 : template <class T>
219 281844 : unsigned FunctionOfMatrix<T>::getNumberOfColumns() const {
220 281844 : if( getConstPntrToComponent(0)->getRank()==2 ) {
221 : unsigned argstart=myfunc.getArgStart();
222 281844 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
223 281844 : if( getPntrToArgument(i)->getRank()==2 ) {
224 281844 : ActionWithMatrix* am=dynamic_cast<ActionWithMatrix*>( getPntrToArgument(i)->getPntrToAction() );
225 281844 : if( am ) return am->getNumberOfColumns();
226 2238 : return getPntrToArgument(i)->getShape()[1];
227 : }
228 : }
229 : }
230 0 : plumed_error(); return 0;
231 : }
232 :
233 : template <class T>
234 4209 : void FunctionOfMatrix<T>::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
235 11667 : for(unsigned i=0; i<getNumberOfArguments(); ++i) plumed_assert( getPntrToArgument(i)->getRank()==2 );
236 4209 : unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(0)->getShape()[1];
237 4209 : if( indices.size()!=size_v+1 ) indices.resize( size_v+1 );
238 642613 : for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + i;
239 : myvals.setSplitIndex( size_v + 1 );
240 4209 : }
241 :
242 : // template <class T>
243 : // void FunctionOfMatrix<T>::buildTaskListFromArgumentRequests( const unsigned& ntasks, bool& reduce, std::set<AtomNumber>& otasks ) {
244 : // // Check if this is the first element in a chain
245 : // if( actionInChain() ) return;
246 : // // If it is computed outside a chain get the tassks the daughter chain needs
247 : // propegateTaskListsForValue( 0, ntasks, reduce, otasks );
248 : // }
249 :
250 : template <class T>
251 2525 : void FunctionOfMatrix<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
252 2525 : if( firststep ) {
253 489 : stored_arguments.resize( getNumberOfArguments() );
254 489 : update_arguments.resize( getNumberOfArguments(), true );
255 489 : std::string control = getFirstActionInChain()->getLabel();
256 1484 : for(unsigned i=0; i<stored_arguments.size(); ++i) {
257 995 : stored_arguments[i] = !getPntrToArgument(i)->ignoreStoredValue( control );
258 995 : if( !stored_arguments[i] ) update_arguments[i] = true;
259 164 : else update_arguments[i] = !argumentDependsOn( headstr, this, getPntrToArgument(i) );
260 : }
261 489 : firststep=false;
262 : }
263 2525 : ActionWithMatrix::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
264 2525 : }
265 :
266 : template <class T>
267 27928651 : void FunctionOfMatrix<T>::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
268 27928651 : unsigned argstart=myfunc.getArgStart(); std::vector<double> args( getNumberOfArguments() - argstart );
269 27928651 : unsigned ind2 = index2;
270 27928651 : if( getConstPntrToComponent(0)->getRank()==2 && index2>=getConstPntrToComponent(0)->getShape()[0] ) ind2 = index2 - getConstPntrToComponent(0)->getShape()[0];
271 24292268 : else if( index2>=getPntrToArgument(0)->getShape()[0] ) ind2 = index2 - getPntrToArgument(0)->getShape()[0];
272 27928651 : if( actionInChain() ) {
273 85619946 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
274 58329699 : if( getPntrToArgument(i)->getRank()==0 ) args[i-argstart] = getPntrToArgument(i)->get();
275 58193979 : else if( !getPntrToArgument(i)->valueHasBeenSet() ) args[i-argstart] = myvals.get( getPntrToArgument(i)->getPositionInStream() );
276 1188593 : else args[i-argstart] = getPntrToArgument(i)->get( getPntrToArgument(i)->getShape()[1]*index1 + ind2 );
277 : }
278 : } else {
279 1727072 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
280 1088668 : if( getPntrToArgument(i)->getRank()==2 ) args[i-argstart]=getPntrToArgument(i)->get( getPntrToArgument(i)->getShape()[1]*index1 + ind2 );
281 0 : else args[i-argstart] = getPntrToArgument(i)->get();
282 : }
283 : }
284 : // Calculate the function and its derivatives
285 27928651 : std::vector<double> vals( getNumberOfComponents() ); Matrix<double> derivatives( getNumberOfComponents(), getNumberOfArguments()-argstart );
286 27928651 : myfunc.calc( this, args, vals, derivatives );
287 : // And set the values
288 99634355 : for(unsigned i=0; i<vals.size(); ++i) myvals.addValue( getConstPntrToComponent(i)->getPositionInStream(), vals[i] );
289 : // Return if we are not computing derivatives
290 27928651 : if( doNotCalculateDerivatives() ) return;
291 :
292 5399619 : if( actionInChain() ) {
293 33385311 : for(int i=0; i<getNumberOfComponents(); ++i) {
294 27990552 : unsigned ostrn=getConstPntrToComponent(i)->getPositionInStream();
295 131996523 : for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
296 104005971 : if( getPntrToArgument(j)->getRank()==2 ) {
297 103890411 : unsigned istrn = getPntrToArgument(j)->getPositionInStream();
298 103890411 : if( stored_arguments[j] ) {
299 395048 : unsigned task_index = getPntrToArgument(i)->getShape()[1]*index1 + ind2;
300 395048 : myvals.clearDerivatives(istrn); myvals.addDerivative( istrn, task_index, 1.0 ); myvals.updateIndex( istrn, task_index );
301 : }
302 470717695 : for(unsigned k=0; k<myvals.getNumberActive(istrn); ++k) {
303 366827284 : unsigned kind=myvals.getActiveIndex(istrn,k);
304 366827284 : myvals.addDerivative( ostrn, arg_deriv_starts[j] + kind, derivatives(i,j)*myvals.getDerivative( istrn, kind ) );
305 : }
306 : }
307 : }
308 : }
309 : // If we are computing a matrix we need to update the indices here so that derivatives are calcualted correctly in functions of these
310 5394759 : if( getConstPntrToComponent(0)->getRank()==2 ) {
311 32784527 : for(int i=0; i<getNumberOfComponents(); ++i) {
312 27690160 : unsigned ostrn=getConstPntrToComponent(i)->getPositionInStream();
313 131395739 : for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
314 103705579 : if( !update_arguments[j] || getPntrToArgument(j)->getRank()==0 ) continue ;
315 : // Ensure we only store one lot of derivative indices
316 : bool found=false;
317 105009550 : for(unsigned k=0; k<j; ++k) {
318 76601003 : if( arg_deriv_starts[k]==arg_deriv_starts[j] ) { found=true; break; }
319 : }
320 103589995 : if( found ) continue;
321 : unsigned istrn = getPntrToArgument(j)->getPositionInStream();
322 138447375 : for(unsigned k=0; k<myvals.getNumberActive(istrn); ++k) {
323 110038828 : unsigned kind=myvals.getActiveIndex(istrn,k);
324 110038828 : myvals.updateIndex( ostrn, arg_deriv_starts[j] + kind );
325 : }
326 : }
327 : }
328 : }
329 : } else {
330 4860 : unsigned base=0; ind2 = index2;
331 4860 : for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
332 4860 : if( getPntrToArgument(j)->getRank()!=2 ) continue ;
333 4860 : if( index2>=getPntrToArgument(j)->getShape()[0] ) ind2 = index2 - getPntrToArgument(j)->getShape()[0];
334 : break;
335 : }
336 13965 : for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
337 9105 : if( getPntrToArgument(j)->getRank()==2 ) {
338 18210 : for(int i=0; i<getNumberOfComponents(); ++i) {
339 9105 : unsigned ostrn=getConstPntrToComponent(i)->getPositionInStream();
340 9105 : unsigned myind = base + getPntrToArgument(j)->getShape()[1]*index1 + ind2;
341 9105 : myvals.addDerivative( ostrn, myind, derivatives(i,j) );
342 9105 : myvals.updateIndex( ostrn, myind );
343 : }
344 : } else {
345 0 : for(int i=0; i<getNumberOfComponents(); ++i) {
346 0 : unsigned ostrn=getConstPntrToComponent(i)->getPositionInStream();
347 0 : myvals.addDerivative( ostrn, base, derivatives(i,j) );
348 0 : myvals.updateIndex( ostrn, base );
349 : }
350 : }
351 9105 : base += getPntrToArgument(j)->getNumberOfValues();
352 : }
353 : }
354 : }
355 :
356 : template <class T>
357 226389 : void FunctionOfMatrix<T>::runEndOfRowJobs( const unsigned& ind, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
358 226389 : if( doNotCalculateDerivatives() ) return;
359 :
360 : unsigned argstart=myfunc.getArgStart();
361 71237 : if( actionInChain() && getConstPntrToComponent(0)->getRank()==2 ) {
362 : // This is triggered if we are outputting a matrix
363 624578 : for(int vv=0; vv<getNumberOfComponents(); ++vv) {
364 558183 : unsigned nmat = getConstPntrToComponent(vv)->getPositionInMatrixStash();
365 : std::vector<unsigned>& mat_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); unsigned ntot_mat=0;
366 558183 : if( mat_indices.size()<nderivatives ) mat_indices.resize( nderivatives );
367 2701544 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
368 2143361 : if( !update_arguments[i] || getPntrToArgument(i)->getRank()==0 ) continue ;
369 : // Ensure we only store one lot of derivative indices
370 : bool found=false;
371 2168017 : for(unsigned j=0; j<i; ++j) {
372 1591516 : if( arg_deriv_starts[j]==arg_deriv_starts[i] ) { found=true; break; }
373 : }
374 2142277 : if( found ) continue;
375 :
376 576501 : if( stored_arguments[i] ) {
377 15483 : unsigned tbase = getPntrToArgument(i)->getShape()[1]*ind;
378 410507 : for(unsigned k=1; k<indices.size(); ++k) {
379 395024 : unsigned ind2 = indices[k] - getConstPntrToComponent(0)->getShape()[0];
380 395024 : mat_indices[ntot_mat + k - 1] = arg_deriv_starts[i] + tbase + ind2;
381 : }
382 15483 : ntot_mat += indices.size()-1;
383 : } else {
384 : unsigned istrn = getPntrToArgument(i)->getPositionInMatrixStash();
385 : std::vector<unsigned>& imat_indices( myvals.getMatrixRowDerivativeIndices( istrn ) );
386 31848426 : for(unsigned k=0; k<myvals.getNumberOfMatrixRowDerivatives( istrn ); ++k) mat_indices[ntot_mat + k] = arg_deriv_starts[i] + imat_indices[k];
387 561018 : ntot_mat += myvals.getNumberOfMatrixRowDerivatives( istrn );
388 : }
389 : }
390 : myvals.setNumberOfMatrixRowDerivatives( nmat, ntot_mat );
391 : }
392 4842 : } else if( actionInChain() ) {
393 : // This is triggered if we are calculating a single scalar in the function
394 8822 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
395 : bool found=false;
396 4411 : for(unsigned j=0; j<i; ++j) {
397 0 : if( arg_deriv_starts[j]==arg_deriv_starts[i] ) { found=true; break; }
398 : }
399 4411 : if( found ) continue;
400 : unsigned istrn = getPntrToArgument(i)->getPositionInMatrixStash();
401 : std::vector<unsigned>& mat_indices( myvals.getMatrixRowDerivativeIndices( istrn ) );
402 926766 : for(unsigned k=0; k<myvals.getNumberOfMatrixRowDerivatives( istrn ); ++k) {
403 1844710 : for(int j=0; j<getNumberOfComponents(); ++j) {
404 922355 : unsigned ostrn = getConstPntrToComponent(j)->getPositionInStream();
405 922355 : myvals.updateIndex( ostrn, arg_deriv_starts[i] + mat_indices[k] );
406 : }
407 : }
408 : }
409 431 : } else if( getConstPntrToComponent(0)->getRank()==2 ) {
410 760 : for(int vv=0; vv<getNumberOfComponents(); ++vv) {
411 380 : unsigned nmat = getConstPntrToComponent(vv)->getPositionInMatrixStash();
412 : std::vector<unsigned>& mat_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); unsigned ntot_mat=0;
413 380 : if( mat_indices.size()<nderivatives ) mat_indices.resize( nderivatives ); unsigned matderbase = 0;
414 986 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
415 606 : if( getPntrToArgument(i)->getRank()==0 ) continue ;
416 606 : unsigned ss = getPntrToArgument(i)->getShape()[1]; unsigned tbase = matderbase + ss*myvals.getTaskIndex();
417 9558 : for(unsigned k=0; k<ss; ++k) mat_indices[ntot_mat + k] = tbase + k;
418 606 : ntot_mat += ss; matderbase += getPntrToArgument(i)->getNumberOfValues();
419 : }
420 : myvals.setNumberOfMatrixRowDerivatives( nmat, ntot_mat );
421 : }
422 : }
423 : }
424 :
425 : template <class T>
426 495 : std::vector<unsigned> FunctionOfMatrix<T>::getValueShapeFromArguments() {
427 495 : unsigned argstart=myfunc.getArgStart(); std::vector<unsigned> shape(2); shape[0]=shape[1]=0;
428 1508 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
429 1013 : plumed_assert( getPntrToArgument(i)->getRank()==2 || getPntrToArgument(i)->getRank()==0 );
430 1013 : if( getPntrToArgument(i)->getRank()==2 ) {
431 948 : if( shape[0]>0 && (getPntrToArgument(i)->getShape()[0]!=shape[0] || getPntrToArgument(i)->getShape()[1]!=shape[1]) ) error("all matrices input should have the same shape");
432 948 : else if( shape[0]==0 ) { shape[0]=getPntrToArgument(i)->getShape()[0]; shape[1]=getPntrToArgument(i)->getShape()[1]; }
433 948 : plumed_assert( !getPntrToArgument(i)->hasDerivatives() );
434 : }
435 : }
436 495 : myfunc.setPrefactor( this, 1.0 ); return shape;
437 : }
438 :
439 : }
440 : }
441 : #endif
|
