From 514196213912d7819f965562012f053f7085d316 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=BCdiger=20Franke?= Date: Wed, 7 Mar 2018 20:40:55 +0100 Subject: [PATCH] Reuse factors of linear equation systems in fmi2GetDirectionalDerivative Belonging to [master]: - OpenModelica/OMCompiler#2260 --- Compiler/Template/CodegenCpp.tpl | 6 ++ .../System/SystemDefaultImplementation.cpp | 14 +++++ .../cpp/Include/Core/System/ILinearAlgLoop.h | 1 + .../Core/System/SystemDefaultImplementation.h | 24 +++++++ .../cpp/Include/FMU2/FMU2Wrapper.cpp | 5 ++ .../cpp/Include/FMU2/FMU2Wrapper.h | 3 +- .../cpp/Include/Solver/Dgesv/DgesvSolver.h | 4 +- .../cpp/Solver/Dgesv/DgesvSolver.cpp | 63 ++++++++++++------- 8 files changed, 96 insertions(+), 24 deletions(-) diff --git a/Compiler/Template/CodegenCpp.tpl b/Compiler/Template/CodegenCpp.tpl index 5ddfc90c89..3b126a2873 100644 --- a/Compiler/Template/CodegenCpp.tpl +++ b/Compiler/Template/CodegenCpp.tpl @@ -3818,6 +3818,11 @@ case SIMCODE(modelInfo = MODELINFO(__)) then } + bool <%modelname%>Algloop<%ls.index%>::getFreeVariablesLock() + { + return _system->getFreeVariablesLock(); + } + bool <%modelname%>Algloop<%ls.index%>::getUseSparseFormat() { return LinearAlgLoopDefaultImplementation::getUseSparseFormat(); @@ -7080,6 +7085,7 @@ case SIMCODE(modelInfo = MODELINFO(__)) then virtual const matrix_t& getAMatrix() ; virtual sparsematrix_t& getSparseAMatrix() ; + bool getFreeVariablesLock(); bool getUseSparseFormat(); void setUseSparseFormat(bool value); diff --git a/SimulationRuntime/cpp/Core/System/SystemDefaultImplementation.cpp b/SimulationRuntime/cpp/Core/System/SystemDefaultImplementation.cpp index 0af21cc817..3c80498111 100644 --- a/SimulationRuntime/cpp/Core/System/SystemDefaultImplementation.cpp +++ b/SimulationRuntime/cpp/Core/System/SystemDefaultImplementation.cpp @@ -71,6 +71,7 @@ SystemDefaultImplementation::SystemDefaultImplementation(IGlobalSettings *global , _conditions0(NULL) , _event_system(NULL) , _modelName(modelName) + , _freeVariablesLock(false) { } @@ -111,6 +112,7 @@ SystemDefaultImplementation::SystemDefaultImplementation(SystemDefaultImplementa , _conditions0(NULL) , _event_system(NULL) , _modelName(instance.getModelName()) + , _freeVariablesLock(false) { } @@ -284,6 +286,18 @@ double SystemDefaultImplementation::getTime() return _simTime; } +/// Set status of independent variables +void SystemDefaultImplementation::setFreeVariablesLock(bool freeVariablesLock) +{ + _freeVariablesLock = freeVariablesLock; +} + +// Get status of independent variables +bool SystemDefaultImplementation::getFreeVariablesLock() +{ + return _freeVariablesLock; +} + /// getter for variables of different types void SystemDefaultImplementation::getBoolean(bool* z) { diff --git a/SimulationRuntime/cpp/Include/Core/System/ILinearAlgLoop.h b/SimulationRuntime/cpp/Include/Core/System/ILinearAlgLoop.h index 472fcaac92..a0b005b3e0 100644 --- a/SimulationRuntime/cpp/Include/Core/System/ILinearAlgLoop.h +++ b/SimulationRuntime/cpp/Include/Core/System/ILinearAlgLoop.h @@ -61,6 +61,7 @@ class ILinearAlgLoop virtual sparsematrix_t& getSparseAMatrix() = 0; virtual bool isLinearTearing() = 0; virtual bool isConsistent() = 0; + virtual bool getFreeVariablesLock() = 0; virtual bool getUseSparseFormat() = 0; virtual void setUseSparseFormat(bool value) = 0; virtual float queryDensity() = 0; diff --git a/SimulationRuntime/cpp/Include/Core/System/SystemDefaultImplementation.h b/SimulationRuntime/cpp/Include/Core/System/SystemDefaultImplementation.h index c1a0f46ce5..9ed5ac9fc2 100644 --- a/SimulationRuntime/cpp/Include/Core/System/SystemDefaultImplementation.h +++ b/SimulationRuntime/cpp/Include/Core/System/SystemDefaultImplementation.h @@ -124,6 +124,11 @@ class BOOST_EXTENSION_SYSTEM_DECL SystemDefaultImplementation void setTime(const double& t); double getTime(); + /// Set modification status of independent variables + // (exploited by linear equation systems in Jacobians) + void setFreeVariablesLock(bool freeVariablesLock); + bool getFreeVariablesLock(); + IGlobalSettings* getGlobalSettings(); shared_ptr getSimObjects() const; @@ -220,5 +225,24 @@ class BOOST_EXTENSION_SYSTEM_DECL SystemDefaultImplementation bool _sparse; bool _useAnalyticalJacobian; + + bool _freeVariablesLock; ///< modification status of independent variables +}; + +/// Mark free variables unchanged. +/// Automatically release lock upon destruction, covering exceptions as well. +class SystemLockFreeVariables { + public: + SystemLockFreeVariables(SystemDefaultImplementation *system, bool lock = true) { + system->setFreeVariablesLock(lock); + _system = system; + } + + ~SystemLockFreeVariables() { + _system->setFreeVariablesLock(false); + } + + private: + SystemDefaultImplementation *_system; }; /** @} */ // end of coreSystem diff --git a/SimulationRuntime/cpp/Include/FMU2/FMU2Wrapper.cpp b/SimulationRuntime/cpp/Include/FMU2/FMU2Wrapper.cpp index 48f62081bb..3306dbdddb 100644 --- a/SimulationRuntime/cpp/Include/FMU2/FMU2Wrapper.cpp +++ b/SimulationRuntime/cpp/Include/FMU2/FMU2Wrapper.cpp @@ -125,6 +125,7 @@ FMU2Wrapper::FMU2Wrapper(fmi2String instanceName, fmi2String GUID, _model = createSystemFMU(&_global_settings); _model->initializeMemory(); _model->initializeFreeVariables(); + _need_update = true; _string_buffer.resize(_model->getDimString()); _clockTick = new bool[_model->getDimClock()]; _clockSubactive = new bool[_model->getDimClock()]; @@ -232,6 +233,7 @@ fmi2Status FMU2Wrapper::terminate() fmi2Status FMU2Wrapper::reset() { _model->initializeFreeVariables(); + _need_update = true; return fmi2OK; } @@ -243,6 +245,7 @@ void FMU2Wrapper::updateModel() } _model->evaluateAll(); // derivatives and algebraic variables _need_update = false; + _needJacUpdate = true; } fmi2Status FMU2Wrapper::setTime(fmi2Real time) @@ -486,9 +489,11 @@ fmi2Status FMU2Wrapper::getDirectionalDerivative(const fmi2ValueReference vrUnkn { if (_need_update) updateModel(); + SystemLockFreeVariables slfv(_model, !_needJacUpdate); _model->getDirectionalDerivative(vrUnknown, nUnknown, vrKnown, nKnown, dvKnown, dvUnknown); + _needJacUpdate = false; return fmi2OK; } diff --git a/SimulationRuntime/cpp/Include/FMU2/FMU2Wrapper.h b/SimulationRuntime/cpp/Include/FMU2/FMU2Wrapper.h index b08ad554fd..23496ffa0d 100644 --- a/SimulationRuntime/cpp/Include/FMU2/FMU2Wrapper.h +++ b/SimulationRuntime/cpp/Include/FMU2/FMU2Wrapper.h @@ -178,7 +178,8 @@ class FMU2Wrapper bool *_clockTick; bool *_clockSubactive; int _nclockTick; - double _need_update; + bool _need_update; + bool _needJacUpdate; void updateModel(); typedef enum { diff --git a/SimulationRuntime/cpp/Include/Solver/Dgesv/DgesvSolver.h b/SimulationRuntime/cpp/Include/Solver/Dgesv/DgesvSolver.h index 48c079a832..e9a351cd35 100644 --- a/SimulationRuntime/cpp/Include/Solver/Dgesv/DgesvSolver.h +++ b/SimulationRuntime/cpp/Include/Solver/Dgesv/DgesvSolver.h @@ -38,7 +38,9 @@ class DgesvSolver : public IAlgLoopSolver *_jHelp; ///< Pivot indices for LAPACK routines bool - _firstCall; ///< Temp - Denotes the first call to the solver, init() is called + _firstCall, ///< Temp - Denotes the first call to the solver, init() is called + _hasDgesvFactors, ///< =true if previous dgesv was called + _hasDgetc2Factors; ///< =true if previous dgetc2 was called const char* *_yNames; ///< Names of variables diff --git a/SimulationRuntime/cpp/Solver/Dgesv/DgesvSolver.cpp b/SimulationRuntime/cpp/Solver/Dgesv/DgesvSolver.cpp index f557617e24..3f85a30caa 100644 --- a/SimulationRuntime/cpp/Solver/Dgesv/DgesvSolver.cpp +++ b/SimulationRuntime/cpp/Solver/Dgesv/DgesvSolver.cpp @@ -18,7 +18,6 @@ DgesvSolver::DgesvSolver(ILinearAlgLoop* algLoop, ILinSolverSettings* settings) : _algLoop (algLoop) , _dimSys (0) - , _yNames (NULL) , _yNominal (NULL) , _y (NULL) @@ -32,6 +31,8 @@ DgesvSolver::DgesvSolver(ILinearAlgLoop* algLoop, ILinSolverSettings* settings) , _zeroVec (NULL) , _iterationStatus (CONTINUE) , _firstCall (true) + , _hasDgesvFactors (false) + , _hasDgetc2Factors (false) , _fNominal (NULL) { } @@ -130,7 +131,8 @@ void DgesvSolver::solve() //use lapack long int dimRHS = 1; // Dimension of right hand side of linear system (=_b) - long int irtrn = 0; // Return-flag of Fortran code + long int info = 0; // Return flag of dgesv + double scale = 0.0; // Scale factor of dgesc2 if (_algLoop->isLinearTearing()) _algLoop->setReal(_zeroVec); //if the system is linear Tearing it means that the system is of the form Ax-b=0, so plugging in x=0 yields -b for the left hand side @@ -138,40 +140,57 @@ void DgesvSolver::solve() _algLoop->evaluate(); _algLoop->getb(_b); - const matrix_t& A = _algLoop->getAMatrix(); - const double* Atemp = A.data().begin(); + if (!_algLoop->getFreeVariablesLock()) { + const matrix_t& A = _algLoop->getAMatrix(); + const double* Atemp = A.data().begin(); - memcpy(_A, Atemp, _dimSys*_dimSys*sizeof(double)); + memcpy(_A, Atemp, _dimSys*_dimSys*sizeof(double)); + _hasDgesvFactors = false; + _hasDgetc2Factors = false; - std::fill(_fNominal, _fNominal + _dimSys, 1e-6); - for (int j = 0, idx = 0; j < _dimSys; j++) - for (int i = 0; i < _dimSys; i++, idx++) - _fNominal[i] = std::max(std::abs(Atemp[idx]), _fNominal[i]); + std::fill(_fNominal, _fNominal + _dimSys, 1e-6); + for (int j = 0, idx = 0; j < _dimSys; j++) + for (int i = 0; i < _dimSys; i++, idx++) + _fNominal[i] = std::max(std::abs(Atemp[idx]), _fNominal[i]); - LOGGER_WRITE_VECTOR("fNominal", _fNominal, _dimSys, LC_LS, LL_DEBUG); + LOGGER_WRITE_VECTOR("fNominal", _fNominal, _dimSys, LC_LS, LL_DEBUG); - for (int j = 0, idx = 0; j < _dimSys; j++) - for (int i = 0; i < _dimSys; i++, idx++) - _A[idx] /= _fNominal[i]; + for (int j = 0, idx = 0; j < _dimSys; j++) + for (int i = 0; i < _dimSys; i++, idx++) + _A[idx] /= _fNominal[i]; + } for (int i = 0; i < _dimSys; i++) _b[i] /= _fNominal[i]; - dgesv_(&_dimSys, &dimRHS, _A, &_dimSys, _iHelp, _b, &_dimSys, &irtrn); + if (!_hasDgesvFactors && !_hasDgetc2Factors) { + dgesv_(&_dimSys, &dimRHS, _A, &_dimSys, _iHelp, _b, &_dimSys, &info); + _hasDgesvFactors = true; + } + else if (_hasDgesvFactors) { + // solve using previously obtained dgesv factors + char trans = 'N'; + dgetrs_(&trans, &_dimSys, &dimRHS, _A, &_dimSys, _iHelp, _b, &_dimSys, &info); + } + else { + // solve using previously obtained dgetc2 factors + dgesc2_(&_dimSys, _A, &_dimSys, _b, _iHelp, _jHelp, &scale); + info = 0; + } - if (irtrn > 0) { - long int info = 0; - double scale = 0.0; - dgetc2_(&_dimSys, _A, &_dimSys, _iHelp, _jHelp, &info); + if (info > 0) { + long int info2 = 0; + dgetc2_(&_dimSys, _A, &_dimSys, _iHelp, _jHelp, &info2); dgesc2_(&_dimSys, _A, &_dimSys, _b, _iHelp, _jHelp, &scale); - LOGGER_WRITE("DgesvSolver: using complete pivoting (dgesv info: " + to_string(irtrn) + ", dgesc2 scale: " + to_string(scale) + ")", LC_LS, LL_DEBUG); + _hasDgetc2Factors = true; + LOGGER_WRITE("DgesvSolver: using complete pivoting (dgesv info: " + to_string(info) + ", dgesc2 scale: " + to_string(scale) + ")", LC_LS, LL_DEBUG); } - else if (irtrn < 0) { + else if (info < 0) { _iterationStatus = SOLVERERROR; if (_algLoop->isLinearTearing()) - throw ModelicaSimulationError(ALGLOOP_SOLVER, "error solving linear tearing system (dgesv info: " + to_string(irtrn) + ")"); + throw ModelicaSimulationError(ALGLOOP_SOLVER, "error solving linear tearing system (dgesv info: " + to_string(info) + ")"); else - throw ModelicaSimulationError(ALGLOOP_SOLVER, "error solving linear system (dgesv info: " + to_string(irtrn) + ")"); + throw ModelicaSimulationError(ALGLOOP_SOLVER, "error solving linear system (dgesv info: " + to_string(info) + ")"); } _iterationStatus = DONE;