Reuse factors of linear equation systems in fmi2GetDirectionalDerivative

Belonging to [master]:
  - #2260

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);

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)
 {

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;

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<ISimObjects> 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

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;
 }
 

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 {

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

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,48 +131,66 @@ 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
 
   _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;