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        1194 : void ActionWithMatrix::registerKeywords( Keywords& keys ) {
      28        1194 :   ActionWithVector::registerKeywords( keys );
      29        1194 : }
      30             : 
      31         428 : ActionWithMatrix::ActionWithMatrix(const ActionOptions&ao):
      32             :   Action(ao),
      33             :   ActionWithVector(ao),
      34         428 :   diagzero(false) {
      35             : 
      36         428 :   if( keywords.exists("ELEMENTS_ON_DIAGONAL_ARE_ZERO") ) {
      37         127 :     parseFlag("ELEMENTS_ON_DIAGONAL_ARE_ZERO",diagzero);
      38         127 :     if( diagzero ) {
      39           4 :       log.printf("  setting diagonal elements equal to zero\n");
      40             :     }
      41             :   }
      42         428 : }
      43             : 
      44             : 
      45             : class RequiredMatrixElementsUpdater {
      46             :   RequiredMatrixElements& outmat;
      47             : public:
      48             :   RequiredMatrixElementsUpdater( RequiredMatrixElements& mat ) : outmat(mat) {}
      49             :   ~RequiredMatrixElementsUpdater() {
      50             :     outmat.update();
      51             :   }
      52             : };
      53             : 
      54        1003 : void ActionWithMatrix::updateBookeepingArrays( RequiredMatrixElements& outmat ) {
      55             :   RequiredMatrixElementsUpdater updater(outmat);
      56        1003 :   Value* myval = getPntrToComponent(0);
      57        1003 :   unsigned lstart = myval->getShape()[0];
      58        1003 :   if( getNumberOfMasks()>0 ) {
      59         168 :     Value* maskarg = getPntrToArgument(getNumberOfArguments()-getNumberOfMasks());
      60         390 :     for(unsigned i=0; i<getNumberOfComponents(); ++i) {
      61         222 :       getPntrToComponent(i)->reshapeMatrixStore( maskarg->getNumberOfColumns() );
      62             :     }
      63      262338 :     for(unsigned i=0; i<lstart; ++i) {
      64      262170 :       unsigned rstart = i*(1+myval->getNumberOfColumns());
      65      262170 :       myval->setMatrixBookeepingElement( rstart, maskarg->getRowLength(i) );
      66    16996161 :       for(unsigned j=0; j<maskarg->getRowLength(i); ++j) {
      67    16471821 :         myval->setMatrixBookeepingElement( rstart + 1 + j, maskarg->getRowIndex(i, j) );
      68             :       }
      69             :     }
      70         835 :   } else if ( diagzero ) {
      71          40 :     for(unsigned i=0; i<getNumberOfComponents(); ++i) {
      72          20 :       getPntrToComponent(i)->reshapeMatrixStore( myval->getShape()[1]-1 );
      73             :     }
      74        2000 :     for(unsigned i=0; i<lstart; ++i) {
      75        1980 :       unsigned k=0, rstart = i*(1+myval->getNumberOfColumns());
      76        1980 :       myval->setMatrixBookeepingElement( rstart, myval->getShape()[1]-1 );
      77      198000 :       for(unsigned j=0; j<myval->getShape()[1]; ++j) {
      78      196020 :         if( i!=j ) {
      79      194040 :           myval->setMatrixBookeepingElement( rstart + 1 + k, j );
      80      194040 :           k++;
      81             :         }
      82             :       }
      83             :     }
      84             :   } else {
      85        1648 :     for(unsigned i=0; i<getNumberOfComponents(); ++i) {
      86         833 :       getPntrToComponent(i)->reshapeMatrixStore( myval->getShape()[1] );
      87             :     }
      88             :   }
      89        1003 :   Value* mycomp = getPntrToComponent(0);
      90        1003 :   outmat.ncols = mycomp->getNumberOfColumns();
      91             :   outmat.resize( mycomp->matrix_bookeeping.size() );
      92    35014890 :   for(unsigned i=0; i<outmat.size(); ++i) {
      93    35013887 :     outmat[i] = mycomp->matrix_bookeeping[i];
      94             :   }
      95        1075 :   for(unsigned i=1; i<getNumberOfComponents(); ++i) {
      96          72 :     getPntrToComponent(i)->copyBookeepingArrayFromArgument( myval );
      97             :   }
      98        1003 : }
      99             : 
     100       12510 : void ActionWithMatrix::transferStashToValues( const std::vector<unsigned>& partialTaskList, const std::vector<double>& stash ) {
     101       12510 :   unsigned ncomp = getNumberOfComponents();
     102       12510 :   unsigned ncols = getPntrToComponent(0)->getNumberOfColumns();
     103       12510 :   unsigned nrows = partialTaskList.size();
     104     1017792 :   for(unsigned i=0; i<nrows; ++i) {
     105     1005282 :     unsigned ncr = getPntrToComponent(0)->getRowLength(partialTaskList[i]);
     106             : #ifndef NDEBUG
     107             :     for(unsigned k=1; k<ncomp; ++k) {
     108             :       plumed_assert( ncr == getPntrToComponent(k)->getRowLength(partialTaskList[i]) );
     109             :     }
     110             : #endif
     111    70602539 :     for(unsigned j=0; j<ncr; ++j) {
     112   202229416 :       for(unsigned k=0; k<ncomp; ++k) {
     113   132632159 :         getPntrToComponent(k)->set( partialTaskList[i]*ncols+j, stash[ncomp*ncols*partialTaskList[i]+j*ncomp+k] );
     114             :       }
     115             :     }
     116             :   }
     117       12510 : }
     118             : 
     119       11626 : void ActionWithMatrix::transferForcesToStash( const std::vector<unsigned>& partialTaskList, std::vector<double>& stash ) const {
     120       11626 :   unsigned ncomp = getNumberOfComponents();
     121       11626 :   unsigned ncols = getConstPntrToComponent(0)->getNumberOfColumns();
     122       11626 :   unsigned nrows = partialTaskList.size();
     123      401973 :   for(unsigned i=0; i<nrows; ++i) {
     124      390347 :     unsigned ncr = getConstPntrToComponent(0)->getRowLength(partialTaskList[i]);
     125             : #ifndef NDEBUG
     126             :     for(unsigned k=1; k<ncomp; ++k) {
     127             :       plumed_assert( ncr == getConstPntrToComponent(k)->getRowLength(partialTaskList[i]) );
     128             :     }
     129             : #endif
     130    15096527 :     for(unsigned j=0; j<ncr; ++j) {
     131    47547231 :       for(unsigned k=0; k<ncomp; ++k) {
     132    32841051 :         stash[ncomp*ncols*partialTaskList[i]+j*ncomp+k] = getConstPntrToComponent(k)->getForce( partialTaskList[i]*ncols+j );
     133             :       }
     134             :     }
     135             :   }
     136       11626 : }
     137             : 
     138             : }