Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "ActionWithMatrix.h"
      23             : #include "tools/Communicator.h"
      24             : 
      25             : namespace PLMD {
      26             : 
      27        3176 : void ActionWithMatrix::registerKeywords( Keywords& keys ) {
      28        3176 :   ActionWithVector::registerKeywords( keys ); keys.use("ARG");
      29        3176 : }
      30             : 
      31        1472 : ActionWithMatrix::ActionWithMatrix(const ActionOptions&ao):
      32             :   Action(ao),
      33             :   ActionWithVector(ao),
      34        1472 :   next_action_in_chain(NULL),
      35        1472 :   matrix_to_do_before(NULL),
      36        1472 :   matrix_to_do_after(NULL),
      37        1472 :   clearOnEachCycle(true)
      38             : {
      39        1472 : }
      40             : 
      41        1472 : ActionWithMatrix::~ActionWithMatrix() {
      42        1472 :   if( matrix_to_do_before ) { matrix_to_do_before->matrix_to_do_after=NULL; matrix_to_do_before->next_action_in_chain=NULL; }
      43        1472 : }
      44             : 
      45          16 : void ActionWithMatrix::getAllActionLabelsInMatrixChain( std::vector<std::string>& mylabels ) const {
      46             :   bool found=false;
      47          24 :   for(unsigned i=0; i<mylabels.size(); ++i) {
      48           8 :     if( getLabel()==mylabels[i] ) { found=true; }
      49             :   }
      50          16 :   if( !found ) mylabels.push_back( getLabel() );
      51          16 :   if( matrix_to_do_after ) matrix_to_do_after->getAllActionLabelsInMatrixChain( mylabels );
      52          16 : }
      53             : 
      54       55934 : void ActionWithMatrix::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
      55       55934 :   ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
      56             : 
      57      121740 :   for(int i=0; i<getNumberOfComponents(); ++i) {
      58       65806 :     Value* myval=getPntrToComponent(i);
      59       65806 :     if( myval->getRank()!=2 || myval->hasDerivatives() ) continue;
      60       35389 :     myval->setPositionInMatrixStash(nmat); nmat++;
      61       35389 :     if( !myval->valueIsStored() ) continue;
      62        7794 :     if( myval->getShape()[1]>maxcol ) maxcol=myval->getShape()[1];
      63             :     myval->setMatrixBookeepingStart(nbookeeping);
      64        7794 :     nbookeeping += myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
      65             :   }
      66             :   // Turn off clearning of derivatives after each matrix run if there are no matrices in the output of this action
      67       55934 :   clearOnEachCycle = false;
      68       86351 :   for(int i=0; i<getNumberOfComponents(); ++i) {
      69       60992 :     const Value* myval=getConstPntrToComponent(i);
      70       60992 :     if( myval->getRank()==2 && !myval->hasDerivatives() ) { clearOnEachCycle = true; break; }
      71             :   }
      72             :   // Turn off clearing of derivatives if we have only the values of adjacency matrices
      73       55934 :   if( doNotCalculateDerivatives() && isAdjacencyMatrix() ) clearOnEachCycle = false;
      74       55934 : }
      75             : 
      76        3317 : void ActionWithMatrix::finishChainBuild( ActionWithVector* act ) {
      77        3317 :   ActionWithMatrix* am=dynamic_cast<ActionWithMatrix*>(act); if( !am || act==this ) return;
      78             :   // Build the list that contains everything we are going to loop over in getTotalMatrixBookeepgin and updateAllNeighbourLists
      79        2269 :   if( next_action_in_chain ) next_action_in_chain->finishChainBuild( act );
      80             :   else {
      81         887 :     next_action_in_chain=am;
      82             :     // Build the list of things we are going to loop over in runTask
      83         887 :     if( am->isAdjacencyMatrix() || act->getName()=="VSTACK" ) return ;
      84         828 :     plumed_massert( !matrix_to_do_after, "cannot add " + act->getLabel() + " in " + getLabel() + " as have already added " + matrix_to_do_after->getLabel() );
      85         828 :     matrix_to_do_after=am; am->matrix_to_do_before=this;
      86             :   }
      87             : }
      88             : 
      89        1263 : const ActionWithMatrix* ActionWithMatrix::getFirstMatrixInChain() const {
      90        1263 :   if( !actionInChain() ) return this;
      91         410 :   return matrix_to_do_before->getFirstMatrixInChain();
      92             : }
      93             : 
      94       30968 : void ActionWithMatrix::getTotalMatrixBookeeping( unsigned& nbookeeping ) {
      95       67917 :   for(int i=0; i<getNumberOfComponents(); ++i) {
      96       36949 :     Value* myval=getPntrToComponent(i);
      97       36949 :     if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() ) continue;
      98        4781 :     myval->reshapeMatrixStore( getNumberOfColumns() );
      99        4781 :     nbookeeping += myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
     100             :   }
     101       30968 :   if( next_action_in_chain ) next_action_in_chain->getTotalMatrixBookeeping( nbookeeping );
     102       30968 : }
     103             : 
     104       31010 : void ActionWithMatrix::calculate() {
     105       31010 :   if( actionInChain() ) return ;
     106             :   // Update all the neighbour lists
     107       17481 :   updateAllNeighbourLists();
     108             :   // Setup the matrix indices
     109       17481 :   unsigned nbookeeping=0; getTotalMatrixBookeeping( nbookeeping );
     110       17481 :   if( matrix_bookeeping.size()!=nbookeeping ) matrix_bookeeping.resize( nbookeeping );
     111             :   std::fill( matrix_bookeeping.begin(), matrix_bookeeping.end(), 0 );
     112             :   // And run all the tasks
     113       17481 :   runAllTasks();
     114             : }
     115             : 
     116       30968 : void ActionWithMatrix::updateAllNeighbourLists() {
     117       30968 :   updateNeighbourList();
     118       30968 :   if( next_action_in_chain ) next_action_in_chain->updateAllNeighbourLists();
     119       30968 : }
     120             : 
     121     1379763 : void ActionWithMatrix::performTask( const unsigned& task_index, MultiValue& myvals ) const {
     122             :   std::vector<unsigned> & indices( myvals.getIndices() );
     123     1379763 :   if( matrix_to_do_before ) {
     124             :     plumed_dbg_assert( myvals.inVectorCall() );
     125      714734 :     runEndOfRowJobs( task_index, indices, myvals );
     126      714734 :     return;
     127             :   }
     128      665029 :   setupForTask( task_index, indices, myvals );
     129             : 
     130             :   // Now loop over the row of the matrix
     131             :   unsigned ntwo_atoms = myvals.getSplitIndex();
     132    29790804 :   for(unsigned i=1; i<ntwo_atoms; ++i) {
     133             :     // This does everything in the stream that is done with single matrix elements
     134    29125775 :     runTask( getLabel(), task_index, indices[i], myvals );
     135             :     // Now clear only elements that are not accumulated over whole row
     136    29125775 :     clearMatrixElements( myvals );
     137             :   }
     138             :   // This updates the jobs that need to be completed when we get to the end of a row of the matrix
     139      665029 :   runEndOfRowJobs( task_index, indices, myvals );
     140             : }
     141             : 
     142    88130811 : void ActionWithMatrix::runTask( const std::string& controller, const unsigned& current, const unsigned colno, MultiValue& myvals ) const {
     143    88130811 :   double outval=0; myvals.setTaskIndex(current); myvals.setSecondTaskIndex( colno );
     144    88130811 :   if( isActive() ) performTask( controller, current, colno, myvals );
     145    88130811 :   bool hasval = !isAdjacencyMatrix();
     146   172328138 :   for(int i=0; i<getNumberOfComponents(); ++i) {
     147   145388947 :     if( fabs(myvals.get( getConstPntrToComponent(i)->getPositionInStream()) )>0 ) { hasval=true; break; }
     148             :   }
     149             : 
     150    88130811 :   if( hasval ) {
     151   284622928 :     for(int i=0; i<getNumberOfComponents(); ++i) {
     152   201333077 :       const Value* myval=getConstPntrToComponent(i);
     153   201333077 :       if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() ) continue;
     154    14169812 :       unsigned matindex = myval->getPositionInMatrixStash(), matbook_start = myval->getMatrixBookeepingStart(), col_stash_index = colno;
     155    14169812 :       if( colno>=myval->getShape()[0] ) col_stash_index = colno - myval->getShape()[0];
     156    14169812 :       unsigned rowstart = matbook_start+current*(1+myval->getNumberOfColumns());
     157    14169812 :       if( myval->forcesWereAdded() ) {
     158     1551066 :         unsigned sind = myval->getPositionInStream(), find = myvals.getMatrixBookeeping()[rowstart];
     159     1551066 :         double fforce = myval->getForce( myvals.getTaskIndex()*myval->getNumberOfColumns() + find );
     160     1551066 :         if( getNumberOfColumns()>=myval->getShape()[1] ) fforce = myval->getForce( myvals.getTaskIndex()*myval->getShape()[1] + col_stash_index );
     161    20795223 :         for(unsigned j=0; j<myvals.getNumberActive(sind); ++j) {
     162    19244157 :           unsigned kindex = myvals.getActiveIndex(sind,j); myvals.addMatrixForce( matindex, kindex, fforce*myvals.getDerivative(sind,kindex ) );
     163             :         }
     164             :       }
     165    14169812 :       double finalval = myvals.get( myval->getPositionInStream() );
     166    14169812 :       if( fabs(finalval)>0 ) myvals.stashMatrixElement( matindex, rowstart, col_stash_index, finalval );
     167             :     }
     168             :   }
     169    88130811 :   if( matrix_to_do_after ) matrix_to_do_after->runTask( controller, current, colno, myvals );
     170    88130811 : }
     171             : 
     172       19981 : void ActionWithMatrix::gatherThreads( const unsigned& nt, const unsigned& bufsize, const std::vector<double>& omp_buffer, std::vector<double>& buffer, MultiValue& myvals ) {
     173       19981 :   ActionWithVector::gatherThreads( nt, bufsize, omp_buffer, buffer, myvals );
     174    64708661 :   for(unsigned i=0; i<matrix_bookeeping.size(); ++i) matrix_bookeeping[i] += myvals.getMatrixBookeeping()[i];
     175       19981 : }
     176             : 
     177       17481 : void ActionWithMatrix::gatherProcesses( std::vector<double>& buffer ) {
     178       17481 :   ActionWithVector::gatherProcesses( buffer );
     179       17481 :   if( matrix_bookeeping.size()>0 && !runInSerial() ) comm.Sum( matrix_bookeeping );
     180       17481 :   unsigned nval=0; transferNonZeroMatrixElementsToValues( nval, matrix_bookeeping );
     181       17481 : }
     182             : 
     183       30968 : void ActionWithMatrix::transferNonZeroMatrixElementsToValues( unsigned& nval, const std::vector<unsigned>& matbook ) {
     184       67917 :   for(int i=0; i<getNumberOfComponents(); ++i) {
     185       36949 :     Value* myval=getPntrToComponent(i);
     186       36949 :     if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() || getNumberOfColumns()>=myval->getShape()[1] ) continue;
     187          99 :     unsigned nelements = myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
     188     5060331 :     for(unsigned j=0; j<nelements; ++j) myval->setMatrixBookeepingElement( j, matbook[nval+j] );
     189          99 :     nval += nelements;
     190             :   }
     191       30968 :   if( next_action_in_chain ) next_action_in_chain->transferNonZeroMatrixElementsToValues( nval, matbook );
     192       30968 : }
     193             : 
     194      389280 : void ActionWithMatrix::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
     195             :     const unsigned& bufstart, std::vector<double>& buffer ) const {
     196      389280 :   if( getConstPntrToComponent(valindex)->getRank()==1 ) { ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer ); return; }
     197      229528 :   const Value* myval=getConstPntrToComponent(valindex);
     198             :   unsigned ncols = myval->getNumberOfColumns(), matind = myval->getPositionInMatrixStash();
     199      229528 :   unsigned matbook_start = myval->getMatrixBookeepingStart(), vindex = bufstart + code*myval->getNumberOfColumns();
     200      229528 :   const std::vector<unsigned> & matbook( myvals.getMatrixBookeeping() ); unsigned nelements = matbook[matbook_start+code*(1+ncols)];
     201      229528 :   if( ncols>=myval->getShape()[1] ) {
     202             :     // In this case we store the full matrix
     203     6470140 :     for(unsigned j=0; j<nelements; ++j) {
     204     6275153 :       unsigned jind = matbook[matbook_start+code*(1+ncols)+1+j];
     205             :       plumed_dbg_massert( vindex+j<buffer.size(), "failing in " + getLabel() + " on value " + myval->getName() );
     206     6275153 :       buffer[vindex + jind] += myvals.getStashedMatrixElement( matind, jind );
     207             :     }
     208             :   } else {
     209             :     // This is for storing sparse matrices when we can
     210      965403 :     for(unsigned j=0; j<nelements; ++j) {
     211      930862 :       unsigned jind = matbook[matbook_start+code*(1+ncols)+1+j];
     212             :       plumed_dbg_massert( vindex+j<buffer.size(), "failing in " + getLabel() + " on value " + myval->getName() );
     213      930862 :       buffer[vindex + j] += myvals.getStashedMatrixElement( matind, jind );
     214             :     }
     215             :   }
     216             : }
     217             : 
     218      721141 : bool ActionWithMatrix::checkForTaskForce( const unsigned& itask, const Value* myval ) const {
     219      721141 :   if( myval->getRank()<2 ) return ActionWithVector::checkForTaskForce( itask, myval );
     220      479889 :   unsigned nelements = myval->getRowLength(itask), startr = itask*myval->getNumberOfColumns();
     221    10439083 :   for(unsigned j=0; j<nelements; ++j ) {
     222    10051254 :     if( fabs( myval->getForce( startr + j ) )>epsilon ) return true;
     223             :   }
     224             :   return false;
     225             : }
     226             : 
     227      210505 : void ActionWithMatrix::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
     228      210505 :   if( myval->getRank()==1 ) { ActionWithVector::gatherForcesOnStoredValue( myval, itask, myvals, forces ); return; }
     229             :   unsigned matind = myval->getPositionInMatrixStash();
     230   842701414 :   for(unsigned j=0; j<forces.size(); ++j) forces[j] += myvals.getStashedMatrixForce( matind, j );
     231             : }
     232             : 
     233    88130811 : void ActionWithMatrix::clearMatrixElements( MultiValue& myvals ) const {
     234    88130811 :   if( isActive() && clearOnEachCycle ) {
     235   153880620 :     for(int i=0; i<getNumberOfComponents(); ++i) {
     236   103870566 :       const Value* myval=getConstPntrToComponent(i);
     237   103870566 :       if( myval->getRank()==2 && !myval->hasDerivatives() ) myvals.clearDerivatives( myval->getPositionInStream() );
     238             :     }
     239             :   }
     240    88130811 :   if( matrix_to_do_after ) matrix_to_do_after->clearMatrixElements( myvals );
     241    88130811 : }
     242             : 
     243             : }