diff --git a/src/core/ActionWithVector.cpp b/src/core/ActionWithVector.cpp index 0e6b9b7212..ceb13ec7ac 100644 --- a/src/core/ActionWithVector.cpp +++ b/src/core/ActionWithVector.cpp @@ -20,6 +20,7 @@ along with plumed. If not, see . +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ #include "ActionWithVector.h" +#include "ActionWithMatrix.h" #include "PlumedMain.h" #include "ActionSet.h" #include "tools/OpenMP.h" @@ -721,11 +722,14 @@ bool ActionWithVector::getNumberOfStoredValues( Value* startat, unsigned& nvals, } void ActionWithVector::runTask( const unsigned& current, MultiValue& myvals ) const { + const ActionWithMatrix* am = dynamic_cast(this); myvals.setTaskIndex(current); myvals.vector_call=true; performTask( current, myvals ); for(unsigned i=0; igetRank()!=1 || myval->hasDerivatives() || !myval->valueIsStored() ) continue; - Value* myv = const_cast( myval ); myv->set( current, myvals.get( myval->getPositionInStream() ) ); + if( am || myval->hasDerivatives() || !myval->valueIsStored() ) continue; + Value* myv = const_cast( myval ); + if( getName()=="RMSD_VECTOR" && myv->getRank()==2 ) continue; + myv->set( current, myvals.get( myval->getPositionInStream() ) ); } if( action_to_do_after ) action_to_do_after->runTask( current, myvals ); } diff --git a/src/core/Value.cpp b/src/core/Value.cpp index 5cd438cf89..d91692c80c 100644 --- a/src/core/Value.cpp +++ b/src/core/Value.cpp @@ -306,6 +306,13 @@ void Value::reshapeMatrixStore( const unsigned& n ) { } } +void Value::copyBookeepingArrayFromArgument( Value* myarg ) { + plumed_dbg_assert( shape.size()==2 && !hasDeriv ); + ncols = myarg->getNumberOfColumns(); matrix_bookeeping.resize( myarg->matrix_bookeeping.size() ); + for(unsigned i=0; imatrix_bookeeping[i]; + data.resize( shape[0]*ncols ); inputForce.resize( shape[0]*ncols ); +} + void Value::setPositionInMatrixStash( const unsigned& p ) { plumed_dbg_assert( shape.size()==2 && !hasDeriv ); matpos=p; diff --git a/src/core/Value.h b/src/core/Value.h index ceb2fabfa7..dfbe885a4a 100644 --- a/src/core/Value.h +++ b/src/core/Value.h @@ -191,6 +191,8 @@ class Value { void buildDataStore( const bool forprint=false ); /// Reshape the storage for sparse matrices void reshapeMatrixStore( const unsigned& n ); +/// Copy the matrix bookeeping stuff + void copyBookeepingArrayFromArgument( Value* myarg ); /// Set the symmetric flag equal true for this matrix void setSymmetric( const bool& sym ); /// Get the total number of scalars that are stored here diff --git a/src/crystdistrib/QuaternionBondProductMatrix.cpp b/src/crystdistrib/QuaternionBondProductMatrix.cpp index dfada2bfe2..897a0ab046 100644 --- a/src/crystdistrib/QuaternionBondProductMatrix.cpp +++ b/src/crystdistrib/QuaternionBondProductMatrix.cpp @@ -143,11 +143,12 @@ void QuaternionBondProductMatrix::performTask( const std::string& controller, co //[dit/dw1 dit/di1 dit/dj1 dit/dk1] etc, and dqt[1] is w.r.t the vector-turned-quaternion called bond // Retrieve the quaternion - for(unsigned i=0; i<4; ++i) quat[i] = getArgumentElement( i, index1, myvals ); + for(unsigned i=0; i<4; ++i) quat[i] = getPntrToArgument(i)->get(index1); // Retrieve the components of the matrix - double weight = getElementOfMatrixArgument( 4, index1, ind2, myvals ); - for(unsigned i=1; i<4; ++i) bond[i] = getElementOfMatrixArgument( 4+i, index1, ind2, myvals ); + unsigned find = getPntrToArgument(4)->getIndexInStore( index1*getPntrToArgument(4)->getShape()[1] + ind2 ); + double weight = getPntrToArgument(4)->get(find, false ); + for(unsigned i=1; i<4; ++i) bond[i] = getPntrToArgument(4+i)->get(find, false ); // calculate normalization factor bond[0]=0.0; @@ -156,16 +157,11 @@ void QuaternionBondProductMatrix::performTask( const std::string& controller, co double normFac3 = normFac*normFac*normFac; //I hold off on normalizing because this can be done at the very end, and it makes the derivatives with respect to 'bond' more simple - - std::vector quat_conj(4); quat_conj[0] = quat[0]; quat_conj[1] = -1*quat[1]; quat_conj[2] = -1*quat[2]; quat_conj[3] = -1*quat[3]; //make a conjugate of q1 my own sanity - - - -//q1_conj * r first, while keep track of derivs + //q1_conj * r first, while keep track of derivs double pref=1; double conj=1; double pref2=1; @@ -236,21 +232,20 @@ void QuaternionBondProductMatrix::performTask( const std::string& controller, co pref2=1; unsigned base=0; for(unsigned i=0; i<4; ++i) { if( i>0 ) {pref=-1; pref2=-1;} - myvals.addValue( getConstPntrToComponent(0)->getPositionInStream(), normFac*pref*quatTemp[i]*quat[i] ); + myvals.addValue( 0, normFac*pref*quatTemp[i]*quat[i] ); wf+=normFac*pref*quatTemp[i]*quat[i]; if( doNotCalculateDerivatives() ) continue ; tempDot=(dotProduct(Vector4d(quat[0],-quat[1],-quat[2],-quat[3]), dqt[0].getCol(i)) + pref2*quatTemp[i])*normFac; - addDerivativeOnVectorArgument( stored[i], 0, i, index1, tempDot, myvals); + myvals.addDerivative( 0, base + index1, tempDot ); myvals.updateIndex( 0, base + index1 ); base += getPntrToArgument(i)->getNumberOfStoredValues(); } //had to split because bond's derivatives depend on the value of the overall quaternion component if( !doNotCalculateDerivatives() ) { - unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream(); for(unsigned i=0; i<4; ++i) { tempDot=dotProduct(Vector4d(quat[0],-quat[1],-quat[2],-quat[3]), dqt[1].getCol(i))*normFac; if( i>0 ) { plumed_assert( !stored[4+i] ); unsigned find = getPntrToArgument(4+i)->getIndexInStore( index1*getPntrToArgument(4+i)->getShape()[1] + ind2 ); - myvals.addDerivative( ostrn, base + find, tempDot ); myvals.updateIndex( ostrn, base + find ); + myvals.addDerivative( 0, base + find, tempDot ); myvals.updateIndex( 0, base + find ); } base += getPntrToArgument(4+i)->getNumberOfStoredValues(); } @@ -262,23 +257,22 @@ void QuaternionBondProductMatrix::performTask( const std::string& controller, co for (unsigned i=0; i<4; i++) { if(i==3) pref=-1; else pref=1; - myvals.addValue( getConstPntrToComponent(1)->getPositionInStream(), normFac*pref*quatTemp[i]*quat[(5-i)%4]); + myvals.addValue( 1, normFac*pref*quatTemp[i]*quat[(5-i)%4]); xf+=normFac*pref*quatTemp[i]*quat[(5-i)%4]; if(i==2) pref2=-1; else pref2=1; if( doNotCalculateDerivatives() ) continue ; tempDot=(dotProduct(Vector4d(quat[1],quat[0],quat[3],-quat[2]), dqt[0].getCol(i)) + pref2*quatTemp[(5-i)%4])*normFac; - addDerivativeOnVectorArgument( stored[i], 1, i, index1, tempDot, myvals); + myvals.addDerivative( 1, base + index1, tempDot ); myvals.updateIndex( 1, base + index1 ); base += getPntrToArgument(i)->getNumberOfStoredValues(); } if( !doNotCalculateDerivatives() ) { - unsigned ostrn = getConstPntrToComponent(1)->getPositionInStream(); for(unsigned i=0; i<4; ++i) { tempDot=dotProduct(Vector4d(quat[1],quat[0],quat[3],-quat[2]), dqt[1].getCol(i))*normFac; if( i>0 ) { plumed_assert( !stored[4+i] ); unsigned find = getPntrToArgument(4+i)->getIndexInStore( index1*getPntrToArgument(4+i)->getShape()[1] + ind2 ); - myvals.addDerivative( ostrn, base + find, tempDot+(-bond[i]*normFac*normFac*xf) ); myvals.updateIndex( ostrn, base + find ); + myvals.addDerivative( 1, base + find, tempDot+(-bond[i]*normFac*normFac*xf) ); myvals.updateIndex( 1, base + find ); } base += getPntrToArgument(4+i)->getNumberOfStoredValues(); } @@ -294,21 +288,20 @@ void QuaternionBondProductMatrix::performTask( const std::string& controller, co if (i==3) pref2=-1; else pref2=1; - myvals.addValue( getConstPntrToComponent(2)->getPositionInStream(), normFac*pref*quatTemp[i]*quat[(i+2)%4]); + myvals.addValue( 2, normFac*pref*quatTemp[i]*quat[(i+2)%4]); yf+=normFac*pref*quatTemp[i]*quat[(i+2)%4]; if( doNotCalculateDerivatives() ) continue ; tempDot=(dotProduct(Vector4d(quat[2],-quat[3],quat[0],quat[1]), dqt[0].getCol(i)) + pref2*quatTemp[(i+2)%4])*normFac; - addDerivativeOnVectorArgument( stored[i], 2, i, index1, tempDot, myvals); + myvals.addDerivative( 2, base + index1, tempDot ); myvals.updateIndex( 2, base + index1 ); base += getPntrToArgument(i)->getNumberOfStoredValues(); } if( !doNotCalculateDerivatives() ) { - unsigned ostrn = getConstPntrToComponent(2)->getPositionInStream(); for(unsigned i=0; i<4; ++i) { tempDot=dotProduct(Vector4d(quat[2],-quat[3],quat[0],quat[1]), dqt[1].getCol(i))*normFac; if( i>0 ) { plumed_assert( !stored[4+i] ); unsigned find = getPntrToArgument(4+i)->getIndexInStore( index1*getPntrToArgument(4+i)->getShape()[1] + ind2 ); - myvals.addDerivative( ostrn, base + find, tempDot+(-bond[i]*normFac*normFac*yf) ); myvals.updateIndex( ostrn, base + find ); + myvals.addDerivative( 2, base + find, tempDot+(-bond[i]*normFac*normFac*yf) ); myvals.updateIndex( 2, base + find ); } base += getPntrToArgument(4+i)->getNumberOfStoredValues(); } @@ -323,22 +316,21 @@ void QuaternionBondProductMatrix::performTask( const std::string& controller, co if(i==1) pref2=-1; else pref2=1; - myvals.addValue( getConstPntrToComponent(3)->getPositionInStream(), normFac*pref*quatTemp[i]*quat[(3-i)]); + myvals.addValue( 3, normFac*pref*quatTemp[i]*quat[(3-i)]); zf+=normFac*pref*quatTemp[i]*quat[(3-i)]; if( doNotCalculateDerivatives() ) continue ; tempDot=(dotProduct(Vector4d(quat[3],quat[2],-quat[1],quat[0]), dqt[0].getCol(i)) + pref2*quatTemp[(3-i)])*normFac; - addDerivativeOnVectorArgument( stored[i], 3, i, index1, tempDot, myvals); + myvals.addDerivative( 3, base + index1, tempDot ); myvals.updateIndex( 3, base + index1 ); base += getPntrToArgument(i)->getNumberOfStoredValues(); } if( doNotCalculateDerivatives() ) return ; - unsigned ostrn = getConstPntrToComponent(3)->getPositionInStream(); for(unsigned i=0; i<4; ++i) { tempDot=dotProduct(Vector4d(quat[3],quat[2],-quat[1],quat[0]), dqt[1].getCol(i))*normFac; if( i>0 ) { plumed_assert( !stored[4+i] ); unsigned find = getPntrToArgument(4+i)->getIndexInStore( index1*getPntrToArgument(4+i)->getShape()[1] + ind2 ); - myvals.addDerivative( ostrn, base + find, tempDot+(-bond[i]*normFac*normFac*zf) ); myvals.updateIndex( ostrn, base + find ); + myvals.addDerivative( 3, base + find, tempDot+(-bond[i]*normFac*normFac*zf) ); myvals.updateIndex( 3, base + find ); } base += getPntrToArgument(4+i)->getNumberOfStoredValues(); } diff --git a/src/matrixtools/MatrixTimesVector.cpp b/src/matrixtools/MatrixTimesVector.cpp index 8af9a84550..a9f860d29d 100644 --- a/src/matrixtools/MatrixTimesVector.cpp +++ b/src/matrixtools/MatrixTimesVector.cpp @@ -19,7 +19,7 @@ You should have received a copy of the GNU Lesser General Public License along with plumed. If not, see . +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#include "core/ActionWithMatrix.h" +#include "core/ActionWithVector.h" #include "core/ActionRegister.h" //+PLUMEDOC MCOLVAR MATRIX_VECTOR_PRODUCT @@ -34,32 +34,25 @@ Calculate the product of the matrix and the vector namespace PLMD { namespace matrixtools { -class MatrixTimesVector : public ActionWithMatrix { +class MatrixTimesVector : public ActionWithVector { private: bool sumrows; - unsigned nderivatives; - std::vector stored_arg; public: static void registerKeywords( Keywords& keys ); explicit MatrixTimesVector(const ActionOptions&); std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ; - unsigned getNumberOfColumns() const override { plumed_error(); } void getNumberOfStreamedDerivatives( unsigned& nderivatives, Value* stopat ) override ; unsigned getNumberOfDerivatives() override ; void prepare() override ; + void calculate() override ; void performTask( const unsigned& task_index, MultiValue& myvals ) const override ; int checkTaskIsActive( const unsigned& itask ) const override ; - bool isInSubChain( unsigned& nder ) override { nder = arg_deriv_starts[0]; return true; } - void setupForTask( const unsigned& task_index, std::vector& indices, MultiValue& myvals ) const ; - void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override; - void runEndOfRowJobs( const unsigned& ind, const std::vector & indices, MultiValue& myvals ) const override ; - void updateAdditionalIndices( const unsigned& ostrn, MultiValue& myvals ) const override ; }; PLUMED_REGISTER_ACTION(MatrixTimesVector,"MATRIX_VECTOR_PRODUCT") void MatrixTimesVector::registerKeywords( Keywords& keys ) { - ActionWithMatrix::registerKeywords(keys); keys.use("ARG"); + ActionWithVector::registerKeywords(keys); keys.use("ARG"); keys.setValueDescription("the vector that is obtained by taking the product between the matrix and the vector that were input"); keys.add("hidden","MASKED_INPUT_ALLOWED","turns on that you are allowed to use masked inputs "); ActionWithValue::useCustomisableComponents(keys); @@ -83,7 +76,7 @@ std::string MatrixTimesVector::getOutputComponentDescription( const std::string& MatrixTimesVector::MatrixTimesVector(const ActionOptions&ao): Action(ao), - ActionWithMatrix(ao), + ActionWithVector(ao), sumrows(false) { if( getNumberOfArguments()<2 ) error("Not enough arguments specified"); @@ -112,9 +105,6 @@ MatrixTimesVector::MatrixTimesVector(const ActionOptions&ao): } getPntrToArgument(n)->buildDataStore(); - ActionWithVector* av=dynamic_cast( getPntrToArgument(0)->getPntrToAction() ); - if( av ) done_in_chain=canBeAfterInChain( av ); - if( getNumberOfArguments()==2 ) { addValue( shape ); setNotPeriodic(); } else { @@ -144,9 +134,6 @@ MatrixTimesVector::MatrixTimesVector(const ActionOptions&ao): getPntrToArgument(i)->buildDataStore(); } - ActionWithVector* av=dynamic_cast( getPntrToArgument(0)->getPntrToAction() ); - if( av ) done_in_chain=canBeAfterInChain( av ); - for(unsigned i=1; igetName(); if( name.find_first_of(".")!=std::string::npos ) { std::size_t dot=name.find_first_of("."); name = name.substr(dot+1); } @@ -154,12 +141,12 @@ MatrixTimesVector::MatrixTimesVector(const ActionOptions&ao): } } else error("You should either have one vector or one matrix in input"); - nderivatives = buildArgumentStore(0); - std::string headstr=getFirstActionInChain()->getLabel(); stored_arg.resize( getNumberOfArguments() ); - for(unsigned i=0; iignoreStoredValue( headstr ); + unsigned nderivatives = buildArgumentStore(0); } unsigned MatrixTimesVector::getNumberOfDerivatives() { + unsigned nderivatives=0; + for(unsigned i=0; igetNumberOfStoredValues(); return nderivatives; } @@ -170,8 +157,6 @@ void MatrixTimesVector::prepare() { } void MatrixTimesVector::getNumberOfStreamedDerivatives( unsigned& nderivatives, Value* stopat ) { - if( actionInChain() ) { ActionWithVector::getNumberOfStreamedDerivatives( nderivatives, stopat ); return; } - nderivatives = 0; for(unsigned i=0; igetNumberOfColumns(); unsigned nmat = mymat->getRowLength(task_index); double val=0; for(unsigned j=0; jget( task_index*ncol + j, false ); - unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream(); - myvals.setValue( ostrn, val ); + myvals.setValue( i, val ); // And the derivatives if( doNotCalculateDerivatives() ) continue; unsigned dloc = arg_deriv_starts[i] + task_index*ncol; for(unsigned j=0; jgetRank()==1 ) { @@ -224,8 +210,7 @@ void MatrixTimesVector::performTask( const unsigned& task_index, MultiValue& myv for(unsigned i=0; iget( task_index*ncol + j, false )*myvec->get( mymat->getRowIndex( task_index, j ) ); - unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream(); - myvals.setValue( ostrn, val ); + myvals.setValue( i, val ); // And the derivatives if( doNotCalculateDerivatives() ) continue; @@ -234,8 +219,8 @@ void MatrixTimesVector::performTask( const unsigned& task_index, MultiValue& myv unsigned kind = mymat->getRowIndex( task_index, j ); double vecval = myvec->get( kind ); double matval = mymat->get( task_index*ncol + j, false ); - myvals.addDerivative( ostrn, dloc + j, vecval ); myvals.updateIndex( ostrn, dloc + j ); - myvals.addDerivative( ostrn, arg_deriv_starts[i+1] + kind, matval ); myvals.updateIndex( ostrn, arg_deriv_starts[i+1] + kind ); + myvals.addDerivative( i, dloc + j, vecval ); myvals.updateIndex( i, dloc + j ); + myvals.addDerivative( i, arg_deriv_starts[i+1] + kind, matval ); myvals.updateIndex( i, arg_deriv_starts[i+1] + kind ); } } } else { @@ -245,8 +230,7 @@ void MatrixTimesVector::performTask( const unsigned& task_index, MultiValue& myv unsigned ncol = mymat->getNumberOfColumns(); unsigned nmat = mymat->getRowLength(task_index); double val=0; for(unsigned j=0; jget( task_index*ncol + j, false )*myvec->get( mymat->getRowIndex( task_index, j ) ); - unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream(); - myvals.setValue( ostrn, val ); + myvals.setValue( i, val ); // And the derivatives if( doNotCalculateDerivatives() ) continue; @@ -256,95 +240,12 @@ void MatrixTimesVector::performTask( const unsigned& task_index, MultiValue& myv unsigned kind = mymat->getRowIndex( task_index, j ); double vecval = myvec->get( kind ); double matval = mymat->get( task_index*ncol + j, false ); - myvals.addDerivative( ostrn, dloc + j, vecval ); myvals.updateIndex( ostrn, dloc + j ); - myvals.addDerivative( ostrn, arg_deriv_starts[n] + kind, matval ); myvals.updateIndex( ostrn, arg_deriv_starts[n] + kind ); + myvals.addDerivative( i, dloc + j, vecval ); myvals.updateIndex( i, dloc + j ); + myvals.addDerivative( i, arg_deriv_starts[n] + kind, matval ); myvals.updateIndex( i, arg_deriv_starts[n] + kind ); } } } } -void MatrixTimesVector::setupForTask( const unsigned& task_index, std::vector& indices, MultiValue& myvals ) const { - unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(0)->getRowLength(task_index); - if( indices.size()!=size_v+1 ) indices.resize( size_v + 1 ); - for(unsigned i=0; igetLabel(); - unsigned ind2 = index2; if( index2>=getPntrToArgument(0)->getShape()[0] ) ind2 = index2 - getPntrToArgument(0)->getShape()[0]; - if( sumrows ) { - unsigned n=getNumberOfArguments()-1; double matval = 0; - for(unsigned i=0; igetPositionInStream(); - Value* myarg = getPntrToArgument(i); - if( !myarg->valueHasBeenSet() || myarg->ignoreStoredValue(headstr) ) myvals.addValue( ostrn, myvals.get( myarg->getPositionInStream() ) ); - else myvals.addValue( ostrn, myarg->get( index1*myarg->getNumberOfColumns() + ind2, true ) ); - // Now lets work out the derivatives - if( doNotCalculateDerivatives() ) continue; - addDerivativeOnMatrixArgument( stored_arg[i], i, i, index1, ind2, 1.0, myvals ); - } - } else if( getPntrToArgument(1)->getRank()==1 ) { - double matval = 0; Value* myarg = getPntrToArgument(0); unsigned vcol = ind2; - if( !myarg->valueHasBeenSet() || myarg->ignoreStoredValue(headstr) ) matval = myvals.get( myarg->getPositionInStream() ); - else { - matval = myarg->get( index1*myarg->getNumberOfColumns() + ind2, true ); - vcol = getPntrToArgument(0)->getRowIndex( index1, ind2 ); - } - for(unsigned i=0; igetPositionInStream(); - double vecval=getArgumentElement( i+1, vcol, myvals ); - // And add this part of the product - myvals.addValue( ostrn, matval*vecval ); - // Now lets work out the derivatives - if( doNotCalculateDerivatives() ) continue; - addDerivativeOnMatrixArgument( stored_arg[0], i, 0, index1, ind2, vecval, myvals ); addDerivativeOnVectorArgument( stored_arg[i+1], i, i+1, vcol, matval, myvals ); - } - } else { - unsigned n=getNumberOfArguments()-1; double matval = 0; unsigned vcol = ind2; - for(unsigned i=0; igetPositionInStream(); - Value* myarg = getPntrToArgument(i); - if( !myarg->valueHasBeenSet() || myarg->ignoreStoredValue(headstr) ) matval = myvals.get( myarg->getPositionInStream() ); - else { - matval = myarg->get( index1*myarg->getNumberOfColumns() + ind2, true ); - vcol = getPntrToArgument(i)->getRowIndex( index1, ind2 ); - } - double vecval=getArgumentElement( n, vcol, myvals ); - // And add this part of the product - myvals.addValue( ostrn, matval*vecval ); - // Now lets work out the derivatives - if( doNotCalculateDerivatives() ) continue; - addDerivativeOnMatrixArgument( stored_arg[i], i, i, index1, ind2, vecval, myvals ); addDerivativeOnVectorArgument( stored_arg[n], i, n, vcol, matval, myvals ); - } - } -} - -void MatrixTimesVector::runEndOfRowJobs( const unsigned& ind, const std::vector & indices, MultiValue& myvals ) const { - if( doNotCalculateDerivatives() ) return ; - - if( getPntrToArgument(1)->getRank()==1 ) { - unsigned istrn = getPntrToArgument(0)->getPositionInMatrixStash(); - std::vector& mat_indices( myvals.getMatrixRowDerivativeIndices( istrn ) ); - for(unsigned j=0; jgetPositionInStream(); - for(unsigned i=0; igetPositionInMatrixStash(); - unsigned ostrn = getConstPntrToComponent(j)->getPositionInStream(); - std::vector& mat_indices( myvals.getMatrixRowDerivativeIndices( istrn ) ); - for(unsigned i=0; igetRank()==1 ) n = 1; - unsigned nvals = getPntrToArgument(n)->getNumberOfValues(); - for(unsigned i=0; i