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