1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2    Copyright (c) 2014-2021 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 "MultiValue.h"
23 
24 namespace PLMD {
25 
MultiValue(const unsigned & nvals,const unsigned & nder)26 MultiValue::MultiValue( const unsigned& nvals, const unsigned& nder ):
27   values(nvals),
28   nderivatives(nder),
29   derivatives(nvals*nder),
30   tmpval(0),
31   tmpder(nder),
32   atLeastOneSet(false)
33 {
34   std::vector<unsigned> myind( nder );
35   for(unsigned i=0; i<nder; ++i) myind[i]=i;
36   hasDerivatives.createIndexListFromVector( myind );
37 }
38 
resize(const unsigned & nvals,const unsigned & nder)39 void MultiValue::resize( const unsigned& nvals, const unsigned& nder ) {
40   values.resize(nvals); nderivatives=nder; derivatives.resize( nvals*nder );
41   tmpder.resize( nder ); hasDerivatives.clear(); std::vector<unsigned> myind( nder );
42   for(unsigned i=0; i<nder; ++i) myind[i]=i;
43   hasDerivatives.createIndexListFromVector( myind );
44   atLeastOneSet=false;
45 }
46 
clearAll()47 void MultiValue::clearAll() {
48   if( atLeastOneSet && !hasDerivatives.updateComplete() ) hasDerivatives.updateActiveMembers();
49   for(unsigned i=0; i<values.size(); ++i) clear(i);
50   clearTemporyDerivatives(); hasDerivatives.deactivateAll(); atLeastOneSet=false;
51 }
52 
clear(const unsigned & ival)53 void MultiValue::clear( const unsigned& ival ) {
54   values[ival]=0;
55   unsigned base=ival*nderivatives, ndert=hasDerivatives.getNumberActive();
56   for(unsigned i=0; i<ndert; ++i) derivatives[ base+hasDerivatives[i] ]=0.;
57 }
58 
clearTemporyDerivatives()59 void MultiValue::clearTemporyDerivatives() {
60   unsigned ndert=hasDerivatives.getNumberActive(); tmpval=0.;
61   for(unsigned i=0; i<ndert; ++i) tmpder[ hasDerivatives[i] ]=0.;
62 }
63 
chainRule(const unsigned & ival,const unsigned & iout,const unsigned & stride,const unsigned & off,const double & df,const unsigned & bufstart,std::vector<double> & buffer)64 void MultiValue::chainRule( const unsigned& ival, const unsigned& iout, const unsigned& stride, const unsigned& off,
65                             const double& df, const unsigned& bufstart, std::vector<double>& buffer ) {
66   if( !hasDerivatives.updateComplete() ) hasDerivatives.updateActiveMembers();
67 
68   plumed_dbg_assert( off<stride );
69   unsigned base=nderivatives*ival, ndert=hasDerivatives.getNumberActive();
70   unsigned start=bufstart+stride*(nderivatives+1)*iout + stride;
71   for(unsigned i=0; i<ndert; ++i) {
72     unsigned jder=hasDerivatives[i];
73     buffer[start+jder*stride] += df*derivatives[base+jder];
74   }
75 }
76 
copyValues(MultiValue & outvals) const77 void MultiValue::copyValues( MultiValue& outvals ) const {
78   plumed_dbg_assert( values.size()<=outvals.getNumberOfValues() );
79   for(unsigned i=0; i<values.size(); ++i) outvals.setValue( i, values[i] );
80 
81 }
82 
copyDerivatives(MultiValue & outvals)83 void MultiValue::copyDerivatives( MultiValue& outvals ) {
84   plumed_dbg_assert( values.size()<=outvals.getNumberOfValues() && nderivatives<=outvals.getNumberOfDerivatives() );
85   if( !hasDerivatives.updateComplete() ) hasDerivatives.updateActiveMembers();
86 
87   outvals.atLeastOneSet=true; unsigned ndert=hasDerivatives.getNumberActive();
88   for(unsigned j=0; j<ndert; ++j) {
89     unsigned jder=hasDerivatives[j]; outvals.hasDerivatives.activate(jder);
90   }
91 
92   unsigned base=0, obase=0;
93   for(unsigned i=0; i<values.size(); ++i) {
94     for(unsigned j=0; j<ndert; ++j) {
95       unsigned jder=hasDerivatives[j];
96       outvals.derivatives[obase+jder] += derivatives[base+jder];
97     }
98     obase+=outvals.nderivatives; base+=nderivatives;
99   }
100 }
101 
quotientRule(const unsigned & nder,const unsigned & oder)102 void MultiValue::quotientRule( const unsigned& nder, const unsigned& oder ) {
103   plumed_dbg_assert( nder<values.size() && oder<values.size() );
104   if( !hasDerivatives.updateComplete() ) hasDerivatives.updateActiveMembers();
105 
106   unsigned ndert=hasDerivatives.getNumberActive(); double wpref;
107   unsigned obase=oder*nderivatives, nbase=nder*nderivatives;
108 
109   if( fabs(tmpval)>epsilon ) { wpref=1.0/tmpval; }
110   else { wpref=1.0; }
111 
112   double pref = values[nder]*wpref*wpref;
113   for(unsigned j=0; j<ndert; ++j) {
114     unsigned jder=hasDerivatives[j];
115     derivatives[obase+jder] = wpref*derivatives[nbase+jder]  - pref*tmpder[jder];
116   }
117   values[oder] = wpref*values[nder];
118 }
119 
120 }
121