diff --git a/src/adjmat/AdjacencyMatrixBase.cpp b/src/adjmat/AdjacencyMatrixBase.cpp index e658218742..f04cdc294c 100644 --- a/src/adjmat/AdjacencyMatrixBase.cpp +++ b/src/adjmat/AdjacencyMatrixBase.cpp @@ -206,25 +206,6 @@ void AdjacencyMatrixBase::updateNeighbourList() { if( read_one_group && maxcol==getConstPntrToComponent(0)->getShape()[1] ) maxcol = maxcol-1; } -void AdjacencyMatrixBase::getAdditionalTasksRequired( ActionWithVector* action, std::vector& atasks ) { - if( action==this ) return; - // Update the neighbour list - updateNeighbourList(); - - unsigned nactive = atasks.size(); - std::vector indlist( 1 + ablocks.size() + threeblocks.size() ); - for(unsigned i=0; i & indices ) const { unsigned natoms=nlist[current]; indices[0]=current; unsigned lstart = getConstPntrToComponent(0)->getShape()[0] + current*(1+natoms_per_list); plumed_dbg_assert( nlist[lstart]==current ); @@ -278,10 +259,9 @@ void AdjacencyMatrixBase::performTask( const std::string& controller, const unsi // Update dynamic list indices for virial unsigned base = 3*getNumberOfAtoms(); for(unsigned j=0; j<9; ++j) myvals.updateIndex( 0, base+j ); // And the indices for the derivatives of the row of the matrix - unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat ); - std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); + unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives(); std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices() ); matrix_indices[nmat_ind+0]=3*index2+0; matrix_indices[nmat_ind+1]=3*index2+1; matrix_indices[nmat_ind+2]=3*index2+2; - myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+3 ); + myvals.setNumberOfMatrixRowDerivatives( nmat_ind+3 ); } // Calculate the components if we need them @@ -319,12 +299,6 @@ void AdjacencyMatrixBase::performTask( const std::string& controller, const unsi myvals.addDerivative( 3, base+1, 0 ); myvals.addDerivative( 3, base+4, 0 ); myvals.addDerivative( 3, base+7, 0 ); myvals.addDerivative( 3, base+2, -atom[0] ); myvals.addDerivative( 3, base+5, -atom[1] ); myvals.addDerivative( 3, base+8, -atom[2] ); for(unsigned k=0; k<9; ++k) { myvals.updateIndex( 1, base+k ); myvals.updateIndex( 2, base+k ); myvals.updateIndex( 3, base+k ); } - for(unsigned k=1; k<4; ++k) { - unsigned nmat = getConstPntrToComponent(k)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat ); - std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); - matrix_indices[nmat_ind+0]=3*index2+0; matrix_indices[nmat_ind+1]=3*index2+1; matrix_indices[nmat_ind+2]=3*index2+2; - myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+3 ); - } } } } @@ -332,25 +306,23 @@ void AdjacencyMatrixBase::performTask( const std::string& controller, const unsi void AdjacencyMatrixBase::runEndOfRowJobs( const unsigned& ind, const std::vector & indices, MultiValue& myvals ) const { if( doNotCalculateDerivatives() ) return; - for(int k=0; kgetPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat ); - std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); - plumed_assert( nmat_ind& matrix_indices( myvals.getMatrixRowDerivativeIndices() ); + plumed_assert( nmat_ind& atasks ) override ; void setupForTask( const unsigned& current, std::vector & indices, MultiValue& myvals ) const override; // void setupCurrentTaskList() override; void updateNeighbourList() override ; diff --git a/src/adjmat/Neighbors.cpp b/src/adjmat/Neighbors.cpp index 709e1714d1..4dab5c918a 100644 --- a/src/adjmat/Neighbors.cpp +++ b/src/adjmat/Neighbors.cpp @@ -42,7 +42,6 @@ class Neighbors : public ActionWithMatrix { explicit Neighbors(const ActionOptions&); unsigned getNumberOfDerivatives() override; unsigned getNumberOfColumns() const override { return number; } - bool canBeAfterInChain( ActionWithVector* av ) { return av->getLabel()!=(getPntrToArgument(0)->getPntrToAction())->getLabel(); } void setupForTask( const unsigned& task_index, std::vector& indices, MultiValue& myvals ) const ; void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override; void runEndOfRowJobs( const unsigned& ival, const std::vector & indices, MultiValue& myvals ) const override {} @@ -64,7 +63,6 @@ Neighbors::Neighbors(const ActionOptions&ao): { if( getNumberOfArguments()!=1 ) error("found wrong number of arguments in input"); if( getPntrToArgument(0)->getRank()!=2 ) error("input argument should be a matrix"); - getPntrToArgument(0)->buildDataStore(); unsigned nlow; parse("NLOWEST",nlow); unsigned nhigh; parse("NHIGHEST",nhigh); diff --git a/src/adjmat/TorsionsMatrix.cpp b/src/adjmat/TorsionsMatrix.cpp index 563c6e7cfa..87737fec44 100644 --- a/src/adjmat/TorsionsMatrix.cpp +++ b/src/adjmat/TorsionsMatrix.cpp @@ -36,9 +36,6 @@ namespace PLMD { namespace adjmat { class TorsionsMatrix : public ActionWithMatrix { -private: - unsigned nderivatives; - bool stored_matrix1, stored_matrix2; public: static void registerKeywords( Keywords& keys ); explicit TorsionsMatrix(const ActionOptions&); @@ -85,10 +82,7 @@ TorsionsMatrix::TorsionsMatrix(const ActionOptions&ao): log.printf("\n"); requestAtoms( atoms_a, false ); std::vector shape(2); shape[0]=getPntrToArgument(0)->getShape()[0]; shape[1]=getPntrToArgument(1)->getShape()[1]; - addValue( shape ); setPeriodic("-pi","pi"); nderivatives = buildArgumentStore(0) + 3*getNumberOfAtoms() + 9; - std::string headstr=getFirstActionInChain()->getLabel(); - stored_matrix1 = getPntrToArgument(0)->ignoreStoredValue( headstr ); - stored_matrix2 = getPntrToArgument(1)->ignoreStoredValue( headstr ); + addValue( shape ); setPeriodic("-pi","pi"); if( nm>0 ) { unsigned iarg = getNumberOfArguments()-1; @@ -98,7 +92,7 @@ TorsionsMatrix::TorsionsMatrix(const ActionOptions&ao): } unsigned TorsionsMatrix::getNumberOfDerivatives() { - return nderivatives; + return 3*getNumberOfAtoms() + 9 + getPntrToArgument(0)->getNumberOfStoredValues() + getPntrToArgument(1)->getNumberOfStoredValues(); } unsigned TorsionsMatrix::getNumberOfColumns() const { @@ -131,8 +125,8 @@ void TorsionsMatrix::performTask( const std::string& controller, const unsigned& // Get the two vectors for(unsigned i=0; i<3; ++i) { - v1[i] = getElementOfMatrixArgument( 0, index1, i, myvals ); - v2[i] = getElementOfMatrixArgument( 1, i, ind2, myvals ); + v1[i] = getPntrToArgument(0)->get( index1*getPntrToArgument(0)->getNumberOfColumns() + i, false ); + v2[i] = getPntrToArgument(1)->get( i*getPntrToArgument(1)->getNumberOfColumns() + ind2, false ); } // Evaluate angle Torsion t; double angle = t.compute( v1, conn, v2, dv1, dconn, dv2 ); @@ -141,9 +135,12 @@ void TorsionsMatrix::performTask( const std::string& controller, const unsigned& if( doNotCalculateDerivatives() ) return; // Add the derivatives on the matrices + unsigned base1 = index1*getPntrToArgument(0)->getNumberOfColumns(); + unsigned ncols = getPntrToArgument(1)->getNumberOfColumns(); + unsigned base2 = getPntrToArgument(0)->getNumberOfStoredValues() + ind2; for(unsigned i=0; i<3; ++i) { - addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, i, dv1[i], myvals ); - addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, i, ind2, dv2[i], myvals ); + myvals.addDerivative( 0, base1 + i, dv1[i] ); myvals.updateIndex( 0, base1 + i ); + myvals.addDerivative( 0, base2 + i*ncols, dv2[i] ); myvals.updateIndex( 0, base2 + i*ncols ); } // And derivatives on positions unsigned narg_derivatives = getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues(); @@ -160,20 +157,20 @@ void TorsionsMatrix::runEndOfRowJobs( const unsigned& ival, const std::vectorgetShape()[1]; - unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat ); + unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives(); unsigned narg_derivatives = getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues(); - std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); unsigned ntwo_atoms = myvals.getSplitIndex(); + std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices() ); unsigned ntwo_atoms = myvals.getSplitIndex(); for(unsigned j=0; j<3; ++j) { matrix_indices[nmat_ind] = mat1s + j; nmat_ind++; matrix_indices[nmat_ind] = narg_derivatives + mat1s + j; nmat_ind++; for(unsigned i=1; i=getPntrToArgument(0)->getShape()[0] ) ind2 = indices[i] - getPntrToArgument(0)->getShape()[0]; - matrix_indices[nmat_ind] = arg_deriv_starts[1] + j*ss + ind2; nmat_ind++; + matrix_indices[nmat_ind] = getPntrToArgument(0)->getNumberOfStoredValues() + j*ss + ind2; nmat_ind++; matrix_indices[nmat_ind] = narg_derivatives + 3*indices[i] + j; nmat_ind++; } } unsigned base = narg_derivatives + 3*getNumberOfAtoms(); for(unsigned j=0; j<9; ++j) { matrix_indices[nmat_ind] = base + j; nmat_ind++; } - myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind ); + myvals.setNumberOfMatrixRowDerivatives( nmat_ind ); } } diff --git a/src/clusters/ClusterDistribution.cpp b/src/clusters/ClusterDistribution.cpp index 7a38671a1c..7480fabffa 100644 --- a/src/clusters/ClusterDistribution.cpp +++ b/src/clusters/ClusterDistribution.cpp @@ -127,11 +127,9 @@ ClusterDistribution::ClusterDistribution(const ActionOptions&ao): } // Request the arguments requestArguments( clusters ); - getPntrToArgument(0)->buildDataStore(); - if( getNumberOfArguments()>1 ) getPntrToArgument(1)->buildDataStore(); // Now create the value std::vector shape(1); shape[0]=clusters[0]->getShape()[0]; - addValue( shape ); setNotPeriodic(); getPntrToValue()->buildDataStore(); + addValue( shape ); setNotPeriodic(); } unsigned ClusterDistribution::getNumberOfDerivatives() { diff --git a/src/clusters/ClusterWeights.cpp b/src/clusters/ClusterWeights.cpp index e19bcd0bcb..347427d0c6 100644 --- a/src/clusters/ClusterWeights.cpp +++ b/src/clusters/ClusterWeights.cpp @@ -89,7 +89,7 @@ ClusterWeights::ClusterWeights(const ActionOptions&ao): requestArguments( clusters ); // Now create the value std::vector shape(1); shape[0]=clusters[0]->getShape()[0]; - addValue( shape ); setNotPeriodic(); getPntrToComponent(0)->buildDataStore(); + addValue( shape ); setNotPeriodic(); // Find out which cluster we want parse("CLUSTER",clustr); if( clustr<1 ) error("cannot look for a cluster larger than the largest cluster"); diff --git a/src/clusters/ClusteringBase.cpp b/src/clusters/ClusteringBase.cpp index 357577a3ae..ca54fa6da9 100644 --- a/src/clusters/ClusteringBase.cpp +++ b/src/clusters/ClusteringBase.cpp @@ -42,7 +42,7 @@ ClusteringBase::ClusteringBase(const ActionOptions&ao): // Now create a value - this holds the data on which cluster each guy is in std::vector shape(1); shape[0]=getPntrToArgument(0)->getShape()[0]; // Build the store here to make sure that next action has all data - addValue( shape ); setNotPeriodic(); getPntrToValue()->buildDataStore(); + addValue( shape ); setNotPeriodic(); // Resize local variables which_cluster.resize( getPntrToArgument(0)->getShape()[0] ); cluster_sizes.resize( getPntrToArgument(0)->getShape()[0] ); } diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h index 0811736aa4..bd11210de5 100644 --- a/src/colvar/MultiColvarTemplate.h +++ b/src/colvar/MultiColvarTemplate.h @@ -44,7 +44,6 @@ class MultiColvarTemplate : public ActionWithVector { unsigned getNumberOfDerivatives() override ; void addValueWithDerivatives( const std::vector& shape=std::vector() ) override ; void addComponentWithDerivatives( const std::string& name, const std::vector& shape=std::vector() ) override ; - void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) override ; void performTask( const unsigned&, MultiValue& ) const override ; void calculate() override; }; @@ -114,6 +113,7 @@ unsigned MultiColvarTemplate::getNumberOfDerivatives() { template void MultiColvarTemplate::calculate() { + if( wholemolecules ) makeWhole(); runAllTasks(); } @@ -127,12 +127,6 @@ void MultiColvarTemplate::addComponentWithDerivatives( const std::string& nam std::vector s(1); s[0]=ablocks[0].size(); addComponent( name, s ); } -template -void MultiColvarTemplate::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) { - if( wholemolecules ) makeWhole(); - ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol ); -} - template void MultiColvarTemplate::performTask( const unsigned& task_index, MultiValue& myvals ) const { // Retrieve the positions diff --git a/src/colvar/RMSDVector.cpp b/src/colvar/RMSDVector.cpp index 021fa2574f..4c9c22152c 100644 --- a/src/colvar/RMSDVector.cpp +++ b/src/colvar/RMSDVector.cpp @@ -94,13 +94,13 @@ RMSDVector::RMSDVector(const ActionOptions&ao): if( displacement && (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]<=1) ) { addComponentWithDerivatives("dist"); componentIsNotPeriodic("dist"); std::vector shape( 1, getPntrToArgument(0)->getNumberOfValues() ); - addComponent( "disp", shape ); getPntrToComponent(1)->buildDataStore(); componentIsNotPeriodic("disp"); + addComponent( "disp", shape ); componentIsNotPeriodic("disp"); } else if( displacement ) { std::vector shape( 1, getPntrToArgument(1)->getShape()[0] ); - addComponent( "dist", shape ); getPntrToComponent(0)->buildDataStore(); + addComponent( "dist", shape ); componentIsNotPeriodic("dist"); shape.resize(2); shape[0] = getPntrToArgument(1)->getShape()[0]; shape[1] = getPntrToArgument(0)->getNumberOfValues(); - addComponent( "disp", shape ); getPntrToComponent(1)->buildDataStore(); getPntrToComponent(1)->reshapeMatrixStore( shape[1] ); + addComponent( "disp", shape ); getPntrToComponent(1)->reshapeMatrixStore( shape[1] ); componentIsNotPeriodic("disp"); } else if( (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]==1) ) { addValue(); setNotPeriodic(); @@ -278,11 +278,18 @@ void RMSDVector::gatherStoredValue( const unsigned& valindex, const unsigned& co } } -void RMSDVector::gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const { - if( getConstPntrToComponent(ival)->getRank()==1 ) { ActionWithVector::gatherForcesOnStoredValue( ival, itask, myvals, forces ); return; } - const std::vector& direction( myvals.getConstFirstAtomDerivativeVector()[1] ); unsigned natoms = direction.size(); - for(unsigned i=0; i& forces ) const { + if( checkComponentsForForce() ) { + for(unsigned ival=0; ivalgetRank()==1 ) { + gatherForcesOnStoredValue( ival, itask, myvals, forces ); + } else { + const std::vector& direction( myvals.getConstFirstAtomDerivativeVector()[1] ); unsigned natoms = direction.size(); + for(unsigned i=0; i& buffer ) const override ; - void gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const override ; + void gatherForces( const unsigned& i, const MultiValue& myvals, std::vector& forces ) const override ; void setReferenceConfigurations(); void calculate() override ; bool checkForTaskForce( const unsigned& itask, const Value* myval ) const override ; diff --git a/src/contour/DistanceFromContourBase.cpp b/src/contour/DistanceFromContourBase.cpp index 6444830865..1334a18fca 100644 --- a/src/contour/DistanceFromContourBase.cpp +++ b/src/contour/DistanceFromContourBase.cpp @@ -46,10 +46,8 @@ DistanceFromContourBase::DistanceFromContourBase( const ActionOptions& ao ): nactive(0) { if( getNumberOfArguments()>1 ) error("should only use one argument for this action"); - if( getNumberOfArguments()==1 ) { - if( getPntrToArgument(0)->getRank()!=1 ) error("ARG for distance from contour should be rank one"); - getPntrToArgument(0)->buildDataStore(); - } + if( getNumberOfArguments()==1 && getPntrToArgument(0)->getRank()!=1 ) error("ARG for distance from contour should be rank one"); + // Read in the multicolvar/atoms std::vector atoms; parseAtomList("POSITIONS",atoms); std::vector origin; parseAtomList("ATOM",origin); diff --git a/src/contour/FindContour.cpp b/src/contour/FindContour.cpp index 0ab14ac7b1..4ca0508955 100644 --- a/src/contour/FindContour.cpp +++ b/src/contour/FindContour.cpp @@ -109,11 +109,7 @@ FindContour::FindContour(const ActionOptions&ao): std::vector argn( ag->getGridCoordinateNames() ); std::vector shape(1); shape[0]=0; - for(unsigned i=0; ibuildDataStore(); - } - // Check for task reduction - updateTaskListReductionStatus(); + for(unsigned i=0; i shape( getInputGridObject().getDimension()-1 ); addValueWithDerivatives( shape ); setNotPeriodic(); - getPntrToComponent(0)->buildDataStore(); } void FindContourSurface::setupValuesOnFirstStep() { diff --git a/src/contour/FindSphericalContour.cpp b/src/contour/FindSphericalContour.cpp index 5479c365c3..d606c464f8 100644 --- a/src/contour/FindSphericalContour.cpp +++ b/src/contour/FindSphericalContour.cpp @@ -151,7 +151,6 @@ FindSphericalContour::FindSphericalContour(const ActionOptions&ao): // Now create a value std::vector shape( 3 ); shape[0]=npoints; shape[1]=shape[2]=1; addValueWithDerivatives( shape ); setNotPeriodic(); - getPntrToComponent(0)->buildDataStore(); checkRead(); } unsigned FindSphericalContour::getNumberOfDerivatives() { diff --git a/src/core/ActionToGetData.cpp b/src/core/ActionToGetData.cpp index 8c363b1e61..a7047fb588 100644 --- a/src/core/ActionToGetData.cpp +++ b/src/core/ActionToGetData.cpp @@ -60,7 +60,7 @@ ActionToGetData::ActionToGetData(const ActionOptions&ao): if( getNumberOfArguments()!=1 ) error("python interface works best when you ask for one argument at a time"); if( getPntrToArgument(0)->getNumberOfValues()==0 ) error("cannot get data as shape of value " + getPntrToArgument(0)->getName() + " has not been set"); - getPntrToArgument(0)->buildDataStore(); data.resize( getPntrToArgument(0)->getNumberOfValues() ); + data.resize( getPntrToArgument(0)->getNumberOfValues() ); } void ActionToGetData::get_rank( const TypesafePtr & dims ) { diff --git a/src/core/ActionWithArguments.cpp b/src/core/ActionWithArguments.cpp index 52ef2d1e21..8a271d4d35 100644 --- a/src/core/ActionWithArguments.cpp +++ b/src/core/ActionWithArguments.cpp @@ -306,10 +306,8 @@ void ActionWithArguments::addForcesOnArguments( const unsigned& argstart, const if( av && av->getNumberOfMasks()>0 ) nargs=nargs-av->getNumberOfMasks(); for(unsigned i=0; istoredata || arguments[i]->getRank()==0 || (arguments[i]->getRank()>0 && arguments[i]->hasDerivatives()) ) { - unsigned nvals = arguments[i]->getNumberOfStoredValues(); - for(unsigned j=0; jaddForce( j, forces[ind], false ); ind++; } - } + unsigned nvals = arguments[i]->getNumberOfStoredValues(); + for(unsigned j=0; jaddForce( j, forces[ind], false ); ind++; } } } diff --git a/src/core/ActionWithMatrix.cpp b/src/core/ActionWithMatrix.cpp index a7521d0a2c..4b4443ee6f 100644 --- a/src/core/ActionWithMatrix.cpp +++ b/src/core/ActionWithMatrix.cpp @@ -30,110 +30,33 @@ void ActionWithMatrix::registerKeywords( Keywords& keys ) { ActionWithMatrix::ActionWithMatrix(const ActionOptions&ao): Action(ao), - ActionWithVector(ao), - next_action_in_chain(NULL), - matrix_to_do_before(NULL), - matrix_to_do_after(NULL), - clearOnEachCycle(true) + ActionWithVector(ao) { } -ActionWithMatrix::~ActionWithMatrix() { - if( matrix_to_do_before ) { matrix_to_do_before->matrix_to_do_after=NULL; matrix_to_do_before->next_action_in_chain=NULL; } -} - -void ActionWithMatrix::getAllActionLabelsInMatrixChain( std::vector& mylabels ) const { - bool found=false; - for(unsigned i=0; igetAllActionLabelsInMatrixChain( mylabels ); -} - -void ActionWithMatrix::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) { - ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol ); - +void ActionWithMatrix::calculate() { + // Update all the neighbour lists + updateNeighbourList(); + // Setup the matrix indices for(int i=0; igetRank()!=2 || myval->hasDerivatives() ) continue; - myval->setPositionInMatrixStash(nmat); nmat++; - if( !myval->valueIsStored() ) continue; - if( myval->getShape()[1]>maxcol ) maxcol=myval->getShape()[1]; - } - // 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; igetRank()==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 ) { - ActionWithMatrix* am=dynamic_cast(act); if( !am || act==this ) return; - // Build the list that contains everything we are going to loop over in getTotalMatrixBookeepgin and updateAllNeighbourLists - if( next_action_in_chain ) next_action_in_chain->finishChainBuild( act ); - else { - next_action_in_chain=am; - // Build the list of things we are going to loop over in runTask - if( am->isAdjacencyMatrix() || act->getName()=="VSTACK" ) return ; - plumed_massert( !matrix_to_do_after, "cannot add " + act->getLabel() + " in " + getLabel() + " as have already added " + matrix_to_do_after->getLabel() ); - matrix_to_do_after=am; am->matrix_to_do_before=this; - } -} - -const ActionWithMatrix* ActionWithMatrix::getFirstMatrixInChain() const { - if( !actionInChain() ) return this; - return matrix_to_do_before->getFirstMatrixInChain(); -} - -void ActionWithMatrix::setupMatrixStore() { - for(int i=0; igetRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() ) continue; myval->reshapeMatrixStore( getNumberOfColumns() ); } - if( next_action_in_chain ) next_action_in_chain->setupMatrixStore(); -} - -void ActionWithMatrix::calculate() { - if( actionInChain() ) return ; - // Update all the neighbour lists - updateAllNeighbourLists(); - // Setup the matrix indices - setupMatrixStore(); // And run all the tasks runAllTasks(); } -void ActionWithMatrix::updateAllNeighbourLists() { - updateNeighbourList(); - if( next_action_in_chain ) next_action_in_chain->updateAllNeighbourLists(); -} +void ActionWithMatrix::performTask( const unsigned& task_index, MultiValue& myvals ) const { + std::vector & indices( myvals.getIndices() ); + setupForTask( task_index, indices, myvals ); -void ActionWithMatrix::clearBookeepingBeforeTask( const unsigned& task_index ) const { // Reset the bookeeping elements for storage for(unsigned i=0; i( getConstPntrToComponent(i) ); unsigned ncols = myval->getNumberOfColumns(); - if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() || ncols>=myval->getShape()[1] ) continue; + if( myval->getRank()!=2 || myval->hasDerivatives() || ncols>=myval->getShape()[1] ) continue; myval->matrix_bookeeping[task_index*(1+ncols)]=0; } - if( matrix_to_do_after ) matrix_to_do_after->clearBookeepingBeforeTask( task_index ); -} - -void ActionWithMatrix::performTask( const unsigned& task_index, MultiValue& myvals ) const { - std::vector & indices( myvals.getIndices() ); - if( matrix_to_do_before ) { - plumed_dbg_assert( myvals.inVectorCall() ); - runEndOfRowJobs( task_index, indices, myvals ); - return; - } - setupForTask( task_index, indices, myvals ); - - // Reset the bookeeping elements for storage - clearBookeepingBeforeTask( task_index ); // Now loop over the row of the matrix unsigned ntwo_atoms = myvals.getSplitIndex(); @@ -159,15 +82,15 @@ void ActionWithMatrix::runTask( const std::string& controller, const unsigned& c double checkval = myvals.get( 0 ); for(int i=0; igetNumberOfColumns(); - if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() ) continue; - unsigned matindex = myval->getPositionInMatrixStash(), col_stash_index = colno; + if( myval->getRank()!=2 || myval->hasDerivatives() ) continue; + unsigned col_stash_index = colno; if( colno>=myval->getShape()[0] ) col_stash_index = colno - myval->getShape()[0]; if( myval->forcesWereAdded() ) { double fforce; if( ncolsgetShape()[1] ) fforce = myval->getForce( myvals.getTaskIndex()*ncols + myval->matrix_bookeeping[current*(1+ncols)] ); else fforce = myval->getForce( myvals.getTaskIndex()*myval->getShape()[1] + col_stash_index ); for(unsigned j=0; jrunTask( controller, current, colno, myvals ); } bool ActionWithMatrix::checkForTaskForce( const unsigned& itask, const Value* myval ) const { @@ -192,20 +114,18 @@ bool ActionWithMatrix::checkForTaskForce( const unsigned& itask, const Value* my return false; } -void ActionWithMatrix::gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const { - if( getConstPntrToComponent(ival)->getRank()==1 ) { ActionWithVector::gatherForcesOnStoredValue( ival, itask, myvals, forces ); return; } - unsigned matind = getConstPntrToComponent(ival)->getPositionInMatrixStash(); const std::vector& mat_indices( myvals.getMatrixRowDerivativeIndices( matind ) ); - for(unsigned i=0; i& forces ) const { + if( checkComponentsForForce() ) { + const std::vector& mat_indices( myvals.getMatrixRowDerivativeIndices() ); + for(unsigned i=0; igetRank()==2 && !myval->hasDerivatives() ) myvals.clearDerivatives( i ); - } + for(int i=0; igetRank()==2 && !myval->hasDerivatives() ) myvals.clearDerivatives( i ); } - if( matrix_to_do_after ) matrix_to_do_after->clearMatrixElements( myvals ); } } diff --git a/src/core/ActionWithMatrix.h b/src/core/ActionWithMatrix.h index 3428393a43..7f90d3b7a5 100644 --- a/src/core/ActionWithMatrix.h +++ b/src/core/ActionWithMatrix.h @@ -28,48 +28,17 @@ namespace PLMD { class ActionWithMatrix : public ActionWithVector { private: - ActionWithMatrix* next_action_in_chain; - ActionWithMatrix* matrix_to_do_before; - ActionWithMatrix* matrix_to_do_after; -/// Update all the neighbour lists in the chain - void updateAllNeighbourLists(); -/// This clears all bookeeping arrays before the ith task - void clearBookeepingBeforeTask( const unsigned& task_index ) const ; /// This is used to clear up the matrix elements void clearMatrixElements( MultiValue& myvals ) const ; -/// This is used to find the total amount of space we need for storing matrix elements - void setupMatrixStore(); /// 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 - double getArgumentElement( const unsigned& ic, const unsigned& jelem, const MultiValue& myvals ) const ; -/// This returns an element of a matrix that is passed an argument - double getElementOfMatrixArgument( const unsigned& imat, const unsigned& irow, const unsigned& jcol, const MultiValue& myvals ) const ; -/// Add derivatives given the derivative wrt to the input vector element as input - void addDerivativeOnVectorArgument( const bool& inchain, const unsigned& ival, const unsigned& jarg, const unsigned& jelem, const double& der, MultiValue& myvals ) const ; -/// Add derivatives given the derative wrt to the input matrix element as input - void addDerivativeOnMatrixArgument( const bool& inchain, const unsigned& ival, const unsigned& jarg, const unsigned& irow, const unsigned& jcol, const double& der, MultiValue& myvals ) const ; public: static void registerKeywords( Keywords& keys ); explicit ActionWithMatrix(const ActionOptions&); - virtual ~ActionWithMatrix(); /// virtual bool isAdjacencyMatrix() const { return false; } -/// - void getAllActionLabelsInMatrixChain( std::vector& mylabels ) const override ; -/// Get the first matrix in this chain - const ActionWithMatrix* getFirstMatrixInChain() const ; -/// - void finishChainBuild( ActionWithVector* act ); /// This should return the number of columns to help with sparse storage of matrices virtual unsigned getNumberOfColumns() const = 0; -/// This requires some thought - void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) override; //// This does some setup before we run over the row of the matrix virtual void setupForTask( const unsigned& task_index, std::vector& indices, MultiValue& myvals ) const = 0; /// Run over one row of the matrix @@ -85,37 +54,8 @@ class ActionWithMatrix : public ActionWithVector { /// Check if there are forces we need to account for on this task bool checkForTaskForce( const unsigned& itask, const Value* myval ) const override ; /// This gathers the force on a particular value - virtual void gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const; + void gatherForces( const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const override ; }; -inline -bool ActionWithMatrix::matrixChainContinues() const { - return matrix_to_do_after!=NULL; -} - -inline -double ActionWithMatrix::getArgumentElement( const unsigned& ic, const unsigned& jelem, const MultiValue& myvals ) const { - return getPntrToArgument(ic)->get( jelem ); -} - -inline -double ActionWithMatrix::getElementOfMatrixArgument( const unsigned& imat, const unsigned& irow, const unsigned& jcol, const MultiValue& myvals ) const { - plumed_dbg_assert( imatgetRank()==2 && !getPntrToArgument(imat)->hasDerivatives() ); - return getArgumentElement( imat, irow*getPntrToArgument(imat)->getShape()[1] + jcol, myvals ); -} - -inline -void ActionWithMatrix::addDerivativeOnVectorArgument( const bool& inchain, const unsigned& ival, const unsigned& jarg, const unsigned& jelem, const double& der, MultiValue& myvals ) const { - plumed_dbg_massert( jarggetRank()<2, "failing in action " + getName() + " with label " + getLabel() ); - unsigned vstart=arg_deriv_starts[jarg]; myvals.addDerivative( ival, vstart + jelem, der ); myvals.updateIndex( ival, vstart + jelem ); -} - -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( jarggetRank()==2 && !getPntrToArgument(jarg)->hasDerivatives() ); - unsigned dloc = arg_deriv_starts[jarg] + irow*getPntrToArgument(jarg)->getNumberOfColumns() + jcol; - myvals.addDerivative( ival, dloc, der ); myvals.updateIndex( ival, dloc ); -} - } #endif diff --git a/src/core/ActionWithVector.cpp b/src/core/ActionWithVector.cpp index 6f9c98c368..91ce2da7cd 100644 --- a/src/core/ActionWithVector.cpp +++ b/src/core/ActionWithVector.cpp @@ -28,25 +28,6 @@ namespace PLMD { -enum class Option { no, yes }; - -Option interpretEnvString(const char* env,const char* str) { - if(!str) return Option::yes; - if(!std::strcmp(str,"yes"))return Option::yes; - if(!std::strcmp(str,"no"))return Option::no; - plumed_error()<<"Cannot understand env var "<( getPntrToArgument(i)->getPntrToAction() ); @@ -106,10 +82,6 @@ ActionWithVector::ActionWithVector(const ActionOptions&ao): } } -ActionWithVector::~ActionWithVector() { - if( action_to_do_before ) action_to_do_before->action_to_do_after=NULL; -} - void ActionWithVector::lockRequests() { ActionAtomistic::lockRequests(); ActionWithArguments::lockRequests(); @@ -124,181 +96,8 @@ void ActionWithVector::calculateNumericalDerivatives(ActionWithValue* av) { plumed_merror("cannot calculate numerical derivative for action " + getName() + " with label " + getLabel() ); } -void ActionWithVector::clearDerivatives( const bool& force ) { - if( !force && actionInChain() ) return; - ActionWithValue::clearDerivatives(); - if( action_to_do_after ) action_to_do_after->clearDerivatives( true ); -} - -void ActionWithVector::clearInputForces( const bool& force ) { - if( !force && actionInChain() ) return; - ActionWithValue::clearInputForces(); - if( action_to_do_after ) action_to_do_after->clearInputForces( true ); -} - -const ActionWithVector* ActionWithVector::getFirstActionInChain() const { - if( !actionInChain() ) return this; - return action_to_do_before->getFirstActionInChain(); -} - -ActionWithVector* ActionWithVector::getFirstActionInChain() { - if( !actionInChain() ) return this; - return action_to_do_before->getFirstActionInChain(); -} - -void ActionWithVector::retrieveAtoms( const bool& force ) { - if( !force && actionInChain() || atomsWereRetrieved ) return; - ActionAtomistic::retrieveAtoms(); atomsWereRetrieved = !actionInChain(); - if( action_to_do_after ) action_to_do_after->retrieveAtoms( true ); -} - -bool ActionWithVector::hasStoredArguments() const { - std::string headstr=getFirstActionInChain()->getLabel(); - for(unsigned i=0; iignoreStoredValue(headstr) ) return true; - } - return false; -} - -bool ActionWithVector::argumentDependsOn( const std::string& headstr, ActionWithVector* faction, Value* thearg ) { - for(unsigned i=0; i( getPntrToArgument(i)->getPntrToAction() ); - if( av && (av->getFirstActionInChain())->getLabel()==headstr ) { - if( av->argumentDependsOn( headstr, faction, thearg ) ) return true;; - } - } - return false; -} - -unsigned ActionWithVector::buildArgumentStore( const unsigned& argstart ) { - return reallyBuildArgumentStore( argstart ); -} - -unsigned ActionWithVector::reallyBuildArgumentStore( const unsigned& argstart ) { - for(unsigned i=argstart; igetRank()>0 ) getPntrToArgument(i)->buildDataStore(); } - unsigned nder=0; arg_deriv_starts.resize( getNumberOfArguments() ); - for(unsigned i=0; igetNumberOfValues(); } - return nder; -} - -ActionWithVector* ActionWithVector::getActionWithDerivatives( ActionWithVector* depaction ) { - if( depaction==this || depaction->checkForDependency(this) ) { - if( getNumberOfAtoms()>0 ) return this; - std::string c=getFirstActionInChain()->getLabel(); - for(unsigned i=0; iignoreStoredValue(c) && !getPntrToArgument(i)->isConstant() ) return this; - } - } - plumed_assert( action_to_do_before ); - return action_to_do_before->getActionWithDerivatives(depaction); -} - -bool ActionWithVector::addActionToChain( const std::vector& alabels, ActionWithVector* act ) { - if( action_to_do_after ) { bool state=action_to_do_after->addActionToChain( alabels, act ); return state; } - // Check action is not already in chain - std::vector mylabels; getFirstActionInChain()->getAllActionLabelsInChain( mylabels ); - for(unsigned i=0; igetLabel()==mylabels[i] ) return true; - } - - // Check that everything that is required has been calculated - for(unsigned i=0; i( alabels[i] ); - plumed_massert( av, "could not cast " + alabels[i] ); bool storingall=true; - for(int j=0; jgetNumberOfComponents(); ++j) { - if( !(av->getPntrToComponent(j))->storedata ) storingall=false; - } - if( !storingall ) return false; - } - } - // This checks that there is nothing that will cause problems in the chain - mylabels.resize(0); getFirstActionInChain()->getAllActionLabelsInChain( mylabels ); - for(unsigned i=0; i( mylabels[i] ); - for(unsigned j=0; j( mylabels[j] ); - if( !av1->canBeAfterInChain( av2 ) ) error("must calculate " + mylabels[j] + " before " + mylabels[i] ); - } - } - action_to_do_after=act; act->action_to_do_before=this; updateTaskListReductionStatus(); - ActionWithVector* head = getFirstActionInChain(); - head->broadcastThatTasksAreReduced( head ); head->finishChainBuild( act ); - return true; -} - -void ActionWithVector::updateTaskListReductionStatus() { - ActionWithVector* head = getFirstActionInChain(); - std::vector task_reducing_actions; head->canReduceTasks( task_reducing_actions ); - if( task_reducing_actions.size()>0 ) head->reduce_tasks=true; -} - -void ActionWithVector::broadcastThatTasksAreReduced( ActionWithVector* aselect ) { - std::string c=getFirstActionInChain()->getLabel(); - for(unsigned i=0; iignoreStoredValue(c) ) { - ActionWithVector* av = dynamic_cast( getPntrToArgument(i)->getPntrToAction() ); - if( av ) { - bool found=false; - ActionWithVector* av_head = av->getFirstActionInChain(); - for(unsigned i=0; itask_control_list.size(); ++i) { - if( aselect==av_head->task_control_list[i] ) { found=true; break; } - } - if( !found ) av_head->task_control_list.insert( av_head->task_control_list.begin(), aselect ); - - av_head->reduce_tasks=true; av_head->updateTaskReductionFlag( av_head->reduce_tasks ); - } - } - } - if( action_to_do_after ) action_to_do_after->broadcastThatTasksAreReduced( aselect ); -} - -void ActionWithVector::updateTaskReductionFlag( bool& head_reduce_tasks ) { - if( actionInChain() ) { - plumed_assert( task_control_list.size()==0 ); - } else { - for(unsigned i=0; igetFirstActionInChain())->reduce_tasks ) head_reduce_tasks=false; - } - } - broadcastThatTasksAreReduced( getFirstActionInChain() ); - if( action_to_do_after ) action_to_do_after->updateTaskReductionFlag( head_reduce_tasks ); -} - -void ActionWithVector::canReduceTasks( std::vector& task_reducing_actions ) { - areAllTasksRequired( task_reducing_actions ); - if( action_to_do_after ) action_to_do_after->canReduceTasks( task_reducing_actions ); -} - -void ActionWithVector::finishChainBuild( ActionWithVector* act ) { - if( action_to_do_after ) action_to_do_after->finishChainBuild( act ); -} - -void ActionWithVector::getAllActionLabelsInChain( std::vector& mylabels ) const { - bool found = false ; - for(unsigned i=0; igetAllActionLabelsInChain( mylabels ); -} - -void ActionWithVector::taskIsActive( const unsigned& current, int& flag ) const { - flag = checkTaskStatus( current, flag ); - if( flag<=0 && action_to_do_after ) action_to_do_after->taskIsActive( current, flag ); -} - -void ActionWithVector::getAdditionalTasksRequired( ActionWithVector* action, std::vector& atasks ) { - for(unsigned i=0; igetAdditionalTasksRequired( action, atasks ); -} - void ActionWithVector::prepare() { - active_tasks.resize(0); atomsWereRetrieved=false; + active_tasks.resize(0); } int ActionWithVector::checkTaskIsActive( const unsigned& itask ) const { @@ -347,64 +146,15 @@ std::vector& ActionWithVector::getListOfActiveTasks( ActionWithVector* if( active_tasks.size()>0 ) return active_tasks; unsigned ntasks=0; getNumberOfTasks( ntasks ); - if( getenvChainForbidden()==Option::yes ) { - std::vector taskFlags( ntasks, -1 ); - for(unsigned i=0; i0 ) nt++; - } - active_tasks.resize(nt); nt=0; - for(unsigned i=0; i0 ) { active_tasks[nt]=i; nt++; } - } - return active_tasks; + std::vector taskFlags( ntasks, -1 ); + for(unsigned i=0; i0 ) nt++; } - - unsigned stride=comm.Get_size(); - unsigned rank=comm.Get_rank(); - if(serial) { stride=1; rank=0; } - - // Get number of threads for OpenMP - unsigned nt=OpenMP::getNumThreads(); - if( nt*stride*10>ntasks ) nt=ntasks/stride/10; - if( nt==0 ) nt=1; - - if( !never_reduce_tasks && reduce_tasks ) { - if( task_control_list.size()>0 ) { - // Get the list of tasks that are active in the action that uses the output of this action - for(unsigned i=0; iretrieveAtoms(); - active_tasks = task_control_list[i]->getListOfActiveTasks( action ); - } - // Now work out else we need from here to calculate the later action - getAdditionalTasksRequired( action, active_tasks ); - } else { - std::vector taskFlags( ntasks, -1 ); - - #pragma omp parallel num_threads(nt) - { - #pragma omp for nowait - for(unsigned i=rank; i=stride ) nt++; - } - active_tasks.resize(nt); nt=0; - for(unsigned i=0; i=stride ) { active_tasks[nt]=i; nt++; } - } - getAdditionalTasksRequired( this, active_tasks ); - } - } else { - active_tasks.resize( ntasks ); - for(unsigned i=0; i0 ) { active_tasks[nt]=i; nt++; } } return active_tasks; } @@ -416,7 +166,7 @@ bool ActionWithVector::doNotCalculateDerivatives() const { void ActionWithVector::clearMatrixBookeeping() { for(unsigned i=0; istoredata ) continue; + Value* myval = getPntrToComponent(i); if( myval->getRank()==2 && myval->getNumberOfColumns()getShape()[1] ) { std::fill(myval->matrix_bookeeping.begin(), myval->matrix_bookeeping.end(), 0); } @@ -425,9 +175,6 @@ void ActionWithVector::clearMatrixBookeeping() { } void ActionWithVector::runAllTasks() { -// Skip this if this is done elsewhere - if( action_to_do_before ) return; - unsigned stride=comm.Get_size(); unsigned rank=comm.Get_rank(); if(serial) { stride=1; rank=0; } @@ -445,33 +192,33 @@ void ActionWithVector::runAllTasks() { if( nt*stride*10>nactive_tasks ) nt=nactive_tasks/stride/10; if( nt==0 ) nt=1; - // Now do all preparations required to run all the tasks - // prepareForTaskLoop(); - - if( !action_to_do_after ) { - forwardPass=true; - for(unsigned i=0; igetRank()==0 ) { forwardPass=false; break; } - } - } // Get the total number of streamed quantities that we need - unsigned nquants=0, nmatrices=0, maxcol=0; - setupStreamedComponents( getLabel(), nquants, nmatrices, maxcol ); // Get size for buffer - unsigned bufsize=0; getSizeOfBuffer( nactive_tasks, bufsize ); + unsigned bufsize=0, nderivatives = 0; bool gridsInStream=false; forwardPass=true; + for(int i=0; ibufstart=bufsize; + if( getPntrToComponent(i)->hasDerivatives() || getPntrToComponent(i)->getRank()==0 ) { + forwardPass=false; bufsize += getPntrToComponent(i)->data.size(); + } + if( getConstPntrToComponent(i)->getRank()>0 && getConstPntrToComponent(i)->hasDerivatives() ) { + nderivatives=getConstPntrToComponent(i)->getNumberOfGridDerivatives(); gridsInStream=true; + } + } if( buffer.size()!=bufsize ) buffer.resize( bufsize ); // Clear buffer buffer.assign( buffer.size(), 0.0 ); // Recover the number of derivatives we require - unsigned nderivatives = 0; bool gridsInStream=checkForGrids(nderivatives); - if( !doNotCalculateDerivatives() && !gridsInStream ) getNumberOfStreamedDerivatives( nderivatives, NULL ); + if( !doNotCalculateDerivatives() && !gridsInStream ) { + if( getNumberOfAtoms()>0 ) nderivatives += 3*getNumberOfAtoms() + 9; + for(unsigned i=0; igetNumberOfValues(); + } #pragma omp parallel num_threads(nt) { std::vector omp_buffer; if( nt>1 ) omp_buffer.resize( bufsize, 0.0 ); - MultiValue myvals( nquants, nderivatives, nmatrices, maxcol ); + MultiValue myvals( getNumberOfComponents(), nderivatives, 0 ); myvals.clearAll(); #pragma omp for nowait @@ -493,30 +240,14 @@ void ActionWithVector::runAllTasks() { // MPI Gather everything if( !serial ) { if( buffer.size()>0 ) comm.Sum( buffer ); - gatherProcesses(); - } - finishComputations( buffer ); forwardPass=false; -} - -void ActionWithVector::gatherProcesses() { - ActionWithMatrix* am=dynamic_cast(this); - for(unsigned i=0; istoredata && !myval->hasDeriv ) { - comm.Sum( myval->data ); if( am && myval->getRank()==2 && myval->getNumberOfColumns()getShape()[1] ) comm.Sum( myval->matrix_bookeeping ); - } - } - if( action_to_do_after ) action_to_do_after->gatherProcesses(); -} - -bool ActionWithVector::checkForGrids( unsigned& nder ) const { - for(int i=0; igetRank()>0 && getConstPntrToComponent(i)->hasDerivatives() ) { - nder=getConstPntrToComponent(i)->getNumberOfGridDerivatives(); return true; + for(unsigned i=0; ihasDeriv ) { + comm.Sum( myval->data ); if( am && myval->getRank()==2 && myval->getNumberOfColumns()getShape()[1] ) comm.Sum( myval->matrix_bookeeping ); + } } } - if( action_to_do_after ) return action_to_do_after->checkForGrids(nder); - return false; + finishComputations( buffer ); forwardPass=false; } void ActionWithVector::getNumberOfTasks( unsigned& ntasks ) { @@ -538,54 +269,6 @@ void ActionWithVector::getNumberOfTasks( unsigned& ntasks ) { } else if( getPntrToComponent(i)->hasDerivatives() && ntasks!=getPntrToComponent(i)->getNumberOfValues() ) error("mismatched numbers of tasks in streamed quantities"); else if( !getPntrToComponent(i)->hasDerivatives() && ntasks!=getPntrToComponent(i)->getShape()[0] ) error("mismatched numbers of tasks in streamed quantities"); } - if( action_to_do_after ) action_to_do_after->getNumberOfTasks( ntasks ); -} - -void ActionWithVector::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) { - for(int i=0; ibufstart=bufsize; - if( getPntrToComponent(i)->hasDerivatives() || getPntrToComponent(i)->getRank()==0 ) bufsize += getPntrToComponent(i)->data.size(); - } - if( action_to_do_after ) action_to_do_after->getSizeOfBuffer( nactive_tasks, bufsize ); -} - -void ActionWithVector::getNumberOfStreamedDerivatives( unsigned& nderivatives, Value* stopat ) { - std::string c=getFirstActionInChain()->getLabel(); - for(unsigned i=0; iignoreStoredValue(c) ) { - if( getPntrToArgument(i)==stopat ) return; - nderivatives += getPntrToArgument(i)->getNumberOfValues(); - } - } - if( getNumberOfAtoms()>0 ) nderivatives += 3*getNumberOfAtoms() + 9; - // Don't do the whole chain if we have been told to stop early - if( stopat && stopat->getPntrToAction()==this ) return; - - if( action_to_do_after ) action_to_do_after->getNumberOfStreamedDerivatives( nderivatives, stopat ); -} - -bool ActionWithVector::getNumberOfStoredValues( Value* startat, unsigned& nvals, const unsigned& astart, const std::vector& stopat ) { - for(unsigned j=astart; jgetPntrToAction()==this || (stopat[j]->getPntrToAction())->checkForDependency(this)) ) return true; - } - - std::string c=getFirstActionInChain()->getLabel(); - for(unsigned i=0; iignoreStoredValue(c) ) { - for(unsigned j=astart; jgetNumberOfValues(); - } - } - if( startat->getPntrToAction()!=this && getNumberOfAtoms()>0 ) return false; - - if( action_to_do_after ) return action_to_do_after->getNumberOfStoredValues( startat, nvals, astart, stopat ); - return false; } void ActionWithVector::runTask( const unsigned& current, MultiValue& myvals ) const { @@ -593,12 +276,11 @@ void ActionWithVector::runTask( const unsigned& current, MultiValue& myvals ) co myvals.setTaskIndex(current); myvals.vector_call=true; performTask( current, myvals ); for(unsigned i=0; ihasDerivatives() || !myval->valueIsStored() ) continue; + if( am || myval->hasDerivatives() ) continue; Value* myv = const_cast( myval ); if( getName()=="RMSD_VECTOR" && myv->getRank()==2 ) continue; myv->set( current, myvals.get( i ) ); } - if( action_to_do_after ) action_to_do_after->runTask( current, myvals ); } void ActionWithVector::gatherAccumulators( const unsigned& taskCode, const MultiValue& myvals, std::vector& buffer ) const { @@ -616,9 +298,8 @@ void ActionWithVector::gatherAccumulators( const unsigned& taskCode, const Multi } } // This looks after storing of vectors - } else if( getConstPntrToComponent(i)->storedata ) gatherStoredValue( i, taskCode, myvals, bufstart, buffer ); + } else gatherStoredValue( i, taskCode, myvals, bufstart, buffer ); } - if( action_to_do_after ) action_to_do_after->gatherAccumulators( taskCode, myvals, buffer ); } void ActionWithVector::finishComputations( const std::vector& buf ) { @@ -638,20 +319,17 @@ void ActionWithVector::finishComputations( const std::vector& buf ) { for(unsigned j=0; jgetNumberOfDerivatives(); ++j) getPntrToComponent(i)->setDerivative( j, buf[bufstart+1+j] ); } } - if( action_to_do_after ) action_to_do_after->finishComputations( buf ); } bool ActionWithVector::checkChainForNonScalarForces() const { for(unsigned i=0; igetRank()>0 && getConstPntrToComponent(i)->forcesWereAdded() ) return true; } - if( action_to_do_after ) return action_to_do_after->checkChainForNonScalarForces(); return false; } bool ActionWithVector::checkForForces() { if( getPntrToComponent(0)->getRank()==0 ) return ActionWithValue::checkForForces(); - else if( actionInChain() ) return false; // Check if there are any forces if( !checkChainForNonScalarForces() ) return false; @@ -662,13 +340,8 @@ bool ActionWithVector::checkForForces() { if(serial) { stride=1; rank=0; } // Get the number of tasks - std::vector force_tasks; - if( getenvChainForbidden()==Option::yes ) { - std::vector & partialTaskList( getListOfActiveTasks( this ) ); - for(unsigned i=0; i force_tasks; std::vector & partialTaskList( getListOfActiveTasks( this ) ); + for(unsigned i=0; i0 ) { - nderiv = 0; - for(unsigned i=0; igetNumberOfStoredValues(); + unsigned nmatrices=0; ActionWithMatrix* am=dynamic_cast(this); + if(am) { + for(unsigned i=0; igetRank()==2 && !getConstPntrToComponent(i)->hasDerivatives() ) { nmatrices=getNumberOfComponents(); } } - ActionAtomistic* aa=castToActionAtomistic(); - if(aa) nderiv += 3*aa->getNumberOfAtoms() + 9; + } + // Recover the number of derivatives we require (this should be equal to the number of forces) + unsigned nderiv=0; + if( getNumberOfAtoms()>0 ) nderiv += 3*getNumberOfAtoms() + 9; + for(unsigned i=0; igetNumberOfStoredValues(); } if( forcesForApply.size()!=nderiv ) forcesForApply.resize( nderiv ); // Clear force buffer @@ -698,7 +371,7 @@ bool ActionWithVector::checkForForces() { { std::vector omp_forces; if( nt>1 ) omp_forces.resize( forcesForApply.size(), 0.0 ); - MultiValue myvals( nquants, nderiv, nmatrices, maxcol ); + MultiValue myvals( getNumberOfComponents(), nderiv, nmatrices ); myvals.clearAll(); #pragma omp for nowait @@ -731,27 +404,8 @@ bool ActionWithVector::checkForTaskForce( const unsigned& itask, const Value* my return fabs(myval->getForce(itask))>epsilon; } -void ActionWithVector::updateForceTasksFromValue( const Value* myval, std::vector& force_tasks ) const { - if( myval->getRank()>0 && myval->forcesWereAdded() ) { - unsigned nt = myval->getNumberOfValues(); - if( !myval->hasDerivatives() ) nt = myval->getShape()[0]; - for(unsigned i=0; i& force_tasks ) const { - if( checkComponentsForForce() ) { - for(unsigned k=0; kgetForceTasks( force_tasks ); -} - void ActionWithVector::gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const { - const Value* myval = getConstPntrToComponent(ival); plumed_dbg_assert( myval->storedata ); + const Value* myval = getConstPntrToComponent(ival); double fforce = myval->getForce(itask); for(unsigned j=0; jgetRank()>0 && myval->forcesWereAdded() ) gatherForcesOnStoredValue( k, itask, myvals, forces ); } } - if( action_to_do_after ) action_to_do_after->gatherForces( itask, myvals, forces ); } void ActionWithVector::apply() { diff --git a/src/core/ActionWithVector.h b/src/core/ActionWithVector.h index 0f50d62668..bf081de186 100644 --- a/src/core/ActionWithVector.h +++ b/src/core/ActionWithVector.h @@ -47,119 +47,49 @@ class ActionWithVector: std::vector buffer; /// The list of active tasks std::vector active_tasks; - /// Action that must be done before this one - ActionWithVector* action_to_do_before; -/// Actions that must be done after this one - ActionWithVector* action_to_do_after; -/// This is the list of actions that control the tasks that we do here - std::vector task_control_list; -/// Work backwards through the chain to find an action that has either stored arguments or derivatives - ActionWithVector* getActionWithDerivatives( ActionWithVector* depaction ); -/// Gather all the data from the MPI processes - void gatherProcesses(); /// Clear all the bookeeping arrays void clearMatrixBookeeping(); -/// Check if there are any grids in the stream - bool checkForGrids(unsigned& nder) const ; /// Run the task void runTask( const unsigned& taskno, MultiValue& myvals ) const ; /// Gather the values that we intend to store in the buffer void gatherAccumulators( const unsigned& taskCode, const MultiValue& myvals, std::vector& buffer ) const ; -/// Gather the forces on non-scalar quantities - void gatherForces( const unsigned& i, const MultiValue& myvals, std::vector& forces ) const ; -/// Get the size of the buffer array that holds the data we are gathering over the MPI loop - void getSizeOfBuffer( const unsigned& nactive_tasks, unsigned& bufsize ); /// Get the number of stored values in the stream bool getNumberOfStoredValues( Value* startat, unsigned& nvals, const unsigned& astart, const std::vector& stopat ); -/// Add this action to the recursive chain - bool addActionToChain( const std::vector& alabels, ActionWithVector* act ); /// Check the chain for non scalar forces bool checkChainForNonScalarForces() const ; -/// Check if a force has been added on one of the components of this action - bool checkComponentsForForce() const ; -/// Get the tasks that we need for forces - void getForceTasks( std::vector& force_tasks ) const ; -/// Check if this ation can reduce the number of tasks we perform - void canReduceTasks( std::vector& task_reducing_actions ); -/// Send information to arguments that tasks are reduced in depedent actions - void broadcastThatTasksAreReduced( ActionWithVector* aselect ); -/// Turn on task reduction flag in dependent actions - void updateTaskReductionFlag( bool& head_reduce_tasks ); -/// Check if a particular task is active at this time - void taskIsActive( const unsigned& current, int& flag ) const ; -/// This is turned on if there is some action that needs all the tasks - bool never_reduce_tasks; -/// Are we allowed to reduce the number of tasks being performed - bool reduce_tasks; -/// Were the atoms retrieved in some earlier action - bool atomsWereRetrieved; -/// This is used to build the argument store when we cannot use the chain - unsigned reallyBuildArgumentStore( const unsigned& argstart ); protected: -/// A vector that contains the start point for the argument derivatives - std::vector arg_deriv_starts; /// Turn off the flag that says this action has a masked input void ignoreMaskArguments(); -/// This updates whether or not we are using all the task reduction stuff - void updateTaskListReductionStatus(); /// Run all calculations in serial bool runInSerial() const ; -/// Check if the arguments of this action depend on thearg - bool argumentDependsOn( const std::string& headstr, ActionWithVector* faction, Value* thearg ); -/// This sets up the arguments at the start of the calculation - unsigned buildArgumentStore( const unsigned& argstart ); /// Run all the tasks in the list void runAllTasks(); +/// Check if a force has been added on one of the components of this action + bool checkComponentsForForce() const ; /// Accumulate the forces from the Values bool checkForForces(); +/// Gather the forces on a particular value + void gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const ; public: static void registerKeywords( Keywords& keys ); explicit ActionWithVector(const ActionOptions&); - virtual ~ActionWithVector(); + virtual ~ActionWithVector() {} void lockRequests() override; void unlockRequests() override; virtual void prepare() override; - void retrieveAtoms( const bool& force=false ) override; /// Check if a mask has been set int getNumberOfMasks() const ; void calculateNumericalDerivatives(ActionWithValue* av) override; /// Turn off the calculation of the derivatives during the forward pass through a calculation bool doNotCalculateDerivatives() const override ; -/// Are we running this command in a chain - bool actionInChain() const ; -/// This is overwritten within ActionWithMatrix and is used to build the chain of just matrix actions - virtual void finishChainBuild( ActionWithVector* act ); -/// Check if there are any stored values in arguments - bool hasStoredArguments() const ; -/// Return a pointer to the first action in the chain - const ActionWithVector* getFirstActionInChain() const ; - ActionWithVector* getFirstActionInChain(); /// Get the list of tasks that are active virtual std::vector& getListOfActiveTasks( ActionWithVector* action ); -/// This is overridden in ActionWithMatrix - virtual void getAllActionLabelsInMatrixChain( std::vector& matchain ) const {} -/// Get the number of derivatives in the stream - virtual void getNumberOfStreamedDerivatives( unsigned& nderivatives, Value* stopat ); -/// Get every the label of every value that is calculated in this chain - void getAllActionLabelsInChain( std::vector& mylabels ) const ; -/// We override clearInputForces here to ensure that forces are deleted from all values - void clearInputForces( const bool& force=false ) override; -/// We override clearDerivatives here to prevent data in streams from being deleted - void clearDerivatives( const bool& force=false ) override; -/// Check if we can be after another ActionWithVector - virtual bool canBeAfterInChain( ActionWithVector* av ) { return true; } /// Do we always need to do all the tasks for this action virtual void areAllTasksRequired( std::vector& task_reducing_actions ) {} /// Find out how many tasks we need to perform in this loop virtual void getNumberOfTasks( unsigned& ntasks ); /// Check the status of the ith task virtual int checkTaskStatus( const unsigned& taskno, int& flag ) const { return flag; } -/// Check if we are in a subchain - virtual bool isInSubChain( unsigned& nder ) { return false; } -/// Get the additional tasks that are required here - virtual void getAdditionalTasksRequired( ActionWithVector* action, std::vector& atasks ); -/// setup the streamed quantities - virtual void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ); /// Determine if a particular task is active based on the values of the input argument virtual int checkTaskIsActive( const unsigned& itask ) const ; /// This we override to perform each individual task @@ -170,23 +100,16 @@ class ActionWithVector: virtual void switchTaskReduction( const bool& task_reduction, ActionWithVector* aselect ) {} /// Gather the values that we intend to store in the buffer virtual void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals, const unsigned& bufstart, std::vector& buffer ) const {} -/// Get the force tasks that are active for this action - virtual void updateForceTasksFromValue( const Value* myval, std::vector& force_tasks ) const ; /// Check if there is a force that needs to be accumulated on the ith task virtual bool checkForTaskForce( const unsigned& itask, const Value* myval ) const ; -/// Gather the forces on a particular value - virtual void gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const ; +/// Gather the forces on non-scalar quantities + virtual void gatherForces( const unsigned& i, const MultiValue& myvals, std::vector& forces ) const ; /// This is to transfer data from the buffer to the final value void finishComputations( const std::vector& buf ); /// Apply the forces on this data virtual void apply(); }; -inline -bool ActionWithVector::actionInChain() const { - return (action_to_do_before!=NULL); -} - inline bool ActionWithVector::runInSerial() const { return serial; diff --git a/src/core/PbcAction.cpp b/src/core/PbcAction.cpp index ccbca7ebbe..39af80c4f4 100644 --- a/src/core/PbcAction.cpp +++ b/src/core/PbcAction.cpp @@ -52,7 +52,7 @@ PbcAction::PbcAction(const ActionOptions&ao): { std::vector shape(2); shape[0]=shape[1]=3; addValue( shape ); setNotPeriodic(); setUnit( "length", "energy" ); - getPntrToValue()->buildDataStore(); getPntrToValue()->reshapeMatrixStore(3); + getPntrToValue()->reshapeMatrixStore(3); } diff --git a/src/core/Value.cpp b/src/core/Value.cpp index b3539269e4..1a447c2645 100644 --- a/src/core/Value.cpp +++ b/src/core/Value.cpp @@ -36,11 +36,9 @@ Value::Value(): action(NULL), value_set(false), hasForce(false), - storedata(false), shape(std::vector()), hasDeriv(true), bufstart(0), - matpos(0), ngrid_der(0), ncols(0), symmetric(false), @@ -59,12 +57,10 @@ Value::Value(const std::string& name): value_set(false), hasForce(false), name(name), - storedata(false), shape(std::vector()), hasDeriv(true), bufstart(0), ngrid_der(0), - matpos(0), ncols(0), symmetric(false), periodicity(unset), @@ -83,11 +79,9 @@ Value::Value(ActionWithValue* av, const std::string& name, const bool withderiv, value_set(false), hasForce(false), name(name), - storedata(false), hasDeriv(withderiv), bufstart(0), ngrid_der(0), - matpos(0), ncols(0), symmetric(false), periodicity(unset), @@ -100,8 +94,6 @@ Value::Value(ActionWithValue* av, const std::string& name, const bool withderiv, if( action ) { if( action->getName()=="ACCUMULATE" || action->getName()=="COLLECT" ) valtype=average; } - if( action ) storedata=action->getName()=="PUT" || valtype==average; - if( ss.size() && withderiv ) storedata=true; setShape( ss ); } @@ -126,7 +118,7 @@ void Value::setShape( const std::vector&ss ) { } else if( shape.size()==0 ) { // This is for scalars data.resize(1); inputForce.resize(1); - } else if( storedata && shape.size()<2 ) { + } else if( shape.size()<2 ) { // This is for vectors (matrices have special version because we have sparse storage) data.resize( tot ); inputForce.resize( tot ); } @@ -278,17 +270,8 @@ void Value::addForce(const std::size_t& iforce, double f, const bool trueind) { } -void Value::buildDataStore( const bool forprint ) { - if( getRank()==0 ) return; - storedata=true; setShape( shape ); - if( !forprint ) return ; - ActionWithVector* av=dynamic_cast( action ); - if( av ) (av->getFirstActionInChain())->never_reduce_tasks=true; -} - void Value::reshapeMatrixStore( const unsigned& n ) { plumed_dbg_assert( shape.size()==2 && !hasDeriv ); - if( !storedata ) return ; ncols=n; if( ncols>shape[1] ) ncols=shape[1]; unsigned size=shape[0]*ncols; if( matrix_bookeeping.size()!=(size+shape[0]) ) { @@ -310,20 +293,12 @@ void Value::copyBookeepingArrayFromArgument( Value* myarg ) { data.resize( shape[0]*ncols ); inputForce.resize( shape[0]*ncols ); } -void Value::setPositionInMatrixStash( const unsigned& p ) { - plumed_dbg_assert( shape.size()==2 && !hasDeriv ); - matpos=p; -} - bool Value::ignoreStoredValue(const std::string& c) const { - if( !storedata && shape.size()>0 ) return true; - ActionWithVector* av=dynamic_cast(action); - if( av ) return (av->getFirstActionInChain())->getLabel()==c; return false; } void Value::setConstant() { - valtype=constant; storedata=true; setShape( shape ); + valtype=constant; setShape( shape ); if( getRank()==2 && !hasDeriv ) reshapeMatrixStore( shape[1] ); } diff --git a/src/core/Value.h b/src/core/Value.h index 59ac360ed5..6b3e7d4ed5 100644 --- a/src/core/Value.h +++ b/src/core/Value.h @@ -80,14 +80,12 @@ class Value { std::map gradients; /// The name of this quantiy std::string name; -/// Are we storing the data for this value if it is vector or matrix - bool storedata; /// What is the shape of the value (0 dimensional=scalar, n dimensional with derivatives=grid, 1 dimensional no derivatives=vector, 2 dimensional no derivatives=matrix) std::vector shape; /// Does this quanity have derivatives bool hasDeriv; /// Variables for storing data - unsigned bufstart, matpos, ngrid_der, ncols; + unsigned bufstart, ngrid_der, ncols; /// If we are storing a matrix is it symmetric? bool symmetric; /// This is a bookeeping array that holds the non-zero elements of the "sparse" matrix @@ -187,8 +185,6 @@ class Value { unsigned getRank() const ; /// Get the shape of the object that is contained in this value const std::vector& getShape() const ; -/// This turns on storing of vectors/matrices - void buildDataStore( const bool forprint=false ); /// Reshape the storage for sparse matrices void reshapeMatrixStore( const unsigned& n ); /// Copy the matrix bookeeping stuff @@ -213,9 +209,6 @@ class Value { void setDerivativeIsZeroWhenValueIsZero(); /// Return a bool that tells us if the derivative is zero when the value is zero bool isDerivativeZeroWhenValueIsZero() const ; -/// This stuff handles where to look for the start of the row that contains the row of the matrix - void setPositionInMatrixStash( const unsigned& p ); - unsigned getPositionInMatrixStash() const ; /// Convert the input index to its corresponding indices void convertIndexToindices(const std::size_t& index, std::vector& indices ) const ; /// Print out all the values in this Value @@ -228,8 +221,6 @@ class Value { unsigned getRowLength( const unsigned& irow ) const ; /// unsigned getRowIndex( const unsigned& irow, const unsigned& jind ) const ; -/// Are we storing this value - bool valueIsStored() const ; /// unsigned getNumberOfColumns() const ; /// @@ -432,22 +423,12 @@ bool Value::isDerivativeZeroWhenValueIsZero() const { return derivativeIsZeroWhenValueIsZero; } -inline -unsigned Value::getPositionInMatrixStash() const { - return matpos; -} - inline void Value::setMatrixBookeepingElement( const unsigned& i, const unsigned& n ) { plumed_dbg_assert( igetNumberOfValues(); for(unsigned i=0; i<4; ++i) { - Value* myarg=getPntrToArgument(i); myarg->buildDataStore(); + Value* myarg=getPntrToArgument(i); if( myarg->getRank()!=1 ) error("first four arguments to this action should be vectors"); if( (myarg->getPntrToAction())->getName()!="QUATERNION_VECTOR" ) error("first four arguments to this action should be quaternions"); std::string mylab=getPntrToArgument(i)->getName(); std::size_t dot=mylab.find_first_of("."); @@ -110,7 +110,6 @@ QuaternionBondProductMatrix::QuaternionBondProductMatrix(const ActionOptions&ao) addComponent( "i", shape ); componentIsNotPeriodic("i"); addComponent( "j", shape ); componentIsNotPeriodic("j"); addComponent( "k", shape ); componentIsNotPeriodic("k"); - unsigned nderivatives = buildArgumentStore(0); } unsigned QuaternionBondProductMatrix::getNumberOfDerivatives() { diff --git a/src/crystdistrib/QuaternionProductMatrix.cpp b/src/crystdistrib/QuaternionProductMatrix.cpp index de87a7368b..103da7eaef 100644 --- a/src/crystdistrib/QuaternionProductMatrix.cpp +++ b/src/crystdistrib/QuaternionProductMatrix.cpp @@ -36,7 +36,7 @@ namespace crystdistrib { class QuaternionProductMatrix : public ActionWithMatrix { private: - unsigned nderivatives; + void addDerivativeOnVectorArgument( const unsigned& ival, const unsigned& jarg, const unsigned& jelem, const double& der, MultiValue& myvals ) const ; public: static void registerKeywords( Keywords& keys ); explicit QuaternionProductMatrix(const ActionOptions&); @@ -79,11 +79,12 @@ QuaternionProductMatrix::QuaternionProductMatrix(const ActionOptions&ao): addComponent( "i", shape ); componentIsNotPeriodic("i"); addComponent( "j", shape ); componentIsNotPeriodic("j"); addComponent( "k", shape ); componentIsNotPeriodic("k"); - nderivatives = buildArgumentStore(0); } unsigned QuaternionProductMatrix::getNumberOfDerivatives() { - return nderivatives; + unsigned nder=0; + for(unsigned i=0; i<8; ++i) nder += getPntrToArgument(i)->getNumberOfStoredValues(); + return nder; } unsigned QuaternionProductMatrix::getNumberOfColumns() const { @@ -105,6 +106,12 @@ void QuaternionProductMatrix::setupForTask( const unsigned& task_index, std::vec myvals.setSplitIndex( size_v + 1 ); } +void QuaternionProductMatrix::addDerivativeOnVectorArgument( const unsigned& ival, const unsigned& jarg, const unsigned& jelem, const double& der, MultiValue& myvals ) const { + plumed_dbg_massert( jarggetRank()<2, "failing in action " + getName() + " with label " + getLabel() ); + unsigned vstart=0; for(unsigned i=0; igetNumberOfStoredValues(); + myvals.addDerivative( ival, vstart + jelem, der ); myvals.updateIndex( ival, vstart + jelem ); +} + void QuaternionProductMatrix::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const { unsigned ostrn, ind2=index2; if( index2>=getPntrToArgument(0)->getShape()[0] ) ind2 = index2 - getPntrToArgument(0)->getShape()[0]; @@ -112,9 +119,9 @@ void QuaternionProductMatrix::performTask( const std::string& controller, const std::vector quat1(4), quat2(4); // Retrieve the first quaternion - for(unsigned i=0; i<4; ++i) quat1[i] = getArgumentElement( i, index1, myvals ); + for(unsigned i=0; i<4; ++i) quat1[i] = getPntrToArgument(i)->get( index1 ); // Retrieve the second quaternion - for(unsigned i=0; i<4; ++i) quat2[i] = getArgumentElement( 4+i, ind2, myvals ); + for(unsigned i=0; i<4; ++i) quat2[i] = getPntrToArgument(4+i)->get( ind2 ); //make q1 the conjugate quat1[1] *= -1; @@ -131,8 +138,8 @@ void QuaternionProductMatrix::performTask( const std::string& controller, const myvals.addValue( 0, pref*quat1[i]*quat2[i] ); if( doNotCalculateDerivatives() ) continue ; if (i>0) conj=-1; - addDerivativeOnVectorArgument( false, 0, i, index1, conj*pref*quat2[i], myvals ); - addDerivativeOnVectorArgument( false, 0, 4+i, ind2, pref2*quat1[i], myvals ); + addDerivativeOnVectorArgument( 0, i, index1, conj*pref*quat2[i], myvals ); + addDerivativeOnVectorArgument( 0, 4+i, ind2, pref2*quat1[i], myvals ); } //i component pref=1; @@ -146,8 +153,8 @@ void QuaternionProductMatrix::performTask( const std::string& controller, const myvals.addValue( 1, pref*quat1[i]*quat2[(5-i)%4]); if( doNotCalculateDerivatives() ) continue ; if (i>0) conj=-1; - addDerivativeOnVectorArgument( false, 1, i, index1, conj*pref*quat2[(5-i)%4], myvals ); - addDerivativeOnVectorArgument( false, 1, 4+i, ind2, pref2*quat1[(5-i)%4], myvals ); + addDerivativeOnVectorArgument( 1, i, index1, conj*pref*quat2[(5-i)%4], myvals ); + addDerivativeOnVectorArgument( 1, 4+i, ind2, pref2*quat1[(5-i)%4], myvals ); } //j component @@ -162,8 +169,8 @@ void QuaternionProductMatrix::performTask( const std::string& controller, const myvals.addValue( 2, pref*quat1[i]*quat2[(i+2)%4]); if( doNotCalculateDerivatives() ) continue ; if (i>0) conj=-1; - addDerivativeOnVectorArgument( false, 2, i, index1, conj*pref*quat2[(i+2)%4], myvals ); - addDerivativeOnVectorArgument( false, 2, 4+i, ind2, pref2*quat1[(i+2)%4], myvals ); + addDerivativeOnVectorArgument( 2, i, index1, conj*pref*quat2[(i+2)%4], myvals ); + addDerivativeOnVectorArgument( 2, 4+i, ind2, pref2*quat1[(i+2)%4], myvals ); } //k component @@ -178,8 +185,8 @@ void QuaternionProductMatrix::performTask( const std::string& controller, const myvals.addValue( 3, pref*quat1[i]*quat2[(3-i)]); if( doNotCalculateDerivatives() ) continue ; if (i>0) conj=-1; - addDerivativeOnVectorArgument( false, 3, i, index1, conj*pref*quat2[3-i], myvals ); - addDerivativeOnVectorArgument( false, 3, 4+i, ind2, pref2*quat1[3-i], myvals ); + addDerivativeOnVectorArgument( 3, i, index1, conj*pref*quat2[3-i], myvals ); + addDerivativeOnVectorArgument( 3, 4+i, ind2, pref2*quat1[3-i], myvals ); } @@ -189,20 +196,18 @@ void QuaternionProductMatrix::performTask( const std::string& controller, const void QuaternionProductMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector & indices, MultiValue& myvals ) const { if( doNotCalculateDerivatives() ) return ; - for(unsigned j=0; jgetPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat ); - std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); unsigned ntwo_atoms = myvals.getSplitIndex(); - // Quaternion for first molecule - unsigned base = 0; for(unsigned k=0; k<4; ++k) { matrix_indices[nmat_ind] = base + ival; base += getPntrToArgument(k)->getShape()[0]; nmat_ind++; } - // Loop over row of matrix - for(unsigned i=1; i=getPntrToArgument(0)->getShape()[0] ) ind2 = indices[i] - getPntrToArgument(0)->getShape()[0]; - base = 4*getPntrToArgument(0)->getShape()[0]; - // Quaternion of second molecule - for(unsigned k=0; k<4; ++k) { matrix_indices[nmat_ind] = base + ind2; base += getPntrToArgument(4+k)->getShape()[0]; nmat_ind++; } - } - myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind ); + unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives(); + std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices() ); unsigned ntwo_atoms = myvals.getSplitIndex(); + // Quaternion for first molecule + unsigned base = 0; for(unsigned k=0; k<4; ++k) { matrix_indices[nmat_ind] = base + ival; base += getPntrToArgument(k)->getShape()[0]; nmat_ind++; } + // Loop over row of matrix + for(unsigned i=1; i=getPntrToArgument(0)->getShape()[0] ) ind2 = indices[i] - getPntrToArgument(0)->getShape()[0]; + base = 4*getPntrToArgument(0)->getShape()[0]; + // Quaternion of second molecule + for(unsigned k=0; k<4; ++k) { matrix_indices[nmat_ind] = base + ind2; base += getPntrToArgument(4+k)->getShape()[0]; nmat_ind++; } } + myvals.setNumberOfMatrixRowDerivatives( nmat_ind ); } diff --git a/src/dimred/ArrangePoints.cpp b/src/dimred/ArrangePoints.cpp index 357eefe4d0..25f6ead2f7 100644 --- a/src/dimred/ArrangePoints.cpp +++ b/src/dimred/ArrangePoints.cpp @@ -99,7 +99,7 @@ ArrangePoints::ArrangePoints( const ActionOptions& ao ) : for(unsigned i=0; igetNumberOfValues() ) error("mismatch between sizes of input coordinates"); std::string num; Tools::convert( i+1, num ); addComponent( "coord-" + num, shape ); - componentIsNotPeriodic( "coord-" + num ); getPntrToArgument(i)->buildDataStore(); + componentIsNotPeriodic( "coord-" + num ); } std::vector args( getArguments() ), target, weights; std::string sfd, errors; // Read in target "distances" and target weights @@ -154,7 +154,6 @@ ArrangePoints::ArrangePoints( const ActionOptions& ao ) : } void ArrangePoints::checkInputMatrix( const std::string& key, const unsigned& nvals, const std::vector& mat ) const { - mat[0]->buildDataStore(); if( mat.size()!=1 ) error("should only be one value in input to " + key ); if( mat[0]->getRank()!=2 || mat[0]->hasDerivatives() ) error("input to " + key + " keyword should be a matrix"); if( mat[0]->getShape()[0]!=nvals || mat[0]->getShape()[1]!=nvals ) error("input to " + key + " keyword has the wrong size"); diff --git a/src/dimred/ProjectPoints.cpp b/src/dimred/ProjectPoints.cpp index e9f4a39562..c5facbcbd5 100644 --- a/src/dimred/ProjectPoints.cpp +++ b/src/dimred/ProjectPoints.cpp @@ -95,7 +95,7 @@ ProjectPoints::ProjectPoints( const ActionOptions& ao ) : if( weights.size()!=1 ) error("should only be one value in input to WEIGHTS" + inum ); if( weights[0]->getRank()!=1 || weights[0]->hasDerivatives() ) error("input to WEIGHTS" + inum + " keyword should be a vector"); if( weights[0]->getShape()[0]!=nvals ) error("number of weights should match number of input coordinates"); - target[0]->buildDataStore(); weights[0]->buildDataStore(); args.push_back( target[0] ); args.push_back( weights[0] ); + args.push_back( target[0] ); args.push_back( weights[0] ); bool has_sf = parseNumbered("FUNC",i,sfd); switchingFunction.push_back( SwitchingFunction() ); if( !has_sf ) { switchingFunction[i-1].set( "CUSTOM FUNC=1-sqrt(x2) R_0=1.0", errors ); diff --git a/src/function/FunctionOfMatrix.h b/src/function/FunctionOfMatrix.h index 4340d0175e..79315026e4 100644 --- a/src/function/FunctionOfMatrix.h +++ b/src/function/FunctionOfMatrix.h @@ -129,8 +129,6 @@ FunctionOfMatrix::FunctionOfMatrix(const ActionOptions&ao): } // Set the periodicities of the output components myfunc.setPeriodicityForOutputs( this ); - // Now setup the action in the chain if we can - unsigned nderivatives = buildArgumentStore(myfunc.getArgStart()); } template diff --git a/src/function/FunctionOfVector.h b/src/function/FunctionOfVector.h index ced67d804d..f5d180e146 100644 --- a/src/function/FunctionOfVector.h +++ b/src/function/FunctionOfVector.h @@ -97,8 +97,8 @@ FunctionOfVector::FunctionOfVector(const ActionOptions&ao): myfunc.read( this ); // Create the task list if( myfunc.doWithTasks() ) { - doAtEnd=false; - } else { plumed_assert( getNumberOfArguments()==1 ); getPntrToArgument(0)->buildDataStore(); } + doAtEnd=false; + } else if( getNumberOfArguments()!=1 ) error("number of arguments should be equal to one"); // Get the names of the components std::vector components( keywords.getOutputComponents() ); // Create the values to hold the output @@ -133,8 +133,6 @@ FunctionOfVector::FunctionOfVector(const ActionOptions&ao): unsigned argstart=myfunc.getArgStart(); // Set the periodicities of the output components myfunc.setPeriodicityForOutputs( this ); - // Setup the derivatives - unsigned nderivatives = buildArgumentStore(myfunc.getArgStart()); } template diff --git a/src/generic/Collect.cpp b/src/generic/Collect.cpp index fe44c0344f..48fe1ccb98 100644 --- a/src/generic/Collect.cpp +++ b/src/generic/Collect.cpp @@ -98,7 +98,7 @@ Collect::Collect( const ActionOptions& ao ): nvals=(clearstride/getStride()); } - std::vector shape(1); shape[0]=nvals; getPntrToArgument(0)->buildDataStore(); + std::vector shape(1); shape[0]=nvals; if( type=="matrix" ) { shape.resize(2); shape[1] = getPntrToArgument(0)->getNumberOfValues(); } if( type=="vector" ) { shape[0] = nvals*getPntrToArgument(0)->getNumberOfValues(); } addValue( shape ); if( shape.size()==2 ) getPntrToComponent(0)->reshapeMatrixStore( shape[1] ); diff --git a/src/generic/CreateMask.cpp b/src/generic/CreateMask.cpp index 4177567ff3..36e2b1c8a8 100644 --- a/src/generic/CreateMask.cpp +++ b/src/generic/CreateMask.cpp @@ -79,11 +79,11 @@ CreateMask::CreateMask( const ActionOptions& ao ) : } else if( stype=="stride" ) { type=stride; log.printf(" setting every %d equally spaced points in output mask to zero \n", nzeros ); } else if( stype=="random" ) { - unsigned seed=230623; parse("SEED",seed); r.setSeed(-seed); getPntrToArgument(0)->buildDataStore(); + unsigned seed=230623; parse("SEED",seed); r.setSeed(-seed); type=random; log.printf(" choosing %d points to set to non-zero in mask in accordance with input weights \n", nzeros ); } else error( stype + " is not a valid way input for TYPE"); std::vector shape(1); shape[0] = getPntrToArgument(0)->getShape()[0]; - addValue( shape ); setNotPeriodic(); getPntrToComponent(0)->buildDataStore(); + addValue( shape ); setNotPeriodic(); for(unsigned i=0; iset( i, 1.0 ); } diff --git a/src/generic/DumpAtoms.cpp b/src/generic/DumpAtoms.cpp index 3db26f4a9f..c5c9be55e6 100644 --- a/src/generic/DumpAtoms.cpp +++ b/src/generic/DumpAtoms.cpp @@ -246,7 +246,7 @@ DumpAtoms::DumpAtoms(const ActionOptions&ao): for(unsigned i=0; igetRank()!=1 || getPntrToArgument(i)->hasDerivatives() ) error("arguments for xyz output should be vectors"); if( getPntrToArgument(i)->getNumberOfValues()!=atoms.size() ) error("number of elements in vector " + getPntrToArgument(i)->getName() + " is not equal to number of atoms output"); - getPntrToArgument(i)->buildDataStore(true); argnames.push_back( getPntrToArgument(i)->getName() ); + argnames.push_back( getPntrToArgument(i)->getName() ); } std::vector str_upper, str_lower; std::string errors; parseVector("LESS_THAN_OR_EQUAL",str_upper); parseVector("GREATER_THAN_OR_EQUAL",str_lower); diff --git a/src/generic/DumpPDB.cpp b/src/generic/DumpPDB.cpp index d94677cee6..736ddb73d0 100644 --- a/src/generic/DumpPDB.cpp +++ b/src/generic/DumpPDB.cpp @@ -129,7 +129,6 @@ void DumpPDB::buildArgnames() { } else if( getPntrToArgument(i)->getRank()==2 ) { (getPntrToArgument(i)->getPntrToAction())->getMatrixColumnTitles( argnames ); } - getPntrToArgument(i)->buildDataStore(); } } diff --git a/src/generic/DumpVector.cpp b/src/generic/DumpVector.cpp index 51f3b7a963..2805071fbb 100644 --- a/src/generic/DumpVector.cpp +++ b/src/generic/DumpVector.cpp @@ -94,7 +94,6 @@ void DumpVector::buildArgnames() { } else if( getPntrToArgument(i)->getRank()==2 ) { (getPntrToArgument(i)->getPntrToAction())->getMatrixColumnTitles( argnames ); } - getPntrToArgument(i)->buildDataStore(); } } diff --git a/src/generic/Print.cpp b/src/generic/Print.cpp index 4f3ca3bfdf..d057f08138 100644 --- a/src/generic/Print.cpp +++ b/src/generic/Print.cpp @@ -130,7 +130,6 @@ Print::Print(const ActionOptions&ao): log.printf(" with format %s\n",fmt.c_str()); for(unsigned i=0; ibuildDataStore(true); } ///////////////////////////////////////// // these are crazy things just for debug: diff --git a/src/generic/PrintNDX.cpp b/src/generic/PrintNDX.cpp index afe4a25116..3b2b58700a 100644 --- a/src/generic/PrintNDX.cpp +++ b/src/generic/PrintNDX.cpp @@ -97,7 +97,7 @@ PrintNDX::PrintNDX(const ActionOptions&ao): for(unsigned i=0; igetRank()!=1 || getPntrToArgument(i)->hasDerivatives() ) error("arguments for print ndx should be vector"); if( getPntrToArgument(i)->getShape()[0]!=all_atoms.size() ) error("mismatch between number of arguments and number of input atoms"); - getPntrToArgument(i)->buildDataStore(true); argnames[i] = getPntrToArgument(i)->getName(); + argnames[i] = getPntrToArgument(i)->getName(); } log.printf(" printing ndx file containing indices of atoms that have arguments in ranges prescribed below \n"); log.printf(" full set of atom indices investigated are : "); diff --git a/src/gridtools/ActionWithGrid.cpp b/src/gridtools/ActionWithGrid.cpp index 658e28c9c4..5400aa947e 100644 --- a/src/gridtools/ActionWithGrid.cpp +++ b/src/gridtools/ActionWithGrid.cpp @@ -45,7 +45,6 @@ ActionWithGrid::ActionWithGrid(const ActionOptions&ao): } void ActionWithGrid::calculate() { - plumed_assert( !actionInChain() ); if( firststep ) { setupOnFirstStep( true ); firststep=false; } runAllTasks(); diff --git a/src/gridtools/InterpolateGrid.cpp b/src/gridtools/InterpolateGrid.cpp index 73c72bd1a2..ec448a3958 100644 --- a/src/gridtools/InterpolateGrid.cpp +++ b/src/gridtools/InterpolateGrid.cpp @@ -68,7 +68,7 @@ class InterpolateGrid : public ActionWithGrid { void performTask( const unsigned& current, MultiValue& myvals ) const override ; void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals, const unsigned& bufstart, std::vector& buffer ) const ; - void gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const override ; + void gatherForces( const unsigned& i, const MultiValue& myvals, std::vector& forces ) const override ; }; PLUMED_REGISTER_ACTION(InterpolateGrid,"INTERPOLATE_GRID") @@ -167,14 +167,16 @@ void InterpolateGrid::performTask( const unsigned& current, MultiValue& myvals ) void InterpolateGrid::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals, const unsigned& bufstart, std::vector& buffer ) const { - plumed_dbg_assert( valindex==0 ); + plumed_dbg_assert( valindex==0 ); unsigned istart = bufstart + (1+output_grid.getDimension())*code; buffer[istart] += myvals.get( 0 ); for(unsigned i=0; i& forces ) const { - std::vector pos(output_grid.getDimension()); double ff = getConstPntrToComponent(ival)->getForce(itask); - output_grid.getGridPointCoordinates( itask, pos ); input_grid.applyForce( this, pos, ff, forces ); +void InterpolateGrid::gatherForces( const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const { + if( checkComponentsForForce() ) { + std::vector pos(output_grid.getDimension()); double ff = getConstPntrToComponent(0)->getForce(itask); + output_grid.getGridPointCoordinates( itask, pos ); input_grid.applyForce( this, pos, ff, forces ); + } } diff --git a/src/gridtools/KDE.cpp b/src/gridtools/KDE.cpp index 92c4c33cc7..8674f77eb1 100644 --- a/src/gridtools/KDE.cpp +++ b/src/gridtools/KDE.cpp @@ -88,8 +88,7 @@ class KDE : public ActionWithGrid { void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals, const unsigned& bufstart, std::vector& buffer ) const override ; bool checkForTaskForce( const unsigned& itask, const Value* myval ) const override ; - void updateForceTasksFromValue( const Value* myval, std::vector& force_tasks ) const override ; - void gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const override ; + void gatherForces( const unsigned& i, const MultiValue& myvals, std::vector& forces ) const override ; }; PLUMED_REGISTER_ACTION(KDE,"KDE") @@ -248,10 +247,8 @@ KDE::KDE(const ActionOptions&ao): if( ignore_out_of_bounds ) log.printf(" ignoring kernels that are outside of grid \n"); addValueWithDerivatives( shape ); setNotPeriodic(); getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero(); - // Make sure we store all the arguments - for(unsigned i=0; ibuildDataStore(); // Check for task reduction - updateTaskListReductionStatus(); setupOnFirstStep( false ); + setupOnFirstStep( false ); } void KDE::setupOnFirstStep( const bool incalc ) { @@ -534,27 +531,18 @@ void KDE::gatherStoredValue( const unsigned& valindex, const unsigned& code, con } } -void KDE::updateForceTasksFromValue( const Value* myval, std::vector& force_tasks ) const { - if( !myval->forcesWereAdded() ) return ; - if( numberOfKernels==1 ) plumed_error(); - - int flag=1; - for(unsigned i=0; iforcesWereAdded() ) return false; if( numberOfKernels==1 ) plumed_error(); return checkTaskIsActive( itask ); } -void KDE::gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const { +void KDE::gatherForces( const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const { if( numberOfKernels==1 ) { plumed_error(); return; } + if( !checkComponentsForForce() ) return; double height; std::vector args( gridobject.getDimension() ); retrieveArgumentsAndHeight( myvals, args, height ); unsigned num_neigh; std::vector neighbors; @@ -585,7 +573,7 @@ void KDE::gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask for(unsigned i=0; igetForce( neighbors[i] ); double newval = height*von_misses_norm*exp( von_misses_concentration*dot ); + double fforce = getConstPntrToComponent(0)->getForce( neighbors[i] ); double newval = height*von_misses_norm*exp( von_misses_concentration*dot ); if( hasheight && getPntrToArgument(args.size())->getRank()==0 ) forces[ hforce_start ] += newval*fforce / height; else if( hasheight ) forces[ hforce_start + getPntrToArgument(args.size())->getIndexInStore(itask) ] += newval*fforce / height; unsigned n=0; for(unsigned j=0; jgetIndexInStore(itask)] += von_misses_concentration*newval*gpoint[j]*fforce; n += getPntrToArgument(j)->getNumberOfStoredValues(); } diff --git a/src/landmarks/FarthestPointSampling.cpp b/src/landmarks/FarthestPointSampling.cpp index 119a982cce..093d50dac5 100644 --- a/src/landmarks/FarthestPointSampling.cpp +++ b/src/landmarks/FarthestPointSampling.cpp @@ -67,7 +67,7 @@ FarthestPointSampling::FarthestPointSampling( const ActionOptions& ao ): log.printf(" selecting %d landmark points \n", nlandmarks ); std::vector shape(1); shape[0] = getPntrToArgument(0)->getShape()[0]; - addValue( shape ); setNotPeriodic(); getPntrToComponent(0)->buildDataStore(); + addValue( shape ); setNotPeriodic(); } void FarthestPointSampling::prepare() { diff --git a/src/mapping/GeometricPath.cpp b/src/mapping/GeometricPath.cpp index bdd6d14acd..3719108509 100644 --- a/src/mapping/GeometricPath.cpp +++ b/src/mapping/GeometricPath.cpp @@ -63,7 +63,6 @@ GeometricPath::GeometricPath(const ActionOptions&ao): ActionWithVector(ao), path_projector(this) { - plumed_assert( !actionInChain() ); // Get the coordinates in the low dimensional space std::vector pcoord; parseVector("PROPERTY", pcoord ); std::vector theprop; ActionWithArguments::interpretArgumentList( pcoord, plumed.getActionSet(), this, theprop ); diff --git a/src/mapping/PathDisplacements.cpp b/src/mapping/PathDisplacements.cpp index a202dd8f77..797dddf28d 100644 --- a/src/mapping/PathDisplacements.cpp +++ b/src/mapping/PathDisplacements.cpp @@ -98,7 +98,6 @@ PathDisplacements::PathDisplacements(const ActionOptions& ao): // And create a value to hold the displacements std::vector shape(2); shape[0]=nrows; shape[1]=ncols; addValue( shape ); setNotPeriodic(); - getPntrToComponent(0)->buildDataStore(); getPntrToComponent(0)->reshapeMatrixStore( shape[1] ); } diff --git a/src/mapping/PathProjectionCalculator.cpp b/src/mapping/PathProjectionCalculator.cpp index fe22867246..40c22335b3 100644 --- a/src/mapping/PathProjectionCalculator.cpp +++ b/src/mapping/PathProjectionCalculator.cpp @@ -45,7 +45,6 @@ PathProjectionCalculator::PathProjectionCalculator( Action* act ): if( aarg->getNumberOfArguments()!=1 ) act->error("should only have one argument to this function"); } // Ensure that values are stored in base calculation and that PLUMED doesn't try to calculate this in the stream - if( mypath_obj ) mypath_obj->buildDataStore(); // Check that the input is a matrix if( mypath_obj ) if( mypath_obj->getRank()!=2 ) act->error("the input to this action should be a matrix"); // Get the labels for the reference points diff --git a/src/matrixtools/DiagonalizeMatrix.cpp b/src/matrixtools/DiagonalizeMatrix.cpp index 7a8cda7fef..5298c268f6 100644 --- a/src/matrixtools/DiagonalizeMatrix.cpp +++ b/src/matrixtools/DiagonalizeMatrix.cpp @@ -92,7 +92,6 @@ DiagonalizeMatrix::DiagonalizeMatrix(const ActionOptions& ao): addComponent( "vals-" + num, eval_shape ); componentIsNotPeriodic( "vals-" + num ); addComponent( "vecs-" + num, evec_shape ); componentIsNotPeriodic( "vecs-" + num ); // Make sure eigenvalues are always stored - getPntrToComponent( 2*i+1 )->buildDataStore(); } std::vector eigvecs_shape(2); eigvecs_shape[0]=eigvecs_shape[1]=getPntrToArgument(0)->getShape()[0]; diff --git a/src/matrixtools/InvertMatrix.cpp b/src/matrixtools/InvertMatrix.cpp index 72633d261a..8659766d2d 100644 --- a/src/matrixtools/InvertMatrix.cpp +++ b/src/matrixtools/InvertMatrix.cpp @@ -71,7 +71,7 @@ InvertMatrix::InvertMatrix(const ActionOptions& ao): if(as) input_is_constant=true; std::vector shape(2); shape[0]=shape[1]=getPntrToArgument(0)->getShape()[0]; addValue( shape ); - setNotPeriodic(); getPntrToComponent(0)->buildDataStore(); getPntrToComponent(0)->reshapeMatrixStore( shape[1] ); + setNotPeriodic(); getPntrToComponent(0)->reshapeMatrixStore( shape[1] ); mymatrix.resize( shape[0], shape[1] ); inverse.resize( shape[0], shape[1] ); } diff --git a/src/matrixtools/MatrixOperationBase.cpp b/src/matrixtools/MatrixOperationBase.cpp index e2b2a4d7b0..d62c406b1c 100644 --- a/src/matrixtools/MatrixOperationBase.cpp +++ b/src/matrixtools/MatrixOperationBase.cpp @@ -44,7 +44,6 @@ MatrixOperationBase::MatrixOperationBase(const ActionOptions&ao): if (getPntrToArgument(0)->getRank()!=1 || getPntrToArgument(0)->hasDerivatives() ) error("input to this argument should be a matrix or vector"); } else error("input to this argument should be a matrix"); } - getPntrToArgument(0)->buildDataStore(); } void MatrixOperationBase::retrieveFullMatrix( Matrix& mymatrix ) { diff --git a/src/matrixtools/MatrixTimesMatrix.cpp b/src/matrixtools/MatrixTimesMatrix.cpp index 7916dd868c..9d3a275c02 100644 --- a/src/matrixtools/MatrixTimesMatrix.cpp +++ b/src/matrixtools/MatrixTimesMatrix.cpp @@ -47,15 +47,12 @@ class MatrixTimesMatrix : public ActionWithMatrix { private: bool squared; bool diagzero; - unsigned nderivatives; - bool stored_matrix1, stored_matrix2; public: static void registerKeywords( Keywords& keys ); explicit MatrixTimesMatrix(const ActionOptions&); void prepare() override ; unsigned getNumberOfDerivatives(); unsigned getNumberOfColumns() const override ; - void getAdditionalTasksRequired( ActionWithVector* action, std::vector& atasks ) override ; void setupForTask( const unsigned& task_index, std::vector& indices, MultiValue& myvals ) const ; void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override; void runEndOfRowJobs( const unsigned& ival, const std::vector & indices, MultiValue& myvals ) const override ; @@ -81,10 +78,7 @@ MatrixTimesMatrix::MatrixTimesMatrix(const ActionOptions&ao): if( getPntrToArgument(1)->getRank()!=2 || getPntrToArgument(1)->hasDerivatives() ) error("second argument to this action should be a matrix"); if( getPntrToArgument(0)->getShape()[1]!=getPntrToArgument(1)->getShape()[0] ) error("number of columns in first matrix does not equal number of rows in second matrix"); std::vector shape(2); shape[0]=getPntrToArgument(0)->getShape()[0]; shape[1]=getPntrToArgument(1)->getShape()[1]; - addValue( shape ); setNotPeriodic(); nderivatives = buildArgumentStore(0); - std::string headstr=getFirstActionInChain()->getLabel(); - stored_matrix1 = getPntrToArgument(0)->ignoreStoredValue( headstr ); - stored_matrix2 = getPntrToArgument(1)->ignoreStoredValue( headstr ); + addValue( shape ); setNotPeriodic(); parseFlag("ELEMENTS_ON_DIAGONAL_ARE_ZERO",diagzero); if( diagzero ) log.printf(" setting diagonal elements equal to zero\n"); @@ -101,7 +95,7 @@ MatrixTimesMatrix::MatrixTimesMatrix(const ActionOptions&ao): } unsigned MatrixTimesMatrix::getNumberOfDerivatives() { - return nderivatives; + return getPntrToArgument(0)->getNumberOfStoredValues() + getPntrToArgument(1)->getNumberOfStoredValues(); } unsigned MatrixTimesMatrix::getNumberOfColumns() const { @@ -113,14 +107,7 @@ void MatrixTimesMatrix::prepare() { ActionWithVector::prepare(); Value* myval = getPntrToComponent(0); if( myval->getShape()[0]==getPntrToArgument(0)->getShape()[0] && myval->getShape()[1]==getPntrToArgument(1)->getShape()[1] ) return; std::vector shape(2); shape[0]=getPntrToArgument(0)->getShape()[0]; shape[1]=getPntrToArgument(1)->getShape()[1]; - myval->setShape(shape); if( myval->valueIsStored() ) myval->reshapeMatrixStore( shape[1] ); -} - -void MatrixTimesMatrix::getAdditionalTasksRequired( ActionWithVector* action, std::vector& atasks ) { - - ActionWithMatrix* adj=dynamic_cast( getPntrToArgument(0)->getPntrToAction() ); - if( !adj->isAdjacencyMatrix() ) return; - adj->retrieveAtoms(); adj->getAdditionalTasksRequired( action, atasks ); + myval->setShape(shape); myval->reshapeMatrixStore( shape[1] ); } void MatrixTimesMatrix::setupForTask( const unsigned& task_index, std::vector& indices, MultiValue& myvals ) const { @@ -155,103 +142,61 @@ void MatrixTimesMatrix::performTask( const std::string& controller, const unsign Value* myarg = getPntrToArgument(0); unsigned nmult=myarg->getRowLength(index1); double matval=0; std::vector dvec1(nmult), dvec2(nmult); std::vector colno(nmult); - if( !stored_matrix1 && !stored_matrix2 ) { - Value* myarg2 = getPntrToArgument(1); - unsigned ncols = myarg->getNumberOfColumns(); - unsigned ncols2 = myarg2->getNumberOfColumns(); - unsigned base = myarg->getNumberOfStoredValues(); - for(unsigned i=0; igetRowIndex( index1, i ); - if( ncols2getShape()[1] ) { - unsigned nr = myarg2->getRowLength(kind); - for(unsigned j=0; jgetRowIndex( kind, j )==ind2 ) { colno[i]=j; break; } - } - if( colno[i]<0 ) continue; - } else colno[i] = ind2; - double val1 = myarg->get( index1*ncols + i, false ); - double val2 = myarg2->get( kind*ncols2 + colno[i], false ); - if( getName()=="DISSIMILARITIES" ) { - double tmp = getPntrToArgument(0)->difference(val2, val1); matval += tmp*tmp; - if( !squared ) { - dvec1[i] = 2*tmp; dvec2[i] = -2*tmp; continue; - } else { val2 = -2*tmp; val1 = 2*tmp; } - } else matval+= val1*val2; - - if( !squared || doNotCalculateDerivatives() ) continue ; - myvals.addDerivative( 0, index1*ncols + i, val2 ); myvals.updateIndex( 0, index1*ncols + i ); - myvals.addDerivative( 0, base + kind*ncols2 + colno[i], val1 ); myvals.updateIndex( 0, base + kind*ncols2 + colno[i] ); - } - // And add this part of the product - if( !squared ) matval = sqrt(matval); - myvals.addValue( 0, matval ); - if( squared || doNotCalculateDerivatives() ) return; - - for(unsigned i=0; igetRowIndex( index1, i ); if( colno[i]<0 ) continue; - myvals.addDerivative( 0, index1*ncols + i, dvec1[i]/(2*matval) ); myvals.updateIndex( 0, index1*ncols + i ); - myvals.addDerivative( 0, base + i*ncols2 + colno[i], dvec2[i]/(2*matval) ); myvals.updateIndex( 0, base + i*ncols2 + colno[i] ); - } - } else { - for(unsigned i=0; igetRowIndex( index1, i ); - double val1 = getElementOfMatrixArgument( 0, index1, kind, myvals ); - double val2 = getElementOfMatrixArgument( 1, kind, ind2, myvals ); - if( getName()=="DISSIMILARITIES" ) { - double tmp = getPntrToArgument(0)->difference(val2, val1); matval += tmp*tmp; - if( !squared ) { - dvec1[i] = 2*tmp; dvec2[i] = -2*tmp; continue; - } else { val2 = -2*tmp; val1 = 2*tmp; } - } else matval+= val1*val2; - - if( doNotCalculateDerivatives() ) continue; - - addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, kind, val2, myvals ); - addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, kind, ind2, val1, myvals ); - } - // And add this part of the product - if( !squared ) matval = sqrt(matval); - myvals.addValue( 0, matval ); - if( squared || doNotCalculateDerivatives() ) return; - - for(unsigned i=0; igetRowIndex( index1, i ); - addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, kind, dvec1[i]/(2*matval), myvals ); - addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, kind, ind2, dvec2[i]/(2*matval), myvals ); - } + Value* myarg2 = getPntrToArgument(1); + unsigned ncols = myarg->getNumberOfColumns(); + unsigned ncols2 = myarg2->getNumberOfColumns(); + unsigned base = myarg->getNumberOfStoredValues(); + for(unsigned i=0; igetRowIndex( index1, i ); + if( ncols2getShape()[1] ) { + unsigned nr = myarg2->getRowLength(kind); + for(unsigned j=0; jgetRowIndex( kind, j )==ind2 ) { colno[i]=j; break; } + } + if( colno[i]<0 ) continue; + } else colno[i] = ind2; + double val1 = myarg->get( index1*ncols + i, false ); + double val2 = myarg2->get( kind*ncols2 + colno[i], false ); + if( getName()=="DISSIMILARITIES" ) { + double tmp = getPntrToArgument(0)->difference(val2, val1); matval += tmp*tmp; + if( !squared ) { + dvec1[i] = 2*tmp; dvec2[i] = -2*tmp; continue; + } else { val2 = -2*tmp; val1 = 2*tmp; } + } else matval+= val1*val2; + + if( !squared || doNotCalculateDerivatives() ) continue ; + myvals.addDerivative( 0, index1*ncols + i, val2 ); myvals.updateIndex( 0, index1*ncols + i ); + myvals.addDerivative( 0, base + kind*ncols2 + colno[i], val1 ); myvals.updateIndex( 0, base + kind*ncols2 + colno[i] ); + } + // And add this part of the product + if( !squared ) matval = sqrt(matval); + myvals.addValue( 0, matval ); + if( squared || doNotCalculateDerivatives() ) return; + + for(unsigned i=0; igetRowIndex( index1, i ); if( colno[i]<0 ) continue; + myvals.addDerivative( 0, index1*ncols + i, dvec1[i]/(2*matval) ); myvals.updateIndex( 0, index1*ncols + i ); + myvals.addDerivative( 0, base + i*ncols2 + colno[i], dvec2[i]/(2*matval) ); myvals.updateIndex( 0, base + i*ncols2 + colno[i] ); } } void MatrixTimesMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector & indices, MultiValue& myvals ) const { if( doNotCalculateDerivatives() ) return ; - unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat ); - std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); - if( !stored_matrix1 && !stored_matrix2 ) { - unsigned mat1s = ival*getPntrToArgument(0)->getNumberOfColumns(); - unsigned nmult = getPntrToArgument(0)->getRowLength(ival); - for(unsigned j=0; jgetShape()[0]; - unsigned base = getPntrToArgument(0)->getNumberOfStoredValues(); - for(unsigned j=0; jgetRowLength(j); - for(unsigned k=0; kgetNumberOfColumns(); - } - } else { - unsigned ntwo_atoms = myvals.getSplitIndex(); - unsigned mat1s = ival*getPntrToArgument(0)->getShape()[1]; - unsigned nmult = getPntrToArgument(0)->getShape()[1], ss = getPntrToArgument(1)->getShape()[1]; - for(unsigned j=0; j=getPntrToArgument(0)->getShape()[0] ) ind2 = indices[i] - getPntrToArgument(0)->getShape()[0]; - matrix_indices[nmat_ind] = arg_deriv_starts[1] + j*ss + ind2; nmat_ind++; - } - } + unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives(); + std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices() ); + unsigned mat1s = ival*getPntrToArgument(0)->getNumberOfColumns(); + unsigned nmult = getPntrToArgument(0)->getRowLength(ival); + for(unsigned j=0; jgetShape()[0]; + unsigned base = getPntrToArgument(0)->getNumberOfStoredValues(); + for(unsigned j=0; jgetRowLength(j); + for(unsigned k=0; kgetNumberOfColumns(); } - myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind ); + myvals.setNumberOfMatrixRowDerivatives( nmat_ind ); } } diff --git a/src/matrixtools/MatrixTimesVector.cpp b/src/matrixtools/MatrixTimesVector.cpp index 6f1f1ef54f..462542a8fa 100644 --- a/src/matrixtools/MatrixTimesVector.cpp +++ b/src/matrixtools/MatrixTimesVector.cpp @@ -41,7 +41,6 @@ class MatrixTimesVector : public ActionWithVector { static void registerKeywords( Keywords& keys ); explicit MatrixTimesVector(const ActionOptions&); std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ; - void getNumberOfStreamedDerivatives( unsigned& nderivatives, Value* stopat ) override ; unsigned getNumberOfDerivatives() override ; void prepare() override ; void calculate() override ; @@ -103,7 +102,6 @@ MatrixTimesVector::MatrixTimesVector(const ActionOptions&ao): if( getPntrToArgument(n)->getRank()>0 ) { if( getPntrToArgument(n)->getRank()!=1 || getPntrToArgument(n)->hasDerivatives() ) error("last argument to this action should be a vector"); } - getPntrToArgument(n)->buildDataStore(); if( getNumberOfArguments()==2 ) { addValue( shape ); setNotPeriodic(); @@ -131,7 +129,6 @@ MatrixTimesVector::MatrixTimesVector(const ActionOptions&ao): if( getPntrToArgument(i)->getRank()==0 ) { if( getPntrToArgument(0)->getShape()[1]!=1 ) error("number of columns in input matrix does not equal number of elements in vector"); } else if( getPntrToArgument(0)->getShape()[1]!=getPntrToArgument(i)->getShape()[0] ) error("number of columns in input matrix does not equal number of elements in vector"); - getPntrToArgument(i)->buildDataStore(); } for(unsigned i=1; i shape(1); shape[0] = getPntrToArgument(0)->getShape()[0]; myval->setShape(shape); } -void MatrixTimesVector::getNumberOfStreamedDerivatives( unsigned& nderivatives, Value* stopat ) { - nderivatives = 0; - for(unsigned i=0; igetNumberOfStoredValues(); - } -} - void MatrixTimesVector::calculate() { runAllTasks(); } @@ -186,7 +173,7 @@ int MatrixTimesVector::checkTaskIsActive( const unsigned& itask ) const { void MatrixTimesVector::performTask( const unsigned& task_index, MultiValue& myvals ) const { if( sumrows ) { - unsigned n=getNumberOfArguments()-1; Value* myvec = getPntrToArgument(n); + unsigned base=0, n=getNumberOfArguments()-1; Value* myvec = getPntrToArgument(n); for(unsigned i=0; igetNumberOfColumns(); @@ -197,18 +184,20 @@ void MatrixTimesVector::performTask( const unsigned& task_index, MultiValue& myv // And the derivatives if( doNotCalculateDerivatives() ) continue; - unsigned dloc = arg_deriv_starts[i] + task_index*ncol; + unsigned dloc = base + task_index*ncol; for(unsigned j=0; jgetNumberOfStoredValues(); } } else if( getPntrToArgument(1)->getRank()==1 ) { Value* mymat = getPntrToArgument(0); + unsigned base = 0; unsigned ncol = mymat->getNumberOfColumns(); unsigned nmat = mymat->getRowLength(task_index); - unsigned dloc = arg_deriv_starts[0] + task_index*ncol; + unsigned dloc = task_index*ncol; for(unsigned i=0; igetNumberOfStoredValues(); double val=0; for(unsigned j=0; jget( task_index*ncol + j, false )*myvec->get( mymat->getRowIndex( task_index, j ) ); myvals.setValue( i, val ); @@ -220,11 +209,12 @@ void MatrixTimesVector::performTask( const unsigned& task_index, MultiValue& myv double vecval = myvec->get( kind ); double matval = mymat->get( task_index*ncol + j, false ); myvals.addDerivative( i, dloc + j, vecval ); myvals.updateIndex( i, dloc + j ); - myvals.addDerivative( i, arg_deriv_starts[i+1] + kind, matval ); myvals.updateIndex( i, arg_deriv_starts[i+1] + kind ); + myvals.addDerivative( i, base + kind, matval ); myvals.updateIndex( i, base + kind ); } } } else { - unsigned n=getNumberOfArguments()-1; Value* myvec = getPntrToArgument(n); + unsigned base=0, n=getNumberOfArguments()-1; Value* myvec = getPntrToArgument(n); + unsigned nmat_der = 0; for(unsigned i=0; igetNumberOfStoredValues(); for(unsigned i=0; igetNumberOfColumns(); @@ -235,14 +225,15 @@ void MatrixTimesVector::performTask( const unsigned& task_index, MultiValue& myv // And the derivatives if( doNotCalculateDerivatives() ) continue; - unsigned dloc = arg_deriv_starts[i] + task_index*ncol; + unsigned dloc = base + task_index*ncol; for(unsigned j=0; jgetRowIndex( task_index, j ); double vecval = myvec->get( kind ); double matval = mymat->get( task_index*ncol + j, false ); myvals.addDerivative( i, dloc + j, vecval ); myvals.updateIndex( i, dloc + j ); - myvals.addDerivative( i, arg_deriv_starts[n] + kind, matval ); myvals.updateIndex( i, arg_deriv_starts[n] + kind ); + myvals.addDerivative( i, nmat_der + kind, matval ); myvals.updateIndex( i, nmat_der + kind ); } + base += mymat->getNumberOfStoredValues(); } } } diff --git a/src/matrixtools/OuterProduct.cpp b/src/matrixtools/OuterProduct.cpp index ffa141b110..e442a645e3 100644 --- a/src/matrixtools/OuterProduct.cpp +++ b/src/matrixtools/OuterProduct.cpp @@ -39,8 +39,6 @@ class OuterProduct : public ActionWithMatrix { private: bool domin, domax, diagzero; LeptonCall function; - unsigned nderivatives; - bool stored_vector1, stored_vector2; public: static void registerKeywords( Keywords& keys ); explicit OuterProduct(const ActionOptions&); @@ -93,15 +91,12 @@ OuterProduct::OuterProduct(const ActionOptions&ao): if( diagzero ) log.printf(" setting diagonal elements equal to zero\n"); std::vector shape(2); shape[0]=getPntrToArgument(0)->getShape()[0]; shape[1]=getPntrToArgument(1)->getShape()[0]; - addValue( shape ); setNotPeriodic(); nderivatives = buildArgumentStore(0); - std::string headstr=getFirstActionInChain()->getLabel(); - stored_vector1 = getPntrToArgument(0)->ignoreStoredValue( headstr ); - stored_vector2 = getPntrToArgument(1)->ignoreStoredValue( headstr ); + addValue( shape ); setNotPeriodic(); if( getPntrToArgument(0)->isDerivativeZeroWhenValueIsZero() || getPntrToArgument(1)->isDerivativeZeroWhenValueIsZero() ) getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero(); } unsigned OuterProduct::getNumberOfDerivatives() { - return nderivatives; + return getPntrToArgument(0)->getNumberOfStoredValues() + getPntrToArgument(1)->getNumberOfStoredValues(); } unsigned OuterProduct::getNumberOfColumns() const { @@ -146,34 +141,37 @@ void OuterProduct::performTask( const std::string& controller, const unsigned& i if( index2>=getPntrToArgument(0)->getShape()[0] ) ind2 = index2 - getPntrToArgument(0)->getShape()[0]; if( diagzero && index1==ind2 ) return; - double fval; unsigned jarg = 0, kelem = index1; bool jstore=stored_vector1; + double fval; unsigned jarg = 0, kelem = index1; std::vector args(2); - args[0] = getArgumentElement( 0, index1, myvals ); - args[1] = getArgumentElement( 1, ind2, myvals ); + args[0] = getPntrToArgument(0)->get( index1 ); + args[1] = getPntrToArgument(1)->get( ind2 ); if( domin ) { - fval=args[0]; if( args[1]getNumberOfStoredValues(); kelem=ind2; } } else if( domax ) { - fval=args[0]; if( args[1]>args[0] ) { fval=args[1]; jarg=1; kelem=ind2; jstore=stored_vector2; } + fval=args[0]; if( args[1]>args[0] ) { fval=args[1]; jarg=getPntrToArgument(0)->getNumberOfStoredValues(); kelem=ind2; } } else { fval=function.evaluate( args ); } myvals.addValue( 0, fval ); if( doNotCalculateDerivatives() ) return ; if( domin || domax ) { - addDerivativeOnVectorArgument( jstore, 0, jarg, kelem, 1.0, myvals ); + myvals.addDerivative( 0, jarg + kelem, 1.0 ); myvals.updateIndex( 0, jarg + kelem ); } else { - addDerivativeOnVectorArgument( stored_vector1, 0, 0, index1, function.evaluateDeriv( 0, args ), myvals ); - addDerivativeOnVectorArgument( stored_vector2, 0, 1, ind2, function.evaluateDeriv( 1, args ), myvals ); + myvals.addDerivative( 0, index1, function.evaluateDeriv( 0, args ) ); myvals.updateIndex( 0, index1 ); + myvals.addDerivative( 0, getPntrToArgument(0)->getNumberOfStoredValues() + ind2, function.evaluateDeriv( 1, args ) ); + myvals.updateIndex( 0, getPntrToArgument(0)->getNumberOfStoredValues() + ind2 ); } if( doNotCalculateDerivatives() ) return ; - unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat ); - myvals.getMatrixRowDerivativeIndices( nmat )[nmat_ind] = arg_deriv_starts[1] + ind2; myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+1 ); + unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives(); + myvals.getMatrixRowDerivativeIndices()[nmat_ind] = getPntrToArgument(0)->getNumberOfStoredValues() + ind2; + myvals.setNumberOfMatrixRowDerivatives( nmat_ind+1 ); } void OuterProduct::runEndOfRowJobs( const unsigned& ival, const std::vector & indices, MultiValue& myvals ) const { if( doNotCalculateDerivatives() ) return ; - unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat ); - myvals.getMatrixRowDerivativeIndices( nmat )[nmat_ind] = ival; myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+1 ); + unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives(); + myvals.getMatrixRowDerivativeIndices()[nmat_ind] = ival; + myvals.setNumberOfMatrixRowDerivatives( nmat_ind+1 ); } } diff --git a/src/matrixtools/TransposeMatrix.cpp b/src/matrixtools/TransposeMatrix.cpp index 125bd43b31..ffc8c32b8a 100644 --- a/src/matrixtools/TransposeMatrix.cpp +++ b/src/matrixtools/TransposeMatrix.cpp @@ -68,7 +68,7 @@ TransposeMatrix::TransposeMatrix(const ActionOptions& ao): else if( getPntrToArgument(0)->getRank()==1 ) { shape.resize(2); shape[0]=1; shape[1]=getPntrToArgument(0)->getShape()[0]; } else if( getPntrToArgument(0)->getShape()[0]==1 ) { shape.resize(1); shape[0] = getPntrToArgument(0)->getShape()[1]; } else { shape.resize(2); shape[0]=getPntrToArgument(0)->getShape()[1]; shape[1]=getPntrToArgument(0)->getShape()[0]; } - addValue( shape ); setNotPeriodic(); getPntrToComponent(0)->buildDataStore(); + addValue( shape ); setNotPeriodic(); if( shape.size()==2 ) getPntrToComponent(0)->reshapeMatrixStore( shape[1] ); } diff --git a/src/matrixtools/Voronoi.cpp b/src/matrixtools/Voronoi.cpp index 29137b8bd8..04dee78f17 100644 --- a/src/matrixtools/Voronoi.cpp +++ b/src/matrixtools/Voronoi.cpp @@ -62,8 +62,7 @@ Voronoi::Voronoi(const ActionOptions&ao): if( getNumberOfArguments()!=1 ) error("should be one arguments to this action, a matrix"); if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) error("argument to this action should be a matrix"); if( getPntrToArgument(0)->getShape()[1]>getPntrToArgument(0)->getShape()[0] ) warning("would expect number of columns in matrix to exceed number of rows"); - getPntrToArgument(0)->buildDataStore(); std::vector shape( getPntrToArgument(0)->getShape() ); - addValue( shape ); setNotPeriodic(); getPntrToComponent(0)->buildDataStore(); + std::vector shape( getPntrToArgument(0)->getShape() ); addValue( shape ); setNotPeriodic(); } void Voronoi::prepare() { diff --git a/src/metatensor/metatensor.cpp b/src/metatensor/metatensor.cpp index 7052d38f08..15a736e770 100644 --- a/src/metatensor/metatensor.cpp +++ b/src/metatensor/metatensor.cpp @@ -556,17 +556,14 @@ MetatensorPlumedAction::MetatensorPlumedAction(const ActionOptions& options): log.printf(" the output of this model is 1x%d vector\n", n_properties_); this->addValue({this->n_properties_}); - this->getPntrToComponent(0)->buildDataStore(); } else if (n_properties_ == 1) { log.printf(" the output of this model is %dx1 vector\n", n_samples_); this->addValue({this->n_samples_}); - this->getPntrToComponent(0)->buildDataStore(); } else { log.printf(" the output of this model is a %dx%d matrix\n", n_samples_, n_properties_); this->addValue({this->n_samples_, this->n_properties_}); - this->getPntrToComponent(0)->buildDataStore(); this->getPntrToComponent(0)->reshapeMatrixStore(n_properties_); } diff --git a/src/refdist/MatrixProductDiagonal.cpp b/src/refdist/MatrixProductDiagonal.cpp index 3c6ff8f0e9..d360121e2b 100644 --- a/src/refdist/MatrixProductDiagonal.cpp +++ b/src/refdist/MatrixProductDiagonal.cpp @@ -78,7 +78,6 @@ MatrixProductDiagonal::MatrixProductDiagonal(const ActionOptions&ao): std::vector shape(1); shape[0]=getPntrToArgument(0)->getShape()[0]; addValue( shape ); setNotPeriodic(); } - getPntrToArgument(0)->buildDataStore(); getPntrToArgument(1)->buildDataStore(); } unsigned MatrixProductDiagonal::getNumberOfDerivatives() { @@ -128,7 +127,7 @@ void MatrixProductDiagonal::performTask( const unsigned& task_index, MultiValue& void MatrixProductDiagonal::calculate() { if( getPntrToArgument(1)->getRank()==1 ) { unsigned nder = getNumberOfDerivatives(); - MultiValue myvals( 1, nder, 0, 0 ); performTask( 0, myvals ); + MultiValue myvals( 1, nder, 0 ); performTask( 0, myvals ); Value* myval=getPntrToComponent(0); myval->set( myvals.get(0) ); for(unsigned i=0; isetDerivative( i, myvals.getDerivative(0,i) ); diff --git a/src/secondarystructure/SecondaryStructureRMSD.cpp b/src/secondarystructure/SecondaryStructureRMSD.cpp index dc2124a3d4..2ab57b3fef 100644 --- a/src/secondarystructure/SecondaryStructureRMSD.cpp +++ b/src/secondarystructure/SecondaryStructureRMSD.cpp @@ -190,7 +190,7 @@ void SecondaryStructureRMSD::calculate() { void SecondaryStructureRMSD::performTask( const unsigned& current, MultiValue& myvals ) const { // Resize the derivatives if need be unsigned nderi = 3*getNumberOfAtoms()+9; - if( myvals.getNumberOfDerivatives()!=nderi ) myvals.resize( myvals.getNumberOfValues(), nderi, 0, 0 ); + if( myvals.getNumberOfDerivatives()!=nderi ) myvals.resize( myvals.getNumberOfValues(), nderi, 0 ); // Retrieve the positions const unsigned natoms = colvar_atoms[current].size(); std::vector pos( natoms ), deriv( natoms ); diff --git a/src/symfunc/ThreeBodyGFunctions.cpp b/src/symfunc/ThreeBodyGFunctions.cpp index fdc0521eda..3890fbabea 100644 --- a/src/symfunc/ThreeBodyGFunctions.cpp +++ b/src/symfunc/ThreeBodyGFunctions.cpp @@ -73,7 +73,6 @@ ThreeBodyGFunctions::ThreeBodyGFunctions(const ActionOptions&ao): log.printf(" using bond weights from matrix labelled %s \n",wval[0]->getName().c_str() ); // Rerequest the arguments std::vector myargs( getArguments() ); myargs.push_back( wval[0] ); requestArguments( myargs ); - for(unsigned i=0; ibuildDataStore(); std::vector shape(1); shape[0] = getPntrToArgument(0)->getShape()[0]; // And now read the functions to compute diff --git a/src/tools/MultiValue.cpp b/src/tools/MultiValue.cpp index fd26ccea18..f56950ad87 100644 --- a/src/tools/MultiValue.cpp +++ b/src/tools/MultiValue.cpp @@ -24,7 +24,7 @@ namespace PLMD { -MultiValue::MultiValue( const size_t& nvals, const size_t& nder, const size_t& nmat, const size_t& maxcol ): +MultiValue::MultiValue( const size_t& nvals, const size_t& nder, const size_t& nmat ): task_index(0), task2_index(0), values(nvals), @@ -37,31 +37,26 @@ MultiValue::MultiValue( const size_t& nvals, const size_t& nder, const size_t& n vector_call(false), nindices(0), nsplit(0), - nmatrix_cols(maxcol), - matrix_force_stash(nder*nmat), - matrix_row_nderivatives(nmat,0), - matrix_row_derivative_indices(nmat) + matrix_force_stash(0), + matrix_row_nderivatives(0), + matrix_row_derivative_indices(nder) { - for(unsigned i=0; i0 ) matrix_force_stash.resize(nder,0); } -void MultiValue::resize( const size_t& nvals, const size_t& nder, const size_t& nmat, const size_t& maxcol ) { - if( values.size()==nvals && nderivatives==nder && matrix_row_nderivatives.size()==nmat && nmatrix_cols==maxcol ) return; +void MultiValue::resize( const size_t& nvals, const size_t& nder, const size_t& nmat ) { + if( values.size()==nvals && nderivatives==nder ) return; values.resize(nvals); nderivatives=nder; derivatives.resize( nvals*nder ); hasderiv.resize(nvals*nder,false); nactive.resize(nvals); active_list.resize(nvals*nder); - nmatrix_cols=maxcol; matrix_force_stash.resize(nmat*nder,0); - matrix_row_nderivatives.resize(nmat,0); matrix_row_derivative_indices.resize(nmat); atLeastOneSet=false; - for(unsigned i=0; i0 ) matrix_force_stash.resize(nder,0); + matrix_row_nderivatives=0; matrix_row_derivative_indices.resize(nder); atLeastOneSet=false; } void MultiValue::clearAll() { for(unsigned i=0; i matrix_force_stash; /// These are used to store the indices that have derivatives wrt to at least one /// of the elements in a matrix - std::vector matrix_row_nderivatives; - std::vector > matrix_row_derivative_indices; + unsigned matrix_row_nderivatives; + std::vector matrix_row_derivative_indices; /// This is a fudge to save on vector resizing in MultiColvar std::vector indices; std::vector tmp_atoms; @@ -64,8 +62,8 @@ class MultiValue { std::vector tmp_atom_virial; std::vector > tmp_vectors; public: - MultiValue( const std::size_t& nvals, const std::size_t& nder, const std::size_t& nmat=0, const std::size_t& maxcol=0 ); - void resize( const std::size_t& nvals, const std::size_t& nder, const std::size_t& nmat=0, const std::size_t& maxcol=0 ); + MultiValue( const std::size_t& nvals, const std::size_t& nder, const std::size_t& nmat=0 ); + void resize( const std::size_t& nvals, const std::size_t& nder, const std::size_t& nmat=0 ); /// Set the task index prior to the loop void setTaskIndex( const std::size_t& tindex ); /// @@ -131,13 +129,13 @@ class MultiValue { /// unsigned getActiveIndex( const unsigned& ) const ; /// Get the bookeeping stuff for the derivatives wrt to rows of matrix - void setNumberOfMatrixRowDerivatives( const unsigned& nmat, const unsigned& nind ); - unsigned getNumberOfMatrixRowDerivatives( const unsigned& nmat ) const ; - std::vector& getMatrixRowDerivativeIndices( const unsigned& nmat ); - const std::vector& getMatrixRowDerivativeIndices( const unsigned& nmat ) const ; + void setNumberOfMatrixRowDerivatives( const unsigned& nind ); + unsigned getNumberOfMatrixRowDerivatives() const ; + std::vector& getMatrixRowDerivativeIndices(); + const std::vector& getMatrixRowDerivativeIndices() const ; /// Stash the forces on the matrix - void addMatrixForce( const unsigned& imat, const unsigned& jind, const double& f ); - double getStashedMatrixForce( const unsigned& imat, const unsigned& jind ) const ; + void addMatrixForce( const unsigned& jind, const double& f ); + double getStashedMatrixForce( const unsigned& jind ) const ; }; inline @@ -235,7 +233,7 @@ std::size_t MultiValue::getNumberOfIndices() const { inline bool MultiValue::inVectorCall() const { - return (matrix_row_nderivatives.size()>0 && vector_call); + return (matrix_force_stash.size()>0 && vector_call); } inline @@ -294,34 +292,34 @@ std::vector& MultiValue::getFirstAtomVirialVector() { } inline -void MultiValue::setNumberOfMatrixRowDerivatives( const unsigned& nmat, const unsigned& nind ) { - plumed_dbg_assert( nmat& MultiValue::getMatrixRowDerivativeIndices( const unsigned& nmat ) { - plumed_dbg_assert( nmat& MultiValue::getMatrixRowDerivativeIndices() { + return matrix_row_derivative_indices; } inline -const std::vector& MultiValue::getMatrixRowDerivativeIndices( const unsigned& nmat ) const { - plumed_dbg_assert( nmat& MultiValue::getMatrixRowDerivativeIndices() const { + return matrix_row_derivative_indices; } inline -void MultiValue::addMatrixForce( const unsigned& imat, const unsigned& jind, const double& f ) { - matrix_force_stash[imat*nderivatives + jind]+=f; +void MultiValue::addMatrixForce( const unsigned& jind, const double& f ) { + matrix_force_stash[jind]+=f; } inline -double MultiValue::getStashedMatrixForce( const unsigned& imat, const unsigned& jind ) const { - return matrix_force_stash[imat*nderivatives + jind]; +double MultiValue::getStashedMatrixForce( const unsigned& jind ) const { + return matrix_force_stash[jind]; } inline diff --git a/src/valtools/Concatenate.cpp b/src/valtools/Concatenate.cpp index 8235762fdc..117e498275 100644 --- a/src/valtools/Concatenate.cpp +++ b/src/valtools/Concatenate.cpp @@ -71,7 +71,7 @@ Concatenate::Concatenate(const ActionOptions& ao): vectors=true; std::vector shape(1); shape[0]=0; for(unsigned i=0; igetRank()>1 ) error("cannot concatenate matrix with vectors"); - getPntrToArgument(i)->buildDataStore(); shape[0] += getPntrToArgument(i)->getNumberOfValues(); + shape[0] += getPntrToArgument(i)->getNumberOfValues(); } log.printf(" creating vector with %d elements \n", shape[0] ); addValue( shape ); bool period=getPntrToArgument(0)->isPeriodic(); @@ -84,7 +84,6 @@ Concatenate::Concatenate(const ActionOptions& ao): } } if( period ) setPeriodic( min, max ); else setNotPeriodic(); - getPntrToComponent(0)->buildDataStore(); if( getPntrToComponent(0)->getRank()==2 ) getPntrToComponent(0)->reshapeMatrixStore( shape[1] ); } else { unsigned nrows=0, ncols=0; std::vector arglist; vectors=false; @@ -96,7 +95,7 @@ Concatenate::Concatenate(const ActionOptions& ao): if( argn.size()==0 ) break; if( argn.size()>1 ) error("should only be one argument to each matrix keyword"); if( argn[0]->getRank()!=0 && argn[0]->getRank()!=2 ) error("input arguments for this action should be matrices"); - argn[0]->buildDataStore(); arglist.push_back( argn[0] ); nt_cols++; + arglist.push_back( argn[0] ); nt_cols++; if( argn[0]->getRank()==0 ) log.printf(" %d %d component of composed matrix is scalar labelled %s\n", i, j, argn[0]->getName().c_str() ); else log.printf(" %d %d component of composed matrix is %d by %d matrix labelled %s\n", i, j, argn[0]->getShape()[0], argn[0]->getShape()[1], argn[0]->getName().c_str() ); } @@ -125,7 +124,7 @@ Concatenate::Concatenate(const ActionOptions& ao): else shape[0] += arglist[k-1]->getShape()[0]; } // Now request the arguments to make sure we store things we need - requestArguments(arglist); addValue( shape ); setNotPeriodic(); getPntrToComponent(0)->buildDataStore(); + requestArguments(arglist); addValue( shape ); setNotPeriodic(); if( getPntrToComponent(0)->getRank()==2 ) getPntrToComponent(0)->reshapeMatrixStore( shape[1] ); } } diff --git a/src/valtools/Flatten.cpp b/src/valtools/Flatten.cpp index e1f8e52ea8..93f49cbf0c 100644 --- a/src/valtools/Flatten.cpp +++ b/src/valtools/Flatten.cpp @@ -67,10 +67,9 @@ Flatten::Flatten(const ActionOptions& ao): { if( getNumberOfArguments()!=1 ) error("should only be one argument for this action"); if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) error("input to this action should be a matrix"); - getPntrToArgument(0)->buildDataStore(true); std::vector inshape( getPntrToArgument(0)->getShape() ); std::vector shape( 1 ); shape[0]=inshape[0]*inshape[1]; - addValue( shape ); setNotPeriodic(); getPntrToComponent(0)->buildDataStore(); + addValue( shape ); setNotPeriodic(); } void Flatten::calculate() { diff --git a/src/valtools/SelectWithMask.cpp b/src/valtools/SelectWithMask.cpp index ac480f15b7..fffba3b080 100644 --- a/src/valtools/SelectWithMask.cpp +++ b/src/valtools/SelectWithMask.cpp @@ -75,7 +75,7 @@ SelectWithMask::SelectWithMask(const ActionOptions& ao): ActionWithArguments(ao) { if( getNumberOfArguments()!=1 ) error("should only be one argument for this action"); - getPntrToArgument(0)->buildDataStore(); std::vector shape; + std::vector shape; if( getPntrToArgument(0)->getRank()==1 ) { std::vector mask; parseArgumentList("MASK",mask); if( mask.size()!=1 ) error("should only be one input for mask"); @@ -97,13 +97,13 @@ SelectWithMask::SelectWithMask(const ActionOptions& ao): ActionWithArguments::interpretArgumentList( labs, plumed.getActionSet(), this, rmask ); } shape.resize(2); - rmask[0]->buildDataStore(); shape[0] = getOutputVectorLength( rmask[0] ); - cmask[0]->buildDataStore(); shape[1] = getOutputVectorLength( cmask[0] ); + shape[0] = getOutputVectorLength( rmask[0] ); + shape[1] = getOutputVectorLength( cmask[0] ); std::vector args( getArguments() ); args.push_back( rmask[0] ); args.push_back( cmask[0] ); requestArguments( args ); } else error("input should be vector or matrix"); - addValue( shape ); getPntrToComponent(0)->buildDataStore(); + addValue( shape ); if( getPntrToArgument(0)->isPeriodic() ) { std::string min, max; getPntrToArgument(0)->getDomain( min, max ); setPeriodic( min, max ); } else setNotPeriodic(); diff --git a/src/valtools/VStack.cpp b/src/valtools/VStack.cpp index 17f632c4ce..9e965d97fb 100644 --- a/src/valtools/VStack.cpp +++ b/src/valtools/VStack.cpp @@ -35,8 +35,6 @@ namespace PLMD { namespace valtools { class VStack : public ActionWithMatrix { -private: - std::vector stored; public: static void registerKeywords( Keywords& keys ); /// Constructor @@ -57,8 +55,6 @@ class VStack : public ActionWithMatrix { void runEndOfRowJobs( const unsigned& ival, const std::vector & indices, MultiValue& myvals ) const override ; /// void getMatrixColumnTitles( std::vector& argnames ) const override ; -/// - void gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const override ; }; PLUMED_REGISTER_ACTION(VStack,"VSTACK") @@ -93,21 +89,7 @@ VStack::VStack(const ActionOptions& ao): std::vector shape(2); shape[0]=nvals; shape[1]=getNumberOfArguments(); addValue( shape ); if( periodic ) setPeriodic( smin, smax ); else setNotPeriodic(); // And store this value - getPntrToComponent(0)->buildDataStore(); getPntrToComponent(0)->reshapeMatrixStore( shape[1] ); - // Setup everything so we can build the store - ActionWithVector* av=dynamic_cast( getPntrToArgument(0)->getPntrToAction() ); - if( av ) { - const ActionWithVector* head0 = av->getFirstActionInChain(); - for(unsigned i=0; i( getPntrToArgument(i)->getPntrToAction() ); - if( !avv ) continue; - if( head0!=avv->getFirstActionInChain() ) { break; } - } - } - unsigned nder = buildArgumentStore(0); - // This checks which values have been stored - stored.resize( getNumberOfArguments() ); std::string headstr=getFirstActionInChain()->getLabel(); - for(unsigned i=0; iignoreStoredValue( headstr ); + getPntrToComponent(0)->reshapeMatrixStore( shape[1] ); } void VStack::getMatrixColumnTitles( std::vector& argnames ) const { @@ -139,26 +121,22 @@ int VStack::checkTaskIsActive( const unsigned& itask ) const { void VStack::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const { unsigned ind2 = index2; if( index2>=getConstPntrToComponent(0)->getShape()[0] ) ind2 = index2 - getConstPntrToComponent(0)->getShape()[0]; - myvals.addValue( 0, getArgumentElement( ind2, index1, myvals ) ); + myvals.addValue( 0, getPntrToArgument(ind2)->get( index1 ) ); if( doNotCalculateDerivatives() ) return; - addDerivativeOnVectorArgument( stored[ind2], 0, ind2, index1, 1.0, myvals ); -} - -void VStack::gatherForcesOnStoredValue( const unsigned& ival, const unsigned& itask, const MultiValue& myvals, std::vector& forces ) const { - unsigned matind = getConstPntrToComponent(ival)->getPositionInMatrixStash(); - for(unsigned i=0; igetNumberOfStoredValues(); + myvals.addDerivative( 0, vstart + index1, 1.0 ); myvals.updateIndex( 0, vstart + index1 ); } void VStack::runEndOfRowJobs( const unsigned& ival, const std::vector & indices, MultiValue& myvals ) const { if( doNotCalculateDerivatives() ) return ; - unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat ); - std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); + unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives(); + std::vector& matrix_indices( myvals.getMatrixRowDerivativeIndices() ); plumed_assert( nmat_indgetShape()[0]; for(unsigned i=0; isetDerivativeIsZeroWhenValueIsZero(); } -bool ActionVolume::isInSubChain( unsigned& nder ) { - nder = 0; getFirstActionInChain()->getNumberOfStreamedDerivatives( nder, getPntrToComponent(0) ); - nder = nder - getNumberOfDerivatives(); - return true; -} - void ActionVolume::requestAtoms( const std::vector & a ) { std::vector all_atoms( getAbsoluteIndexes() ); for(unsigned i=0; igetRank()==0 ) { setupRegions(); unsigned nref = getNumberOfAtoms() - 1; Vector wdf; Tensor vir; std::vector refders( nref ); diff --git a/src/volumes/ActionVolume.h b/src/volumes/ActionVolume.h index cafd04ff83..048b2007ea 100644 --- a/src/volumes/ActionVolume.h +++ b/src/volumes/ActionVolume.h @@ -57,7 +57,6 @@ class ActionVolume : public ActionWithVector { int checkTaskStatus( const unsigned& taskno, int& flag ) const override; void calculate(); virtual void setupRegions() = 0; - bool isInSubChain( unsigned& nder ) override ; void performTask( const unsigned&, MultiValue& ) const ; virtual double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector& refders ) const=0; };