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 "core/ActionWithVector.h"
      23             : #include "core/ParallelTaskManager.h"
      24             : #include "core/ActionRegister.h"
      25             : #include <limits>
      26             : 
      27             : //+PLUMEDOC FUNCTION MATRIX_PRODUCT_DIAGONAL
      28             : /*
      29             : Calculate the product of two matrices and return a vector that contains the diagonal elements of the ouptut vector
      30             : 
      31             : This action allows you to [multiply](https://en.wikipedia.org/wiki/Matrix_multiplication) two matrices and recover the diagonal of the
      32             : product of the matrices. The following input shows an example where two contact matrices are multiplied together.
      33             : 
      34             : ```plumed
      35             : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.1}
      36             : c2: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.2}
      37             : m: MATRIX_PRODUCT_DIAGONAL ARG=c1,c2
      38             : PRINT ARG=m FILE=colvar
      39             : ```
      40             : 
      41             : As you can see from the documentation for the shortcut [EUCLIDEAN_DISTANCE](EUCLIDEAN_DISTANCE.md) this action is useful for computing distances between
      42             : vectors.
      43             : 
      44             : */
      45             : //+ENDPLUMEDOC
      46             : 
      47             : namespace PLMD {
      48             : namespace refdist {
      49             : 
      50             : class MatrixProductDiagonalInput {
      51             : public:
      52             :   // the assignment is customized (see below): copy construction is the default one
      53             :   MatrixProductDiagonalInput() = default;
      54             :   MatrixProductDiagonalInput(const MatrixProductDiagonalInput&) = default;
      55             :   MatrixProductDiagonalInput& operator=( const MatrixProductDiagonalInput& m ) {
      56             :     // Don't need to copy A and B here as they will be set in calculate
      57             :     return *this;
      58             :   }
      59             : };
      60             : 
      61             : class MatrixProductDiagonal : public ActionWithVector {
      62             : public:
      63             :   using input_type = MatrixProductDiagonalInput;
      64             :   using PTM = ParallelTaskManager<MatrixProductDiagonal>;
      65             :   typedef typename PTM::ParallelActionsInput ParallelActionsInput;
      66             :   typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
      67             : private:
      68             :   PTM taskmanager;
      69             : public:
      70             :   static void registerKeywords( Keywords& keys );
      71             :   explicit MatrixProductDiagonal(const ActionOptions&);
      72             :   unsigned getNumberOfDerivatives() override ;
      73             :   void calculate() override ;
      74             :   void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
      75             :   static void performTask( std::size_t task_index, const MatrixProductDiagonalInput& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output );
      76             :   static int getNumberOfValuesPerTask( std::size_t task_index, const MatrixProductDiagonalInput& actiondata );
      77             :   static void getForceIndices( std::size_t task_index, std::size_t colno, std::size_t ntotal_force, const MatrixProductDiagonalInput& actiondata, const ParallelActionsInput& input, ForceIndexHolder force_indices );
      78             : };
      79             : 
      80             : PLUMED_REGISTER_ACTION(MatrixProductDiagonal,"MATRIX_PRODUCT_DIAGONAL")
      81             : 
      82         114 : void MatrixProductDiagonal::registerKeywords( Keywords& keys ) {
      83         114 :   ActionWithVector::registerKeywords(keys);
      84         228 :   keys.addInputKeyword("compulsory","ARG","vector/matrix","the two vectors/matrices whose product are to be taken");
      85         228 :   keys.setValueDescription("scalar/vector","a vector containing the diagonal elements of the matrix that obtaned by multiplying the two input matrices together");
      86         114 :   PTM::registerKeywords( keys );
      87         114 : }
      88             : 
      89          55 : MatrixProductDiagonal::MatrixProductDiagonal(const ActionOptions&ao):
      90             :   Action(ao),
      91             :   ActionWithVector(ao),
      92          55 :   taskmanager(this) {
      93          55 :   if( getNumberOfArguments()!=2 ) {
      94           0 :     error("should be two vectors or matrices in argument to this action");
      95             :   }
      96             : 
      97         110 :   const unsigned ncols = [&]()->unsigned{
      98          55 :     if( getPntrToArgument(0)->getRank()==1 ) {
      99           2 :       if( getPntrToArgument(0)->hasDerivatives() ) {
     100           0 :         error("first argument to this action should be a vector or matrix");
     101             :       }
     102             :       return 1;
     103          53 :     } else if( getPntrToArgument(0)->getRank()==2 ) {
     104          53 :       if( getPntrToArgument(0)->hasDerivatives() ) {
     105           0 :         error("first argument to this action should be a matrix");
     106             :       }
     107          53 :       return getPntrToArgument(0)->getShape()[1];
     108             :     }
     109             :     //you should not be here!
     110           0 :     error("first argument to this action should be a vector or matrix");
     111             :     //this is a dummy returns(in fact the code won't get here.
     112             :     //this lambda is there due to a -Wmaybe-uninitialized
     113             :     return std::numeric_limits<unsigned>::max();
     114          55 :   }();
     115             : 
     116          55 :   if( getPntrToArgument(1)->getRank()==1 ) {
     117          32 :     if( getPntrToArgument(1)->hasDerivatives() ) {
     118           0 :       error("second argument to this action should be a vector or matrix");
     119             :     }
     120          32 :     if( ncols!=getPntrToArgument(1)->getShape()[0] ) {
     121           0 :       error("number of columns in first matrix does not equal number of elements in vector");
     122             :     }
     123          32 :     if( getPntrToArgument(0)->getShape()[0]!=1 ) {
     124           0 :       error("matrix output by this action must be square");
     125             :     }
     126          32 :     addValueWithDerivatives();
     127          32 :     setNotPeriodic();
     128             :   } else {
     129          23 :     if( getPntrToArgument(1)->getRank()!=2 || getPntrToArgument(1)->hasDerivatives() ) {
     130           0 :       error("second argument to this action should be a vector or a matrix");
     131             :     }
     132          23 :     if( ncols!=getPntrToArgument(1)->getShape()[0] ) {
     133           0 :       error("number of columns in first matrix does not equal number of rows in second matrix");
     134             :     }
     135          23 :     if( getPntrToArgument(0)->getShape()[0]!=getPntrToArgument(1)->getShape()[1] ) {
     136           0 :       error("matrix output by this action must be square");
     137             :     }
     138          23 :     std::vector<std::size_t> shape(1);
     139          23 :     shape[0]=getPntrToArgument(0)->getShape()[0];
     140          23 :     addValue( shape );
     141          23 :     setNotPeriodic();
     142          23 :     taskmanager.setupParallelTaskManager( ncols + getPntrToArgument(1)->getShape()[0], 0 );
     143             :     taskmanager.setActionInput( MatrixProductDiagonalInput() );
     144          23 :   }
     145          55 : }
     146             : 
     147        1239 : unsigned MatrixProductDiagonal::getNumberOfDerivatives() {
     148        1239 :   if( doNotCalculateDerivatives() ) {
     149             :     return 0;
     150             :   }
     151          93 :   return getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();;
     152             : }
     153             : 
     154        2792 : void MatrixProductDiagonal::calculate() {
     155        2792 :   if( getPntrToArgument(0)->getRank()==1 && getPntrToArgument(1)->getRank()==1 ) {
     156           0 :     double val1 = getPntrToArgument(0)->get(0);
     157           0 :     double val2 = getPntrToArgument(1)->get(0);
     158           0 :     Value* myval = getPntrToComponent(0);
     159           0 :     myval->set( val1*val2 );
     160             : 
     161           0 :     if( doNotCalculateDerivatives() ) {
     162             :       return;
     163             :     }
     164             : 
     165             :     myval->setDerivative( 0, val2 );
     166             :     myval->setDerivative( 1, val1 );
     167        2792 :   } else if( getPntrToArgument(1)->getRank()==1 ) {
     168             :     Value* arg1 = getPntrToArgument(0);
     169             :     Value* arg2 = getPntrToArgument(1);
     170             :     unsigned nmult = arg1->getRowLength(0);
     171        1179 :     unsigned nvals1 = arg1->getNumberOfValues();
     172             : 
     173             :     double matval=0;
     174        1179 :     Value* myval = getPntrToComponent(0);
     175        4674 :     for(unsigned i=0; i<nmult; ++i) {
     176        3495 :       unsigned kind = arg1->getRowIndex( 0, i );
     177        3495 :       double val1 = arg1->get( i, false );
     178        3495 :       double val2 = arg2->get( kind );
     179        3495 :       matval += val1*val2;
     180             : 
     181        3495 :       if( doNotCalculateDerivatives() ) {
     182        3420 :         continue;
     183             :       }
     184             : 
     185             :       myval->addDerivative( kind, val2 );
     186          75 :       myval->addDerivative( nvals1 + kind, val1 );
     187             :     }
     188             :     myval->set( matval );
     189             :   } else {
     190        1613 :     taskmanager.runAllTasks();
     191             :   }
     192             : }
     193             : 
     194        1204 : void MatrixProductDiagonal::applyNonZeroRankForces( std::vector<double>& outforces ) {
     195        1204 :   taskmanager.applyForces( outforces );
     196        1204 : }
     197             : 
     198      130004 : void MatrixProductDiagonal::performTask( std::size_t task_index,
     199             :     const MatrixProductDiagonalInput& actiondata,
     200             :     ParallelActionsInput& input,
     201             :     ParallelActionsOutput& output ) {
     202      130004 :   auto arg0=ArgumentBookeepingHolder::create( 0, input );
     203      130004 :   auto arg1=ArgumentBookeepingHolder::create( 1, input );
     204      130004 :   std::size_t fpos = task_index*(1+arg0.ncols);
     205      130004 :   std::size_t nmult = arg0.bookeeping[fpos];
     206      130004 :   std::size_t vstart = task_index*arg0.ncols;
     207      130004 :   output.values[0] = 0;
     208      130004 :   if( arg1.ncols<arg1.shape[1] ) {
     209           0 :     plumed_merror("multiplying by a sparse matrix is not implemented as I don't think it is needed");
     210             :   } else {
     211      130004 :     std::size_t base = arg1.start + task_index;
     212     3905172 :     for(unsigned i=0; i<nmult; ++i) {
     213     3775168 :       double val1 = input.inputdata[vstart+i];
     214     3775168 :       double val2 = input.inputdata[ base + arg1.ncols*arg0.bookeeping[fpos+1+i] ];
     215     3775168 :       output.values[0] += val1*val2;
     216     3775168 :       output.derivatives[i] = val2;
     217     3775168 :       output.derivatives[nmult + i] = val1;
     218             :     }
     219             :   }
     220      130004 : }
     221             : 
     222       56018 : int MatrixProductDiagonal::getNumberOfValuesPerTask( std::size_t task_index, const MatrixProductDiagonalInput& actiondata ) {
     223       56018 :   return 1;
     224             : }
     225             : 
     226       56018 : void MatrixProductDiagonal::getForceIndices( std::size_t task_index,
     227             :     std::size_t colno,
     228             :     std::size_t ntotal_force,
     229             :     const MatrixProductDiagonalInput& actiondata,
     230             :     const ParallelActionsInput& input,
     231             :     ForceIndexHolder force_indices ) {
     232       56018 :   auto arg0=ArgumentBookeepingHolder::create( 0, input );
     233       56018 :   auto arg1=ArgumentBookeepingHolder::create( 1, input );
     234       56018 :   std::size_t fpos = task_index*(1+arg0.ncols);
     235       56018 :   std::size_t nmult = arg0.bookeeping[fpos];
     236       56018 :   std::size_t vstart = task_index*arg0.ncols;
     237       56018 :   if( arg1.ncols<arg1.shape[1] ) {
     238           0 :     plumed_merror("multiplying by a sparse matrix is not implemented as I don't think it is needed");
     239             :   } else {
     240       56018 :     std::size_t base = arg1.start + task_index;
     241     1910866 :     for(unsigned i=0; i<nmult; ++i) {
     242     1854848 :       force_indices.indices[0][i] = vstart + i;
     243     1854848 :       force_indices.indices[0][nmult+i] = base + arg1.ncols*arg0.bookeeping[fpos+1+i];
     244             :     }
     245       56018 :     force_indices.threadsafe_derivatives_end[0] = 2*nmult;
     246       56018 :     force_indices.tot_indices[0] = 2*nmult;
     247             :   }
     248       56018 : }
     249             : 
     250             : }
     251             : }