diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp
index 049d8f948..26d6a33ca 100644
--- a/opm/simulators/wells/MultisegmentWellEquations.cpp
+++ b/opm/simulators/wells/MultisegmentWellEquations.cpp
@@ -130,6 +130,16 @@ apply(const BVector& x, BVector& Ax) const
     duneC_.mmtv(invDBx,Ax);
 }
 
+template<class Scalar, int numWellEq, int numEq>
+void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
+apply(BVector& r) const
+{
+    // invDrw_ = duneD^-1 * resWell_
+    const BVectorWell invDrw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_);
+    // r = r - duneC_^T * invDrw
+    duneC_.mmtv(invDrw, r);
+}
+
 #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 f2c9e73d7..324cc82c0 100644
--- a/opm/simulators/wells/MultisegmentWellEquations.hpp
+++ b/opm/simulators/wells/MultisegmentWellEquations.hpp
@@ -77,6 +77,9 @@ public:
     //! \brief Apply linear operator to vector.
     void apply(const BVector& x, BVector& Ax) const;
 
+    //! \brief Apply linear operator to vector.
+    void apply(BVector& r) 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 7c62470f8..bb263b96c 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -217,14 +217,11 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     apply(BVector& r) const
     {
-        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
+            return;
+        }
 
-        // invDrw_ = duneD^-1 * resWell_
-        const BVectorWell invDrw = mswellhelpers::applyUMFPack(this->linSys_.duneD_,
-                                                               this->linSys_.duneDSolver_,
-                                                               this->linSys_.resWell_);
-        // r = r - duneC_^T * invDrw
-        this->linSys_.duneC_.mmtv(invDrw, r);
+        this->linSys_.apply(r);
     }