Merge branch 'master' into derivatives-from-backpropegation

src/adjmat/AdjacencyMatrixBase.cpp

@@ -52,6 +52,7 @@ AdjacencyMatrixBase::AdjacencyMatrixBase(const ActionOptions& ao):
   neighbour_list_updated(false),
   linkcells(comm),
   threecells(comm),
+  maxcol(0),
   natoms_per_list(0)
 {
   std::vector<unsigned> shape(2); std::vector<AtomNumber> t; parseAtomList("GROUP", t );
@@ -291,9 +292,9 @@ void AdjacencyMatrixBase::setupForTask( const unsigned& current, std::vector<uns
 void AdjacencyMatrixBase::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
   Vector zero; zero.zero(); plumed_dbg_assert( index2<myvals.getAtomVector().size() );
   double weight = calculateWeight( zero, myvals.getAtomVector()[index2], myvals.getNumberOfIndices()-myvals.getSplitIndex(), myvals );
-  if( fabs(weight)<epsilon ) return;
-
   unsigned w_ind = getConstPntrToComponent(0)->getPositionInStream(); myvals.setValue( w_ind, weight );
+  if( fabs(weight)<epsilon ) { myvals.setValue( w_ind, 0 ); return; }
+
   if( !doNotCalculateDerivatives() ) {
     // Update dynamic list indices for central atom
     myvals.updateIndex( w_ind, 3*index1+0 ); myvals.updateIndex( w_ind, 3*index1+1 ); myvals.updateIndex( w_ind, 3*index1+2 );

src/adjmat/ContactMatrix.cpp

@@ -113,6 +113,7 @@ double ContactMatrix::calculateWeight( const Vector& pos1, const Vector& pos2, c
   if( mod2<epsilon ) return 0.0;  // Atoms can't be bonded to themselves
   double dfunc, val = switchingFunction.calculateSqr( mod2, dfunc );
   if( val<epsilon ) return 0.0;
+  if( doNotCalculateDerivatives() ) return val;
   addAtomDerivatives( 0, (-dfunc)*distance, myvals );
   addAtomDerivatives( 1, (+dfunc)*distance, myvals );
   addBoxDerivatives( (-dfunc)*Tensor(distance,distance), myvals );

src/cltools/Benchmark.cpp

@@ -484,7 +484,7 @@ int Benchmark::main(FILE* in, FILE*out,Communicator& pc) {
     log.link(log_dev_null.get());
   }
   log.setLinePrefix("BENCH:  ");
-
+  log <<"Welcome to PLUMED benchmark\n";
   std::vector<Kernel> kernels;
 
   // perform comparative analysis
@@ -589,13 +589,15 @@ int Benchmark::main(FILE* in, FILE*out,Communicator& pc) {
     {
       std::string paths;
       parse("--kernel",paths);
+      log <<"Using --kernel=" << paths << "\n";
       allpaths=Tools::getWords(paths,":");
     }
 
     std::vector<std::string> allplumed;
     {
       std::string paths;
       parse("--plumed",paths);
+      log <<"Using --plumed=" << paths << "\n";
       allplumed=Tools::getWords(paths,":");
     }
 
@@ -628,27 +630,36 @@ int Benchmark::main(FILE* in, FILE*out,Communicator& pc) {
   // read other flags:
   bool shuffled=false;
   parseFlag("--shuffled",shuffled);
+  if (shuffled)
+    log << "Using --shuffled\n";
   int nf; parse("--nsteps",nf);
+  log << "Using --nsteps=" << nf << "\n";
   unsigned natoms; parse("--natoms",natoms);
-
+  log << "Using --natoms=" << natoms << "\n";
   double maxtime; parse("--maxtime",maxtime);
+  log << "Using --maxtime=" << maxtime << "\n";
 
   bool domain_decomposition=false;
   parseFlag("--domain-decomposition",domain_decomposition);
+  if (domain_decomposition)
+    log << "Using --domain-decomposition\n";
+
   if(pc.Get_size()>1) domain_decomposition=true;
   if(domain_decomposition) shuffled=true;
 
   double timeToSleep;
   parse("--sleep",timeToSleep);
+  log << "Using --sleep=" << timeToSleep << "\n";
 
   std::vector<int> shuffled_indexes;
 
   {
     std::string atomicDistr;
     parse("--atom-distribution",atomicDistr);
     distribution = getAtomDistribution(atomicDistr);
+    log << "Using --atom-distribution=" << atomicDistr << "\n";
   }
-
+  log <<"Initializing the setup of the kernel(s)\n";
   const auto initial_time=std::chrono::high_resolution_clock::now();
 
   for(auto & k : kernels) {
@@ -690,7 +701,6 @@ int Benchmark::main(FILE* in, FILE*out,Communicator& pc) {
   // trap signals:
   SignalHandlerGuard sigIntGuard(SIGINT, signalHandler);
 
-
   for(int step=0; nf<0 || step<nf; ++step) {
     std::shuffle(kernels_ptr.begin(),kernels_ptr.end(),rng);
     distribution->positions(pos,step,atomicGenerator);

src/core/ActionWithMatrix.cpp

@@ -33,7 +33,8 @@ ActionWithMatrix::ActionWithMatrix(const ActionOptions&ao):
   ActionWithVector(ao),
   next_action_in_chain(NULL),
   matrix_to_do_before(NULL),
-  matrix_to_do_after(NULL)
+  matrix_to_do_after(NULL),
+  clearOnEachCycle(true)
 {
 }
 
@@ -62,6 +63,14 @@ void ActionWithMatrix::setupStreamedComponents( const std::string& headstr, unsi
     myval->setMatrixBookeepingStart(nbookeeping);
     nbookeeping += myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
   }
+  // Turn off clearning of derivatives after each matrix run if there are no matrices in the output of this action
+  clearOnEachCycle = false;
+  for(int i=0; i<getNumberOfComponents(); ++i) {
+    const Value* myval=getConstPntrToComponent(i);
+    if( myval->getRank()==2 && !myval->hasDerivatives() ) { clearOnEachCycle = true; break; }
+  }
+  // Turn off clearing of derivatives if we have only the values of adjacency matrices
+  if( doNotCalculateDerivatives() && isAdjacencyMatrix() ) clearOnEachCycle = false;
 }
 
 void ActionWithMatrix::finishChainBuild( ActionWithVector* act ) {
@@ -222,7 +231,7 @@ void ActionWithMatrix::gatherForcesOnStoredValue( const Value* myval, const unsi
 }
 
 void ActionWithMatrix::clearMatrixElements( MultiValue& myvals ) const {
-  if( isActive() ) {
+  if( isActive() && clearOnEachCycle ) {
     for(int i=0; i<getNumberOfComponents(); ++i) {
       const Value* myval=getConstPntrToComponent(i);
       if( myval->getRank()==2 && !myval->hasDerivatives() ) myvals.clearDerivatives( myval->getPositionInStream() );

src/core/ActionWithMatrix.h

@@ -44,6 +44,8 @@ class ActionWithMatrix : public ActionWithVector {
 /// This does the calculation of a particular matrix element
   void runTask( const std::string& controller, const unsigned& current, const unsigned colno, MultiValue& myvals ) const ;
 protected:
+/// This turns off derivative clearing for contact matrix if we are not storing derivatives
+  bool clearOnEachCycle;
 /// Does the matrix chain continue on from this action
   bool matrixChainContinues() const ;
 /// This returns the jelem th element of argument ic
@@ -132,14 +134,8 @@ inline
 void ActionWithMatrix::addDerivativeOnMatrixArgument( const bool& inchain, const unsigned& ival, const unsigned& jarg, const unsigned& irow, const unsigned& jcol, const double& der, MultiValue& myvals ) const {
   plumed_dbg_assert( jarg<getNumberOfArguments() && getPntrToArgument(jarg)->getRank()==2 && !getPntrToArgument(jarg)->hasDerivatives() );
   unsigned ostrn = getConstPntrToComponent(ival)->getPositionInStream(), vstart=arg_deriv_starts[jarg];
-  if( !inchain && getPntrToArgument(jarg)->getNumberOfColumns()<getPntrToArgument(jarg)->getShape()[1] ) {
-    unsigned dloc = vstart + irow*getPntrToArgument(jarg)->getNumberOfColumns(); Value* myarg=getPntrToArgument(jarg);
-    for(unsigned i=0; i<myarg->getRowLength(irow); ++i) {
-      if( myarg->getRowIndex(irow,i)==jcol ) { myvals.addDerivative( ostrn, dloc+i, der ); myvals.updateIndex( ostrn, dloc+i ); return; }
-    }
-    plumed_merror("could not find element of sparse matrix to add derivative to");
-  } else if( !inchain ) {
-    unsigned dloc = vstart + irow*getPntrToArgument(jarg)->getShape()[1] + jcol;
+  if( !inchain ) {
+    unsigned dloc = vstart + irow*getPntrToArgument(jarg)->getNumberOfColumns() + jcol;
     myvals.addDerivative( ostrn, dloc, der ); myvals.updateIndex( ostrn, dloc );
   } else {
     unsigned istrn = getPntrToArgument(jarg)->getPositionInStream();

src/matrixtools/MatrixTimesVector.cpp

@@ -152,32 +152,44 @@ void MatrixTimesVector::prepare() {
 void MatrixTimesVector::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
   unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(0)->getRowLength(task_index);
   if( indices.size()!=size_v+1 ) indices.resize( size_v + 1 );
-  for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + getPntrToArgument(0)->getRowIndex( task_index, i );
+  for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + i;
   myvals.setSplitIndex( size_v + 1 );
 }
 
 void MatrixTimesVector::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
   unsigned ind2 = index2; if( index2>=getPntrToArgument(0)->getShape()[0] ) ind2 = index2 - getPntrToArgument(0)->getShape()[0];
   if( getPntrToArgument(1)->getRank()==1 ) {
+    double matval = 0; Value* myarg = getPntrToArgument(0); unsigned vcol = ind2;
+    if( !myarg->valueHasBeenSet() ) matval = myvals.get( myarg->getPositionInStream() );
+    else {
+      matval = myarg->get( index1*myarg->getNumberOfColumns() + ind2, false );
+      vcol = getPntrToArgument(0)->getRowIndex( index1, ind2 );
+    }
     for(unsigned i=0; i<getNumberOfArguments()-1; ++i) {
       unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream();
-      double matval = getElementOfMatrixArgument( 0, index1, ind2, myvals ), vecval=getArgumentElement( i+1, ind2, myvals );
+      double vecval=getArgumentElement( i+1, vcol, myvals );
       // And add this part of the product
       myvals.addValue( ostrn, matval*vecval );
       // Now lets work out the derivatives
       if( doNotCalculateDerivatives() ) continue;
-      addDerivativeOnMatrixArgument( stored_arg[0], i, 0, index1, ind2, vecval, myvals ); addDerivativeOnVectorArgument( stored_arg[i+1], i, i+1, ind2, matval, myvals );
+      addDerivativeOnMatrixArgument( stored_arg[0], i, 0, index1, ind2, vecval, myvals ); addDerivativeOnVectorArgument( stored_arg[i+1], i, i+1, vcol, matval, myvals );
     }
   } else {
-    unsigned n=getNumberOfArguments()-1;
+    unsigned n=getNumberOfArguments()-1; double matval = 0; unsigned vcol = ind2;
     for(unsigned i=0; i<getNumberOfArguments()-1; ++i) {
       unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream();
-      double matval = getElementOfMatrixArgument( i, index1, ind2, myvals ), vecval=getArgumentElement( n, ind2, myvals );
+      Value* myarg = getPntrToArgument(i);
+      if( !myarg->valueHasBeenSet() ) matval = myvals.get( myarg->getPositionInStream() );
+      else {
+        matval = myarg->get( index1*myarg->getNumberOfColumns() + ind2, false );
+        vcol = getPntrToArgument(i)->getRowIndex( index1, ind2 );
+      }
+      double vecval=getArgumentElement( n, vcol, myvals );
       // And add this part of the product
       myvals.addValue( ostrn, matval*vecval );
       // Now lets work out the derivatives
       if( doNotCalculateDerivatives() ) continue;
-      addDerivativeOnMatrixArgument( stored_arg[i], i, i, index1, ind2, vecval, myvals ); addDerivativeOnVectorArgument( stored_arg[n], i, n, ind2, matval, myvals );
+      addDerivativeOnMatrixArgument( stored_arg[i], i, i, index1, ind2, vecval, myvals ); addDerivativeOnVectorArgument( stored_arg[n], i, n, vcol, matval, myvals );
     }
   }
 }