From 5099743b1e874ffd1fcce5b323ecd99ef3c089ad Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Wed, 19 May 2010 17:43:08 +0000
Subject: [PATCH] work on dg functions
---
Numeric/fullMatrix.cpp | 13 ++
Numeric/fullMatrix.h | 14 +-
Solver/function.cpp | 251 +++++++++++++++++++----------------
Solver/function.h | 177 +++++++++---------------
Solver/functionDerivator.cpp | 9 +-
5 files changed, 229 insertions(+), 235 deletions(-)
diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 9875746c20..d0dd089e39 100644
--- a/Numeric/fullMatrix.cpp
+++ b/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
{
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 8688c548ad..32904e38c0 100644
--- a/Numeric/fullMatrix.h
+++ b/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,14 +152,16 @@ class fullMatrix
if (_own_data && _data) delete[] _data;
_data = new scalar[_r * _c];
_own_data = true;
- scale(0.);
+ if(resetValue)
+ scale(0.);
return true;
}
else{
_r = r;
_c = c;
}
- scale(0.);
+ if(resetValue)
+ scale(0.);
return false; // no reallocation
}
void setAsProxy(const fullMatrix<scalar> &original)
@@ -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;
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 0e67868319..5c603d903f 100644
--- a/Solver/function.cpp
+++ b/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 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);
+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));
+ }
+ 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;
- break;
- }
+ // 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);
+ }
+ // no cache found, create a new one
+ if (r==NULL) {
+ r = new dataCacheDouble (this, (function*)(f));
+ 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);
}
- if (r==NULL)
- r = new dataCacheDouble(this, (function*)(f));
- if (caller)
- r->addMeAsDependencyOf(caller);
- return *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)
diff --git a/Solver/function.h b/Solver/function.h
index dfac032bc4..09d433ef4a 100644
--- a/Solver/function.h
+++ b/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();
+ if(_valid) {
+ for(std::set<dataCacheDouble*>::iterator it = _dependOnMe.begin();
+ it!=_dependOnMe.end(); it++)
+ (*it)->_valid=false;
+ }
_valid = true;
return _value;
}
- std::vector<functionReplaceCache*> functionReplaceCaches;
-
- inline const double &operator () (int i, int j)
- {
- if(!_valid) {
- _eval();
- _valid = true;
- }
- return _value(i,j);
- }
//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,13 +155,12 @@ 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++) {
(*it)->_valid=false;
}
- for(std::list<dataCacheMap*>::iterator it = _children.begin(); it!= _children.end(); it++) {
+ for(std::list<dataCacheMap*>::iterator it = _children.begin(); it != _children.end(); it++) {
(*it)->setElement(element);
}
}
@@ -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();
};
diff --git a/Solver/functionDerivator.cpp b/Solver/functionDerivator.cpp
index ea9f7d6150..79e99bccce 100644
--- a/Solver/functionDerivator.cpp
+++ b/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;
};
--
GitLab