Some optimization improvements in histograms

src/analysis/Histogram.cpp

@@ -311,6 +311,8 @@ void Histogram::turnOnDerivatives() {
     if( !mbase ) error("do not know how to get histogram derivatives for actions of type " + myvessels[i]->getName() );
     tmp_atoms = mbase->getAbsoluteIndexes();
     for(unsigned j=0; j<tmp_atoms.size(); ++j) all_atoms.push_back( tmp_atoms[j] );
+    // Make a tempory multi value so we can avoid vector resizing
+    stashes[i]->resizeTemporyMultiValues( 1 );
   }
   ActionAtomistic::requestAtoms( all_atoms );
   finalForces.resize( 3*all_atoms.size() + 9 );
@@ -386,7 +388,10 @@ void Histogram::compute( const unsigned& current, MultiValue& myvals ) const {
     for(unsigned i=2; i<myvessels[0]->getNumberOfQuantities(); ++i) myvals.setValue( i-1, cvals[i] );
     myvals.setValue( 0, cvals[0] ); myvals.setValue( myvessels[0]->getNumberOfQuantities() - 1, ww );
     if( in_apply ) {
-      MultiValue tmpval( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
+      MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
+      if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
+          tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
+        tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
       stashes[0]->retrieveDerivatives( stashes[0]->getTrueIndex(current), true, tmpval );
       for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
         unsigned jder=tmpval.getActiveIndex(j); myvals.addDerivative( 0, jder, tmpval.getDerivative(0, jder) );
@@ -406,7 +411,10 @@ void Histogram::compute( const unsigned& current, MultiValue& myvals ) const {
         stashes[j]->retrieveSequentialValue( current, false, cvals ); totweight *= cvals[0];
       }
       // And this bit the derivatives
-      MultiValue tmpval( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
+      MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
+      if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
+          tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
+        tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
       stashes[0]->retrieveDerivatives( stashes[0]->getTrueIndex(current), false, tmpval );
       for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
         unsigned jder=tmpval.getActiveIndex(j);
@@ -421,7 +429,10 @@ void Histogram::compute( const unsigned& current, MultiValue& myvals ) const {
       tnorm *= cvals[0]; myvals.setValue( 1+i, cvals[1] );
       // Get the derivatives as well if we are in apply
       if( in_apply ) {
-        MultiValue tmpval( myvessels[i]->getNumberOfQuantities(), myvessels[i]->getNumberOfDerivatives() );
+        MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
+        if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
+            tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
+          tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
         stashes[i]->retrieveDerivatives( stashes[i]->getTrueIndex(current), false, tmpval );
         for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
           unsigned jder=tmpval.getActiveIndex(j);

src/gridtools/GridVessel.cpp

@@ -191,9 +191,13 @@ void GridVessel::getIndices( const unsigned& index, std::vector<unsigned>& indic
 }
 
 void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
-  plumed_dbg_assert( bounds_set && x.size()==dimension && ipoint<npoints );
+  std::vector<unsigned> tindices( dimension ); getGridPointCoordinates( ipoint, tindices, x );
+}
+
+void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
+  plumed_dbg_assert( bounds_set && x.size()==dimension && tindices.size()==dimension && ipoint<npoints );
   if( gtype==flat ) {
-    std::vector<unsigned> tindices( dimension ); getIndices( ipoint, tindices );
+    getIndices( ipoint, tindices );
     for(unsigned i=0; i<dimension; ++i) x[i] = min[i] + dx[i]*tindices[i];
   } else if( gtype==fibonacci ) {
     x[1] = ((ipoint*fib_offset) - 1) + (fib_offset/2);

src/gridtools/GridVessel.h

@@ -117,6 +117,7 @@ public:
   unsigned getNumberOfPoints() const;
 /// Get the coordinates for a point in the grid
   void getGridPointCoordinates( const unsigned&, std::vector<double>& ) const ;
+  void getGridPointCoordinates( const unsigned&, std::vector<unsigned>&, std::vector<double>& ) const ;
 /// Get the dimensionality of the function
   unsigned getDimension() const ;
 /// Get the number of components in the vector stored on each grid point

src/gridtools/HistogramOnGrid.cpp

@@ -121,11 +121,11 @@ void HistogramOnGrid::calculate( const unsigned& current, MultiValue& myvals, st
       std::vector<double> intforce( 2*dimension, 0.0 );
       std::vector<Value*> vv( getVectorOfValues() );
 
-      double newval; std::vector<double> xx( dimension );
+      double newval; std::vector<unsigned> tindices( dimension ); std::vector<double> xx( dimension );
       for(unsigned i=0; i<num_neigh; ++i) {
         unsigned ineigh=neighbors[i];
         if( inactive( ineigh ) ) continue ;
-        getGridPointCoordinates( ineigh, xx );
+        getGridPointCoordinates( ineigh, tindices, xx );
         if( kernel ) {
           for(unsigned j=0; j<dimension; ++j) vv[j]->set(xx[j]);
           newval = kernel->evaluate( vv, der, true );