diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp
index 1f9a6d405..b881f0595 100644
--- a/opm/simulators/wells/MultisegmentWellEquations.cpp
+++ b/opm/simulators/wells/MultisegmentWellEquations.cpp
@@ -161,6 +161,13 @@ void MultisegmentWellEquations<Scalar,numWellEq,numEq>::createSolver()
 #endif
 }
 
+template<class Scalar, int numWellEq, int numEq>
+typename MultisegmentWellEquations<Scalar,numWellEq,numEq>::BVectorWell
+MultisegmentWellEquations<Scalar,numWellEq,numEq>::solve() const
+{
+    return mswellhelpers::applyUMFPack(*duneDSolver_, resWell_);
+}
+
 #define INSTANCE(numWellEq, numEq) \
 template class MultisegmentWellEquations<double,numWellEq,numEq>;
 
diff --git a/opm/simulators/wells/MultisegmentWellEquations.hpp b/opm/simulators/wells/MultisegmentWellEquations.hpp
index 1818720d2..2025d6b1c 100644
--- a/opm/simulators/wells/MultisegmentWellEquations.hpp
+++ b/opm/simulators/wells/MultisegmentWellEquations.hpp
@@ -83,6 +83,9 @@ public:
     //! \brief Compute the LU-decomposition of D matrix.
     void createSolver();
 
+    //! \brief Apply inverted D matrix to residual and return result.
+    BVectorWell solve() const;
+
     // two off-diagonal matrices
     OffDiagMatWell duneB_;
     OffDiagMatWell duneC_;
diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp
index 25aa96b63..f03affb6a 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -519,8 +519,7 @@ namespace Opm
 
         // We assemble the well equations, then we check the convergence,
         // which is why we do not put the assembleWellEq here.
-        const BVectorWell dx_well = mswellhelpers::applyUMFPack(*this->linSys_.duneDSolver_,
-                                                                this->linSys_.resWell_);
+        const BVectorWell dx_well = this->linSys_.solve();
 
         updateWellState(dx_well, well_state, deferred_logger);
     }
@@ -1495,8 +1494,7 @@ namespace Opm
 
             assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
 
-            const BVectorWell dx_well = mswellhelpers::applyUMFPack(*this->linSys_.duneDSolver_,
-                                                                    this->linSys_.resWell_);
+            const BVectorWell dx_well = this->linSys_.solve();
 
             if (it > this->param_.strict_inner_iter_wells_) {
                 relax_convergence = true;