From c4e5102f474e43f45b6e66c023983df8d235e64e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=BCdiger=20Franke?= Date: Wed, 7 Mar 2018 22:43:28 +0100 Subject: [PATCH] Reuse factors of linear equation systems during Jacobian evaluation Belonging to [master]: - OpenModelica/OMCompiler#2261 --- Compiler/Template/CodegenCpp.tpl | 5 +- .../Solver/LinearSolver/LinearSolver.h | 4 +- .../cpp/Solver/LinearSolver/LinearSolver.cpp | 64 ++++++++++++------- 3 files changed, 48 insertions(+), 25 deletions(-) diff --git a/Compiler/Template/CodegenCpp.tpl b/Compiler/Template/CodegenCpp.tpl index 3b126a28731..e65e68fc21d 100644 --- a/Compiler/Template/CodegenCpp.tpl +++ b/Compiler/Template/CodegenCpp.tpl @@ -13052,11 +13052,12 @@ case _ then let jacvals = if stringEq(eqsCount, "0") then '' else (sparsepattern |> (index,indexes) hasindex index0 => let jaccol = ( indexes |> i_index hasindex index1 => - (match indexColumn case "1" then '_<%matrixName%>jacobian(0,<%index%>) = _<%matrixName%>jac_y(0);/*test1<%index0%>,<%index1%>*/' - else '_<%matrixName%>jacobian(<%i_index%>,<%index%>) = _<%matrixName%>jac_y(<%i_index%>);/*test2<%index0%>,<%index1%>*/' + (match indexColumn case "1" then '_<%matrixName%>jacobian(0,<%index%>) = _<%matrixName%>jac_y(0);' + else '_<%matrixName%>jacobian(<%i_index%>,<%index%>) = _<%matrixName%>jac_y(<%i_index%>);' ) ;separator="\n") << + <%if intEq(index0, 1) then 'SystemLockFreeVariables slfv(this);'%> _<%matrixName%>jac_x(<%index0%>) = 1; calc<%matrixName%>JacobianColumn(); _<%matrixName%>jac_x.clear(); diff --git a/SimulationRuntime/cpp/Include/Solver/LinearSolver/LinearSolver.h b/SimulationRuntime/cpp/Include/Solver/LinearSolver/LinearSolver.h index b0d1863ae63..7b6a4e821c7 100644 --- a/SimulationRuntime/cpp/Include/Solver/LinearSolver/LinearSolver.h +++ b/SimulationRuntime/cpp/Include/Solver/LinearSolver/LinearSolver.h @@ -42,7 +42,9 @@ class LinearSolver : public IAlgLoopSolver _dimSys; ///< Temp - Number of unknowns (=dimension of system of equations) 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 dgesv was called previously + _hasDgetc2Factors; ///< =true if dgetc2 was called previously long int *_ihelpArray, //pivot indices for lapackroutine *_jhelpArray; //pivot indices for lapackroutine diff --git a/SimulationRuntime/cpp/Solver/LinearSolver/LinearSolver.cpp b/SimulationRuntime/cpp/Solver/LinearSolver/LinearSolver.cpp index 5e730572528..8c3c4c7bee2 100644 --- a/SimulationRuntime/cpp/Solver/LinearSolver/LinearSolver.cpp +++ b/SimulationRuntime/cpp/Solver/LinearSolver/LinearSolver.cpp @@ -40,6 +40,8 @@ LinearSolver::LinearSolver(ILinearAlgLoop* algLoop, ILinSolverSettings* settings , _iterationStatus (CONTINUE) , _firstCall (true) + , _hasDgesvFactors (false) + , _hasDgetc2Factors (false) , _scale (NULL) , _generateoutput (false) , _fNominal (NULL) @@ -198,27 +200,31 @@ void LinearSolver::solve() //if !_sparse, we use LAPACK routines, otherwise we use KLU to solve the linear system if (!_sparse) { //use lapack - long int dimRHS = 1; // Dimension of right hand side of linear system (=_b) - long int irtrn = 0; // Return-flag of Fortran code - - const matrix_t& A = _algLoop->getAMatrix(); - const double* Atemp = A.data().begin(); - - memcpy(_A, Atemp, _dimSys*_dimSys*sizeof(double)); - - // scale Jacobian - 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]); + long int dimRHS = 1; // Dimension of right hand side of linear system (=_b) + long int info = 0; // Return-flag of Fortran code + + if (!_algLoop->getFreeVariablesLock()) { + const matrix_t& A = _algLoop->getAMatrix(); + const double* Atemp = A.data().begin(); + + memcpy(_A, Atemp, _dimSys*_dimSys*sizeof(double)); + _hasDgesvFactors = false; + _hasDgetc2Factors = false; + + // scale Jacobian + 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]; @@ -228,7 +234,7 @@ void LinearSolver::solve() std::cout << "We solve a linear system with coefficient matrix" << std::endl; for (int i=0; i<_dimSys; i++) { for (int j=0; j<_dimSys; j++) { - std::cout << Atemp[i+j*_dimSys] << " "; + std::cout << _A[i+j*_dimSys] << " "; } std::cout << std::endl; } @@ -239,11 +245,25 @@ void LinearSolver::solve() std::cout << std::endl; } - dgesv_(&_dimSys, &dimRHS, _A, &_dimSys, _ihelpArray, _b, &_dimSys, &irtrn); + if (!_hasDgesvFactors && !_hasDgetc2Factors) { + dgesv_(&_dimSys, &dimRHS, _A, &_dimSys, _ihelpArray, _b, &_dimSys, &info); + _hasDgesvFactors = true; + } + else if (_hasDgesvFactors) { + // solve using previously obtained dgesv factors + char trans = 'N'; + dgetrs_(&trans, &_dimSys, &dimRHS, _A, &_dimSys, _ihelpArray, _b, &_dimSys, &info); + } + else { + // solve using previously obtained dgetc2 factors + dgesc2_(&_dimSys, _A, &_dimSys, _b, _ihelpArray, _jhelpArray, _scale); + info = 0; + } - if (irtrn != 0) { - dgetc2_(&_dimSys, _A, &_dimSys, _ihelpArray, _jhelpArray, &irtrn); + if (info != 0) { + dgetc2_(&_dimSys, _A, &_dimSys, _ihelpArray, _jhelpArray, &info); dgesc2_(&_dimSys, _A, &_dimSys, _b, _ihelpArray, _jhelpArray, _scale); + _hasDgetc2Factors = true; LOGGER_WRITE("LinearSolver: Linear system singular, using perturbed system matrix.", LC_LS, LL_DEBUG); _iterationStatus = DONE; }