diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp
index 1e0d3621e..8cfab7157 100644
--- a/opm/autodiff/BlackoilModelEbos.hpp
+++ b/opm/autodiff/BlackoilModelEbos.hpp
@@ -98,6 +98,7 @@ namespace Opm {
     /// upwind weighting of mobilities.
     class BlackoilModelEbos
     {
+        typedef BlackoilModelEbos ThisType;
     public:
         // ---------  Types and enums  ---------
         typedef AutoDiffBlock<double> ADB;
@@ -362,6 +363,12 @@ namespace Opm {
             return linsolver_.iterations();
         }
 
+        template <class X, class Y>
+        void applyWellModel(const X& x, Y& y )
+        {
+           wellModel().apply(x, y);
+        }
+
 
         /// Solve the Jacobian system Jx = r where J is the Jacobian and
         /// r is the residual.
@@ -374,17 +381,19 @@ namespace Opm {
             typedef Dune::BCRSMatrix <MatrixBlockType>      Mat;
             typedef Dune::BlockVector<VectorBlockType>      BVector;
 
-            typedef Dune::MatrixAdapter<Mat,BVector,BVector> Operator;
             auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
             auto& ebosResid = ebosSimulator_.model().linearizer().residual();
 
-            Operator opA(ebosJac);
+            typedef WellModelMatrixAdapter<Mat,BVector,BVector, ThisType> Operator;
+            Operator opA(ebosJac, const_cast< ThisType&  > (*this));
             const double relax = 0.9;
             typedef Dune::SeqILU0<Mat, BVector, BVector> SeqPreconditioner;
             SeqPreconditioner precond(opA.getmat(), relax);
             Dune::SeqScalarProduct<BVector> sp;
 
-            wellModel().apply(ebosJac, ebosResid);
+            // apply well residual to the residual.
+            wellModel().apply(ebosResid);
+
             Dune::BiCGSTABSolver<BVector> linsolve(opA, sp, precond,
                                                    0.01,
                                                    100,
@@ -393,10 +402,49 @@ namespace Opm {
             x = 0.0;
             linsolve.apply(x, ebosResid, result);
 
+            // recover wells.
             xw = 0.0;
             wellModel().recoverVariable(x, xw);
         }
 
+        //=====================================================================
+        // Implementation for ISTL-matrix based operator
+        //=====================================================================
+
+        /*!
+           \brief Adapter to turn a matrix into a linear operator.
+
+           Adapts a matrix to the assembled linear operator interface
+         */
+        template<class M, class X, class Y, class WellModel>
+        class WellModelMatrixAdapter : public Dune::MatrixAdapter<M,X,Y>
+        {
+          typedef Dune::MatrixAdapter<M,X,Y> BaseType;
+        public:
+          //! export types
+          typedef M matrix_type;
+          typedef X domain_type;
+          typedef Y range_type;
+          typedef typename X::field_type field_type;
+
+          //! constructor: just store a reference to a matrix
+          explicit WellModelMatrixAdapter (const M& A, WellModel& wellMod ) : BaseType( A ), wellMod_( wellMod ) {}
+
+          //! apply operator to x:  \f$ y = A(x) \f$
+          virtual void apply (const X& x, Y& y) const
+          {
+            BaseType::apply( x, y );
+
+            wellMod_.applyWellModel(x, y );
+          }
+
+        private:
+          WellModel& wellMod_;
+        };
+
+        /** @} end documentation */
+
+
         /// Apply an update to the primary variables, chopped if appropriate.
         /// \param[in]      dx                updates to apply to primary variables
         /// \param[in, out] reservoir_state   reservoir state variables
diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp
index bf74cde17..6e532ca74 100644
--- a/opm/autodiff/StandardWellsDense.hpp
+++ b/opm/autodiff/StandardWellsDense.hpp
@@ -238,87 +238,26 @@ namespace Opm {
                 }
             }
 
-            Mat matAdd( const Mat& A, const Mat& B ) const
-            {
-                return matBinaryOp( A, B, true );
-            }
-
-            Mat matSubstract( const Mat& A, const Mat& B ) const
-            {
-                return matBinaryOp( A, B, false );
-            }
-
-            Mat matBinaryOp( const Mat& A, const Mat& B, const bool add ) const
-            {
-                const int avg_cols_per_row = 20;
-                const double overflow_fraction = 0.4;
-                Mat res( A.N(), A.M(), avg_cols_per_row,overflow_fraction, Mat::implicit );
-                assert( A.N() == B.N() && A.M() == B.M() );
-
-                // res = A
-                for( auto rowit = A.begin(), rowEnd = A.end(); rowit != rowEnd; ++rowit )
-                {
-                    const int rowIdx = rowit.index();
-                    const auto colEnd = rowit->end();
-                    for( auto colit = rowit->begin(); colit != colEnd; ++colit )
-                    {
-                        const int colIdx =  colit.index();
-                        res.entry( rowIdx, colIdx ) = (*colit);
-                    }
-                }
-
-                // res += B
-                for( auto rowit = B.begin(), rowEnd = B.end(); rowit != rowEnd; ++rowit )
-                {
-                    const int rowIdx = rowit.index();
-                    const auto colEnd = rowit->end();
-                    for( auto colit = rowit->begin(); colit != colEnd; ++colit )
-                    {
-                        const int colIdx =  colit.index();
-                        if( add )
-                        {
-                            res.entry( rowIdx, colIdx ) += (*colit);
-                        }
-                        else
-                        {
-                            res.entry( rowIdx, colIdx ) -= (*colit);
-                        }
-                    }
-                }
-
-                res.compress();
-                return res;
-            }
-
             void addRhs(BVector& x, Mat& jac) const {
                 assert(x.size() == rhs.size());
                 x += rhs_;
-                // jac = A + duneA
-                //jac = matAdd( jac, duneA_ );
                 jac += duneA_;
             }
 
-            void apply( Mat& A,
-                        BVector& res) const {
+            // substract Binv(D)rw from r;
+            void apply( BVector& r) const {
+                BVector invDrw(invDuneD_.N());
+                invDuneD_.mv(resWell_,invDrw);
+                duneB_.mmv(invDrw, r);
+            }
 
-                Mat BmultinvD;
-                Mat duneA;
-
-                Dune::matMultMat(BmultinvD, duneB_ , invDuneD_);
-                Dune::matMultMat(duneA, BmultinvD, duneC_);
-                //std::cout << "before" << std::endl;
-                //std::cout << "A" << std::endl;
-
-                //print(A);
-                //std::cout << "duneA" << std::endl;
-
-                //print(duneA);
-                // A = E - duneA
-                A = matSubstract( A, duneA );
-                //A -= duneA;
-                //std::cout << "after" << std::endl;
-                //print(A);
-                BmultinvD.mmv(resWell_, res);
+            // subtract B*inv(D)*C * x from A*x
+            void apply(const BVector& x, BVector& Ax) {
+                BVector Cx(duneC_.N());
+                duneC_.mv(x, Cx);
+                BVector invDCx(invDuneD_.N());
+                invDuneD_.mv(Cx, invDCx);
+                duneB_.mmv(invDCx,Ax);
             }
 
             void recoverVariable(const BVector& x, BVector& xw) const {