Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 : Copyright (c) 2014-2017 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 "MatrixOperationBase.h" 23 : #include "core/ActionSetup.h" 24 : #include "core/ActionRegister.h" 25 : 26 : //+PLUMEDOC MCOLVAR INVERT_MATRIX 27 : /* 28 : Calculate the inverse of the input matrix 29 : 30 : \par Examples 31 : 32 : */ 33 : //+ENDPLUMEDOC 34 : 35 : namespace PLMD { 36 : namespace matrixtools { 37 : 38 : class InvertMatrix : public MatrixOperationBase { 39 : private: 40 : bool input_is_constant; 41 : Matrix<double> mymatrix; 42 : Matrix<double> inverse; 43 : public: 44 : static void registerKeywords( Keywords& keys ); 45 : /// Constructor 46 : explicit InvertMatrix(const ActionOptions&); 47 : /// 48 6 : unsigned getNumberOfDerivatives() override { return 0; } 49 : /// Do the calculation 50 : void calculate() override; 51 : /// 52 : void apply() override; 53 0 : double getForceOnMatrixElement( const unsigned& jrow, const unsigned& krow ) const override { plumed_error(); } 54 : }; 55 : 56 : PLUMED_REGISTER_ACTION(InvertMatrix,"INVERT_MATRIX") 57 : 58 18 : void InvertMatrix::registerKeywords( Keywords& keys ) { 59 18 : MatrixOperationBase::registerKeywords( keys ); 60 18 : keys.setValueDescription("the inverse of the input matrix"); 61 18 : } 62 : 63 9 : InvertMatrix::InvertMatrix(const ActionOptions& ao): 64 : Action(ao), 65 : MatrixOperationBase(ao), 66 9 : input_is_constant(false) 67 : { 68 9 : if( getPntrToArgument(0)->getShape()[0]!=getPntrToArgument(0)->getShape()[1] ) error("input matrix should be square"); 69 : 70 9 : ActionSetup* as = dynamic_cast<ActionSetup*>( getPntrToArgument(0)->getPntrToAction() ); 71 9 : if(as) input_is_constant=true; 72 : 73 9 : std::vector<unsigned> shape(2); shape[0]=shape[1]=getPntrToArgument(0)->getShape()[0]; addValue( shape ); 74 9 : setNotPeriodic(); getPntrToComponent(0)->buildDataStore(); getPntrToComponent(0)->reshapeMatrixStore( shape[1] ); 75 9 : mymatrix.resize( shape[0], shape[1] ); inverse.resize( shape[0], shape[1] ); 76 9 : } 77 : 78 9 : void InvertMatrix::calculate() { 79 : // Retrieve the matrix from input 80 9 : retrieveFullMatrix( mymatrix ); 81 : // Now invert the matrix 82 9 : Invert( mymatrix, inverse ); 83 : // And set the inverse 84 9 : unsigned k = 0; Value* myval=getPntrToComponent(0); 85 30 : for(unsigned i=0; i<mymatrix.nrows(); ++i) { 86 78 : for(unsigned j=0; j<mymatrix.ncols(); ++j) { 87 57 : myval->set( k, inverse(i,j) ); k++; 88 : } 89 : } 90 : 91 9 : if( !doNotCalculateDerivatives() && !input_is_constant ) error("derivatives of inverse matrix have not been implemented"); 92 9 : } 93 : 94 0 : void InvertMatrix::apply() { 95 0 : if( doNotCalculateDerivatives() || input_is_constant ) return; 96 0 : error("derivatives of inverse matrix have not been implemented"); 97 : } 98 : 99 : } 100 : }