work on dg functions

Numeric/fullMatrix.cpp

@@ -37,6 +37,7 @@ extern "C" {
                       std::complex<double> *alpha, std::complex<double> *a, int *lda, 
                       std::complex<double> *x, int *incx, std::complex<double> *beta, 
                       std::complex<double> *y, int *incy);
+  void F77NAME(dscal)(int *n, double *alpha,double *x,  int *incx);
 }
 
 template<> 
@@ -46,6 +47,18 @@ void fullVector<double>::axpy(fullVector<double> &x,double alpha)
   F77NAME(daxpy)(&M, &alpha, x._data,&INCX, _data, &INCY);
 }
 
+template<> 
+void fullMatrix<double>::scale(const double s) {
+  int N = _r*_c;
+  int stride = 1;
+  double ss = s;
+  F77NAME(dscal)(&N, &ss,_data, &stride);
+    /*if(s == 0.)
+      for(int i = 0; i < _r * _c; ++i) _data[i] = 0.;
+    else
+      for(int i = 0; i < _r * _c; ++i) _data[i] *= s;*/
+}
+
 template<> 
 void fullMatrix<double>::mult(const fullMatrix<double> &b, fullMatrix<double> &c) const
 {

Numeric/fullMatrix.h

@@ -144,7 +144,7 @@ class fullMatrix
   }
   fullMatrix() : _own_data(false),_r(0), _c(0), _data(0) {}
   ~fullMatrix() { if(_data && _own_data) delete [] _data; }
-  bool resize(int r, int c) // data will be owned (same as constructor)
+  bool resize(int r, int c, bool resetValue = true) // data will be owned (same as constructor)
   {
     if ((r * c > _r * _c) || !_own_data){
       _r = r;
@@ -152,6 +152,7 @@ class fullMatrix
       if (_own_data && _data) delete[] _data;
       _data = new scalar[_r * _c];
       _own_data = true;
+      if(resetValue)
         scale(0.);
       return true;
     }
@@ -159,6 +160,7 @@ class fullMatrix
       _r = r;
       _c = c;
     }
+    if(resetValue)
       scale(0.);
     return false; // no reallocation
   }
@@ -272,13 +274,17 @@ class fullMatrix
   {
     for(int i = 0; i < _r * _c; i++) _data[i] = m._data[i];
   }
-  inline void scale(const double s)
+  void scale(const double s)
+#if !defined(HAVE_BLAS)
   {
+    
     if(s == 0.)
       for(int i = 0; i < _r * _c; ++i) _data[i] = 0.;
     else
       for(int i = 0; i < _r * _c; ++i) _data[i] *= s;
   }
+#endif
+  ;
   inline void add(const double &a) 
   {
     for(int i = 0; i < _r * _c; ++i) _data[i] += a;

Solver/function.cpp

@@ -9,112 +9,128 @@
   #include "dlfcn.h"
 #endif
 #include "Bindings.h"
+
+
+struct functionReplaceCache {
+  dataCacheMap *map;
+  std::vector <dataCacheDouble*> toReplace;
+  std::vector <dataCacheDouble*> toCompute;
+};
+
+
 function::~function() {
 }
 
 function::function(int nbCol, bool invalidatedOnElement):_nbCol(nbCol), _invalidatedOnElement(invalidatedOnElement){}
 
 void function::addFunctionReplace(functionReplace &fr) {
+  fr._master = this;
   _functionReplaces.push_back(&fr);
 }
-dataCacheDouble::~dataCacheDouble()
-{
-  for(unsigned i = 0; i< functionReplaceCaches.size(); i++) {
-    delete functionReplaceCaches[i];
+void function::setArgument(fullMatrix<double> &v, const function *f, int iMap) {
+  if(f==NULL)
+    throw;
+  arguments.push_back(argument(v, iMap, f));
+  dependencies.insert(dependency(iMap, f));
+  for(std::set<dependency>::iterator it = f->dependencies.begin(); it != f->dependencies.end(); it++) {
+    if (it->iMap> 0 && iMap >0)
+      Msg::Error("consecutive secondary caches");
+    dependencies.insert(dependency(iMap + it->iMap, it->f));
   }
-};
-
-void dataCacheDouble::addMeAsDependencyOf (dataCacheDouble *newDep)
-{
-  _dependOnMe.insert(newDep);
-  newDep->_iDependOn.insert(this);
-  for(std::set<dataCacheDouble*>::iterator it = _iDependOn.begin();
-      it != _iDependOn.end(); it++) {
-    (*it)->_dependOnMe.insert(newDep);
-    newDep->_iDependOn.insert(*it);
+  for (int i = 0; i < _functionReplaces.size(); i++) {
+    functionReplace &replace = *_functionReplaces[i];
+    for (std::set<dependency>::iterator it = replace._fromParent.begin(); it != replace._fromParent.end(); it++) {
+      if (it->iMap> 0 && iMap >0)
+        Msg::Error("consecutive secondary caches");
+      dependencies.insert(dependency(iMap + it->iMap, it->f));
     }
   }
-
-dataCacheDouble::dataCacheDouble(dataCacheMap &m, int nRowByPoint, int nbCol):
-  _cacheMap(m),_value(m.getNbEvaluationPoints()*nRowByPoint,nbCol)
-{
-  _nRowByPoint=nRowByPoint;
-  m.addDataCacheDouble(this,true);
-  _function = NULL;
 }
 
 dataCacheDouble::dataCacheDouble(dataCacheMap *m, function *f):
-  _cacheMap(*m),_value(m->getNbEvaluationPoints(),f->getNbCol())
+  _cacheMap(*m),_value(m->getNbEvaluationPoints(),f->getNbCol()), _function(f)
 {
-  _nRowByPoint=1;
   m->addDataCacheDouble(this, f->isInvalitedOnElement());
-  _function = f;
-  for(unsigned i=0; i<f->_childrenCache.size(); i++) {
-    m->addSecondaryCache(m->newChild()); }
-  _substitutions.resize(f->_substitutedFunctions.size());
-  for(unsigned i=0; i<f->_substitutedFunctions.size(); i++) {
-    function::substitutedFunction s = f->_substitutedFunctions[i];
-    _substitutions[i].first = &m->getSecondaryCache(s.iMap)->substitute(s.f0);
-    _substitutions[i].second = &m->get(s.f1,this);
-  }
-  _dependencies.resize ( _function->arguments.size());
-  for (unsigned int i=0;i<_function->arguments.size();i++) {
-    int iCache = _function->arguments[i].iMap;
-    const function *f = _function->arguments[i].f;
-    _dependencies[i] = &m->getSecondaryCache(iCache)->get(f,this);
-  }
-  for (unsigned i = 0; i < f->_functionReplaces.size(); i++) {
-    functionReplaceCaches.push_back (new functionReplaceCache(m, f->_functionReplaces[i], this)); 
-  }
-  f->registerInDataCacheMap(m, this);
 }
 
-void dataCacheDouble::resize() {
-  _value = fullMatrix<double>(_nRowByPoint==0?1:_nRowByPoint*_cacheMap.getNbEvaluationPoints(),_value.size2());
+void dataCacheDouble::resize(int nrow) {
+  _value.resize(nrow, _value.size2());
+  _valid = false;
 }
+
 void dataCacheDouble::_eval() {
-  for(unsigned int i=0;i<_dependencies.size(); i++){
-    _function->arguments[i].val->setAsProxy((*_dependencies[i])());
+  for(unsigned int i=0;i<_directDependencies.size(); i++){
+    _function->arguments[i].val->setAsProxy(_directDependencies[i]->get());
   }
   for (unsigned i = 0; i < _function->_functionReplaces.size(); i++) {
-    _function->_functionReplaces[i]->currentCache = functionReplaceCaches[i];
-    for (unsigned j = 0; j < functionReplaceCaches[i]->toReplace.size() ; j++){
-      _function->_functionReplaces[i]->_toReplace[j].val->setAsProxy((*functionReplaceCaches[i]->toReplace[j])._value);
+    _function->_functionReplaces[i]->currentCache = &functionReplaceCaches[i];
+    for (unsigned j = 0; j < functionReplaceCaches[i].toReplace.size() ; j++){
+      _function->_functionReplaces[i]->_toReplace[j].val->setAsProxy((*functionReplaceCaches[i].toReplace[j])._value);
     }
   }
   _function->call(&_cacheMap, _value);
+  _valid = true;
 }
 
-//dataCacheMap members
+
+
 dataCacheDouble &dataCacheMap::get(const function *f, dataCacheDouble *caller) {
+  // do I have a cache for this function ?
   dataCacheDouble *&r = _cacheDoubleMap[f];
-  dataCacheMap *cParent = _parent;
-  while (cParent && r==NULL) {
-    /*std::map<const function *, dataCacheDouble *>::iterator it = cParent->_cacheDoubleMap.find(f);
-    if (it != cParent->_cacheDoubleMap.end()) {*/
-     // r = it->second;
-     r = &_parent->get(f,caller);
-      for (std::set<dataCacheDouble*>::iterator dep = r->_iDependOn.begin(); dep != r->_iDependOn.end(); dep++) {
-        //if (_cacheDoubleMap.find((*dep)->_function) != _cacheDoubleMap.end()) {
-        if (&get((*dep)->_function) != *dep) {
-          r = NULL;
+  // can I use the cache of my parent ?
+  if(_parent && r==NULL) {
+    bool okFromParent = true;
+    for (std::set<function::dependency>::iterator it = f->dependencies.begin(); it != f->dependencies.end(); it++) {
+      if (it->iMap > _parent->_secondaryCaches.size())
+        okFromParent=false;
+      dataCacheMap *m = getSecondaryCache(it->iMap);
+      if (m->_cacheDoubleMap.find(it->f) != m->_cacheDoubleMap.end()) {
+        okFromParent = false;
         break;
       }
     }
-    //}
-    cParent = cParent->_parent;
+    if (okFromParent)
+      r = &_parent->get (f,caller);
   }
-  if (r==NULL)
+  // no cache found, create a new one
+  if (r==NULL) {
     r = new dataCacheDouble (this, (function*)(f));
-  if (caller)
-    r->addMeAsDependencyOf(caller);
-  return *r;
+    r->_directDependencies.resize (f->arguments.size());
+    for (unsigned int i = 0; i < f->arguments.size(); i++) {
+      r->_directDependencies[i] = &getSecondaryCache(f->arguments[i].iMap)->get (f->arguments[i].f, r);
+    }
+    for (unsigned i = 0; i < f->_functionReplaces.size(); i++) {
+      functionReplaceCache replaceCache;
+      functionReplace *replace = f->_functionReplaces[i];
+      dataCacheMap *rMap = newChild();
+      for (unsigned i = 0; i < _secondaryCaches.size(); i ++)
+        rMap->addSecondaryCache (getSecondaryCache(i+1)->newChild());
+      for (int i = 0; i < replace->_nChildren; i ++)
+        rMap->addSecondaryCache (rMap->newChild());
+      for (std::vector<function::argument>::iterator it = replace->_toReplace.begin(); it!= replace->_toReplace.end(); it++ ) {
+        dataCacheMap *m = rMap->getSecondaryCache(it->iMap);
+        dataCacheDouble *s = new dataCacheDouble(m, (function*)it->f);
+        m->_cacheDoubleMap[it->f] = s;
+        replaceCache.toReplace.push_back(s);
+      }
+      for (std::vector<function::argument>::iterator it = replace->_toCompute.begin(); it!= replace->_toCompute.end(); it++ ) {
+        replaceCache.toCompute.push_back (&rMap->getSecondaryCache(it->iMap)->get(it->f, r));
+      }
+      replaceCache.map = rMap;
+      r->functionReplaceCaches.push_back (replaceCache); 
+    }
+    ((function*)f)->registerInDataCacheMap(this, r);
   }
 
-dataCacheDouble &dataCacheMap::substitute(const function *f) 
-{
-  dataCacheDouble *&r= _cacheDoubleMap[f];
-  r = new dataCacheDouble(this, (function*)(f));
+  // update the dependency tree
+  if (caller) {
+    r->_dependOnMe.insert(caller);
+    caller->_iDependOn.insert(r);
+    for(std::set<dataCacheDouble*>::iterator it = r->_iDependOn.begin(); it != r->_iDependOn.end(); it++) {
+      (*it)->_dependOnMe.insert(caller);
+      caller->_iDependOn.insert(*it);
+    }
+  }
   return *r;
 }
 
@@ -124,8 +140,50 @@ dataCacheMap::~dataCacheMap()
       it!=_allDataCaches.end(); it++) {
     delete *it;
   }
+  for (std::list<dataCacheMap*>::iterator it = _children.begin(); it != _children.end(); it++) {
+    delete *it;
+  }
 }
 
+void functionReplace::replace(fullMatrix<double> &v, const function *f, int iMap) {
+  _replaced.insert(function::dependency(iMap,f));
+  _toReplace.push_back(function::argument(v, iMap, f));
+}
+
+void functionReplace::get(fullMatrix<double> &v, const function *f, int iMap) {
+  bool allDepFromParent = true;
+  for (std::set<function::dependency>::iterator itDep = f->dependencies.begin(); itDep != f->dependencies.end(); itDep++){
+    bool depFromParent = (_replaced.count(*itDep)==0);
+    for (std::set<function::dependency>::iterator itDep2 = itDep->f->dependencies.begin(); itDep2 != itDep->f->dependencies.end() && depFromParent; itDep2++)
+      depFromParent &= (_replaced.count(*itDep2)==0);
+    if(depFromParent)
+      _master->dependencies.insert(*itDep);
+    else
+      allDepFromParent = false;
+  }
+  function::dependency asDep(iMap, f);
+  if (allDepFromParent && _replaced.count(asDep)==0)
+    _master->dependencies.insert(asDep);
+    
+  _toCompute.push_back(function::argument(v, iMap, f));
+}
+void functionReplace::addChild() {
+  _nChildren++;
+}
+functionReplace::functionReplace(){
+  _nChildren=0;
+}
+
+void functionReplace::compute(){
+  for (unsigned i = 0; i < _toReplace.size(); i++){
+    currentCache->toReplace[i]->set();
+  }
+  for (unsigned i = 0; i < _toCompute.size(); i++) {
+    _toCompute[i].val->setAsProxy(currentCache->toCompute[i]->get());
+  }
+};
+
+
 // now some example of functions
 
 // constant values copied over each line
@@ -407,8 +465,7 @@ class functionLua : public function {
 void dataCacheMap::setNbEvaluationPoints(int nbEvaluationPoints) {
   _nbEvaluationPoints = nbEvaluationPoints;
   for(std::set<dataCacheDouble*>::iterator it = _allDataCaches.begin(); it!= _allDataCaches.end(); it++){
-    (*it)->resize();
-    (*it)->_valid = false;
+    (*it)->resize(nbEvaluationPoints);
   }
     for(std::list<dataCacheMap*>::iterator it = _children.begin(); it!= _children.end(); it++) {
       (*it)->setNbEvaluationPoints(nbEvaluationPoints);
@@ -562,44 +619,6 @@ void function::registerBindings(binding *b){
 #endif
 }
 
-void functionReplace::replace(fullMatrix<double> &v, const function *f, int iMap) {
-  _toReplace.push_back(function::argument(v, iMap, f));
-}
-
-void functionReplace::get(fullMatrix<double> &v, const function *f, int iMap) {
-  _toCompute.push_back(function::argument(v, iMap, f));
-}
-
-void functionReplace::compute(){
-  for (unsigned i = 0; i < _toReplace.size(); i++){
-    currentCache->toReplace[i]->set();
-  }
-  for (unsigned i = 0; i < _toCompute.size(); i++) {
-    currentCache->toCompute[i]->_valid =false;
-    _toCompute[i].val->setAsProxy((*currentCache->toCompute[i])());
-  }
-};
-
-
-functionReplaceCache::functionReplaceCache(dataCacheMap *m, functionReplace *rep, dataCacheDouble *from) {
-  map = m->newChild();
-  for (unsigned i = 0; i < m->_secondaryCaches.size(); i ++) {
-    map->addSecondaryCache (m->getSecondaryCache(i+1)->newChild());
-  }
-  for (unsigned i = 0; i < rep->_toReplace.size(); i++) {
-    toReplace.push_back (&map->getSecondaryCache(rep->_toReplace[i].iMap)->substitute(rep->_toReplace[i].f));
-  }
-  for (unsigned i = 0; i < rep->_toCompute.size(); i++) {
-    dataCacheMap *m2 = map->getSecondaryCache(rep->_toCompute[i].iMap);
-    toCompute.push_back (&m2->get(rep->_toCompute[i].f, from));
-  }
-}
-functionReplaceCache::~functionReplaceCache() {
-  for (unsigned i = 0; i< map->_secondaryCaches.size(); i++) {
-    delete map->_secondaryCaches[i];
-  }
-  delete map;
-}
 functionConstant *function::_timeFunction = NULL;
 functionConstant *function::getTime() {
   if (! _timeFunction)

Solver/function.h

@@ -21,56 +21,46 @@ class functionReplaceCache;
 // An abstract interface to functions 
 // more explanation at the head of this file
 class function {
-  static functionConstant *_timeFunction;
-  static functionConstant *_dtFunction;  
-  public :
-  class argument {
-    //iMap is the id of the dataCacheMap, e.g. on interfaces
   public :
-    int iMap;
-    const function *f;
-    fullMatrix<double> *val;
-    argument(fullMatrix<double> &v, int iMap_, const function *f_){
-      val = &v;
-      iMap = iMap_;
-      f = f_;
+  int _nbCol;
+  bool _invalidatedOnElement;
+  std::vector<functionReplace*> _functionReplaces;
+  class dependency {
+    public : int iMap; const function *f;
+    dependency(int iMap_, const function *f_){iMap = iMap_; f = f_; }
+    bool operator < (const dependency &d) const {
+      return (d.iMap <iMap || d.f < f);
     }
   };
-  std::vector<functionReplace*> _functionReplaces;
-  class substitutedFunction {
+  void printDep() {
+    for(std::set<dependency>::iterator it = dependencies.begin(); it != dependencies.end(); it++)
+      printf("%i %p\n", it->iMap, it->f);
+  }
+  class argument {
+    //iMap is the id of the dataCacheMap, e.g. on interfaces
     public:
-    int iMap;
-    const function *f0, *f1; // f1 replaces f0
+    int iMap; const function *f; fullMatrix<double> *val;
+    argument(fullMatrix<double> &v, int iMap_, const function *f_){ val = &v; iMap = iMap_; f = f_; }
   };
-  int _nbCol;
-  bool _invalidatedOnElement;
-  std::vector<int> _childrenCache;
-  std::vector<substitutedFunction> _substitutedFunctions;
-  virtual void call (dataCacheMap *m, fullMatrix<double> &res)=0;
-  virtual void registerInDataCacheMap(dataCacheMap *m, dataCacheDouble *d) {}
   std::vector<argument> arguments;
-  const void setArgument(fullMatrix<double> &v, const function *f, int iMap = 0) {
-    if(f==NULL)
-      throw;
-    arguments.push_back(argument(v, iMap, f));
-  }
-  void addChildDataCacheMap(int parent) {
-    _childrenCache.push_back(parent);
-  }
-  void substituteFunction( int iMap, const function *f0, const function *f1) {
-    substitutedFunction s;
-    s.iMap= iMap;
-    s.f0 = f0;
-    s.f1 = f1;
-    _substitutedFunctions.push_back(s);
-  }
-  virtual ~function();
-  static void registerBindings(binding *b);
+  std::set<dependency> dependencies;
+  void setArgument(fullMatrix<double> &v, const function *f, int iMap = 0);
+  void addFunctionReplace(functionReplace &fr);
+
   function(int nbCol, bool invalidatedOnElement = true);
-  inline int getNbCol()const {return _nbCol;}
+  virtual ~function();
+  virtual void call (dataCacheMap *m, fullMatrix<double> &res)=0;
+  virtual void registerInDataCacheMap(dataCacheMap *m, dataCacheDouble *d) {}
+
   inline bool isInvalitedOnElement() { return _invalidatedOnElement;}
-  void addFunctionReplace(functionReplace &fr);
+  inline int getNbCol()const {return _nbCol;}
+
+  static void registerBindings(binding *b);
 
+  private:
+  static functionConstant *_timeFunction;
+  static functionConstant *_dtFunction;  
+  public:
   static function *getSolution();
   static function *getSolutionGradient();
   static function *getParametricCoordinates();
@@ -84,41 +74,24 @@ class function {
 // then the size of the matrix is automatically adjusted
 class dataCacheDouble {
   friend class dataCacheMap;
-  public:
+  friend class dgDataCacheMap;
+  friend class functionReplace;
+  std::vector<dataCacheDouble*> _directDependencies;
+  function *_function;
+  dataCacheMap &_cacheMap;
+  std::vector<functionReplaceCache> functionReplaceCaches;
+  protected:
   // pointers to all of the dataCache depending on me
   std::set<dataCacheDouble*> _dependOnMe;
   std::set<dataCacheDouble*> _iDependOn;
-protected :
-  // invalidates all the cached data that depends on me
-  inline void _invalidateDependencies()
-  {
-    // if this is too slow we can keep a C array cache of the _dependOnMe set
-    for(std::set<dataCacheDouble*>::iterator it = _dependOnMe.begin();
-      it!=_dependOnMe.end(); it++)
-        (*it)->_valid=false;
-  }
-public :
+  fullMatrix<double> _value;
   bool _valid;
-  // dataCacheMap is the only one supposed to call this
-  void addMeAsDependencyOf (dataCacheDouble *newDep);
-  inline bool somethingDependOnMe() {
-    return !_dependOnMe.empty();
-  }
-  inline bool doIDependOn(dataCacheDouble &other) {
-    return (_iDependOn.find(&other)!=_iDependOn.end());
-  }
-  std::vector<dataCacheDouble*> _dependencies;
-  std::vector<std::pair<dataCacheDouble*, dataCacheDouble*> > _substitutions;
-
+  dataCacheDouble(dataCacheMap *,function *f);
 
-  int _nRowByPoint;
-  function *_function;
-  dataCacheMap &_cacheMap;
   // do the actual computation and put the result into _value
-  // still virtual because it is overrided by conservation law terms, as soon as conservation law terms will be regular functions, we will remove this
-  virtual void _eval();
+  void _eval();
+  void resize(int nrow);
  public:
-  fullMatrix<double> _value;
   //set the value (without computing it by _eval) and invalidate the dependencies
   // this function is needed to be able to pass the _value to functions like gemm or mult
   // but you cannot keep the reference to the _value, you should always use the set function 
@@ -126,41 +99,33 @@ public :
   // take care if you use this to set a proxy you must ensure that the value pointed to are not modified
   // without further call to set because the dependencies won't be invalidate
   inline fullMatrix<double> &set() {
-    if(_valid)
-      _invalidateDependencies();
-    _valid = true;
-    return _value;
+    if(_valid) {
+      for(std::set<dataCacheDouble*>::iterator it = _dependOnMe.begin();
+          it!=_dependOnMe.end(); it++)
+        (*it)->_valid=false;
     }
-  std::vector<functionReplaceCache*> functionReplaceCaches;
-
-  inline const double &operator () (int i, int j)
-  {
-    if(!_valid) {
-      _eval();
     _valid = true;
-    }
-    return _value(i,j);
+    return _value;
   }
   //access _value and compute it if necessary
-  inline const fullMatrix<double> &operator () ()
-  {
-    if(!_valid) {
+  inline const fullMatrix<double> &get() {
+    if(!_valid)
       _eval();
-      _valid=true;
-    }
     return _value;
   }
-  void resize();
-  dataCacheDouble(dataCacheMap *,function *f);
-  dataCacheDouble(dataCacheMap &m, int nRowByPoint, int nbCol);
-  virtual ~dataCacheDouble();
+  // dataCacheMap is the only one supposed to call this
+  inline bool somethingDependOnMe() {
+    return !_dependOnMe.empty();
+  }
+  inline bool doIDependOn(dataCacheDouble &other) {
+    return (_iDependOn.find(&other)!=_iDependOn.end());
+  }
 };
 
 
 class dgDataCacheMap;
 // more explanation at the head of this file
 class dataCacheMap {
-  friend class dataCacheDouble;
  public:
   dataCacheMap  *_parent;
   std::list<dataCacheMap*> _children;
@@ -177,10 +142,6 @@ class dataCacheMap {
     if(invalidatedOnElement)
       _toInvalidateOnElement.insert(data);
   }
-  void printList() {
-    for(std::set<dataCacheDouble*>::iterator it = _toInvalidateOnElement.begin(); it!= _toInvalidateOnElement.end(); it++)
-      printf("%p\n",*it);
-  }
   virtual dgDataCacheMap *asDgDataCacheMap() {
     Msg::Error("I'm not a dgDataCacheMap\n");
     return NULL;
@@ -194,7 +155,6 @@ class dataCacheMap {
     _secondaryCaches.push_back(s);
   }
   dataCacheDouble &get(const function *f, dataCacheDouble *caller=0);
-  dataCacheDouble &substitute(const function *f);
   virtual void setElement(MElement *element) {
     _element=element;
     for(std::set<dataCacheDouble*>::iterator it = _toInvalidateOnElement.begin(); it!= _toInvalidateOnElement.end(); it++) {
@@ -221,27 +181,22 @@ class dataCacheMap {
   inline int getNbEvaluationPoints(){return _nbEvaluationPoints;}
   ~dataCacheMap();
 };
-
 class functionReplace {
-  friend class functionReplaceCache;
+  friend class dataCacheMap;
   friend class dataCacheDouble;
-  protected:
+  public :
+  function *_master;
+  int _nChildren;
   functionReplaceCache *currentCache;
+  std::set <function::dependency> _replaced;
+  std::set <function::dependency> _fromParent;
   std::vector <function::argument> _toReplace;
   std::vector <function::argument> _toCompute;
-  public :
   void get(fullMatrix<double> &v, const function *, int iMap = 0);
   void replace(fullMatrix<double> &v, const function *, int iMap = 0);
   void compute ();
-};
-
-class functionReplaceCache {
-  public :
-  dataCacheMap *map;
-  std::vector <dataCacheDouble*> toReplace;
-  std::vector <dataCacheDouble*> toCompute;
-  functionReplaceCache(dataCacheMap *m, functionReplace *rep, dataCacheDouble *from);
-  ~functionReplaceCache();
+  functionReplace();
+  void addChild();
 };
 
 

Solver/functionDerivator.cpp

@@ -2,17 +2,18 @@
 #include "function.h"
 
 void functionDerivator::compute() {
-  _xRef = _x();
-  _fRef = _f();
-  _dfdx = fullMatrix<double>(_f().size1(),_f().size2()*_x().size2());
+  _xRef = _x.get();
+  _fRef = _f.get();
+  _dfdx = fullMatrix<double>(_fRef.size1(),_fRef.size2()*_xRef.size2());
   for (int j=0;j<_xRef.size2();j++) {
     _xDx = _xRef;
     for (int i=0;i<_fRef.size1();i++)
       _xDx(i,j) += _epsilon;
     _x.set()=_xDx;
+    const fullMatrix<double> &f = _f.get();
     for (int i=0; i<_fRef.size1(); i++) 
       for (int k=0; k<_fRef.size2(); k++) 
-        _dfdx(i, k*_xRef.size2()+j) = (_f(i,k)-_fRef(i,k))/_epsilon;
+        _dfdx(i, k*_xRef.size2()+j) = (f(i,k)-_fRef(i,k))/_epsilon;
   }
   _x.set()=_xRef;
 };