diff --git a/lib/solvers/OptimalVcycleMultigridSolver.m b/lib/solvers/OptimalVcycleMultigridSolver.m
index 9c3cc0baa13945cc83fef787961aa38e967c5358..496c2a73082af19af51919bf7ebb1c0ae04f8572 100644
--- a/lib/solvers/OptimalVcycleMultigridSolver.m
+++ b/lib/solvers/OptimalVcycleMultigridSolver.m
@@ -1,6 +1,7 @@
 % OptimalVcycleMultigridSolver (subclass of IterativeSolver) Solves linear
-%   equations iteratively, using a geometric multigrid solver with V-cycle and
-%   specified MultiLevelSmoother.
+%   equations iteratively, using a geometric multigrid solver with specified
+%   MultiLevelSmoother.
+%   One iteration = Vcycle(0,1) (no pre-/one post-smoothing step).
 %
 %   To obtain the solver described in [Innerberger, Miraçi, Praetorius,
 %   Streitberger; Optimal computational costs of AFEM with optimal local
@@ -28,7 +29,7 @@ classdef OptimalVcycleMultigridSolver < IterativeSolver
     
     properties (Access=protected)
         smoother
-        matrix
+        projectedMatrix
         lhsNorm
     end
     
@@ -46,8 +47,8 @@ classdef OptimalVcycleMultigridSolver < IterativeSolver
             % override superclass method, because ALL matrices are needed
             setupSystemMatrix@IterativeSolver(obj, A, varargin{:});
 
-            obj.smoother.setup(A, varargin{:});
-            obj.matrix{obj.nLevels} = A;
+            projectedA = obj.smoother.setup(A, varargin{:});
+            obj.projectedMatrix{obj.nLevels} = projectedA;
         end
 
         function setupRhs(obj, b, varargin)
@@ -101,20 +102,26 @@ classdef OptimalVcycleMultigridSolver < IterativeSolver
             end
 
             % exact solve on coarsest level to compute accumulative lifting of x
-            Cx = obj.matrix{1} \ res{1};
-            eta2 = dot(Cx, obj.matrix{1} * Cx, Dim.Vector);
+            Cx = obj.projectedMatrix{1} \ res{1};
+            eta2 = dot(Cx, obj.projectedMatrix{1} * Cx, Dim.Vector);
 
             % ascending cascade: (local) smoothing and compute error estimator
-            for k = 2:L
+            for k = 2:(L-1)
                 Cx = obj.smoother.prolongate(Cx, k);
-                updatedResidual = res{k} - obj.matrix{k} * Cx;
+                updatedResidual = res{k} - obj.projectedMatrix{k} * Cx;
                 rho = obj.smoother.smooth(updatedResidual, k);
-                [etaUpdate2, xUpdate] = computeOptimalUpdate(obj.matrix{k}, updatedResidual, rho);
+                [etaUpdate2, xUpdate] = computeOptimalUpdate(obj.projectedMatrix{k}, updatedResidual, rho);
                 Cx = Cx + xUpdate;
                 eta2 = eta2 + etaUpdate2;
             end
 
-            algEstimator = sqrt(eta2);
+            % final smoothing and residual update (with non-projected matrix)
+            Cx = obj.smoother.prolongate(Cx, L);
+            updatedResidual = res{L} - obj.A * Cx;
+            rho = obj.smoother.smooth(updatedResidual, L);
+            [etaUpdate2, xUpdate] = computeOptimalUpdate(obj.A, updatedResidual, rho);
+            Cx = Cx + xUpdate;
+            algEstimator = sqrt(eta2 + etaUpdate2);
         end
     end
 end
diff --git a/lib/solvers/smoothers/MultilevelSmoother.m b/lib/solvers/smoothers/MultilevelSmoother.m
index 45ff623060c116210a886c73c0587de8027eb30a..0616a0300625badd384be2df0da7087d56af112e 100644
--- a/lib/solvers/smoothers/MultilevelSmoother.m
+++ b/lib/solvers/smoothers/MultilevelSmoother.m
@@ -5,6 +5,9 @@
 %
 % smoother.setup(A) hook to set up the smoother for the matrix A.
 %
+% nonInvertedSmootherMatrix = smoother.setup(A) set up the smoother and also
+%   return the non-inverted matrix, based on which the smoother is computed.
+%
 % Cx = smoother.smooth(x, k) applies the smoother to the vector x on level k.
 %
 % Px = smoother.prolongate(x, k) prolongate vector x from level k to level k+1.
@@ -20,14 +23,32 @@ classdef MultilevelSmoother < handle
         nLevels
     end
 
+    properties (Access=protected)
+        fes
+        blf
+    end
+
     %% methods
     methods (Access=public)
-        function obj = MultilevelSmoother()
+        function obj = MultilevelSmoother(fes, blf)
+            arguments
+                fes FeSpace
+                blf BilinearForm
+            end
+
+            assert(isempty(blf.b), ...
+                'Multilevel smoothers only tested for symmetric problems.')
+
             obj.nLevels = 0;
+            obj.fes = fes;
+            obj.blf = blf;
         end
 
-        function setup(obj, ~)
+        function nonInvertedSmootherMatrix = setup(obj, ~)
             obj.nLevels = obj.nLevels + 1;
+
+            % should be implemented by subclasses
+            nonInvertedSmootherMatrix = [];
         end
     end
 
diff --git a/lib/solvers/smoothers/P1JacobiSmoother.m b/lib/solvers/smoothers/P1JacobiSmoother.m
index 59d4311240877cbe1245f1d580ae5e4aa8e43a86..91a82f955b4f64675e7967b963f90cbba2a2189c 100644
--- a/lib/solvers/smoothers/P1JacobiSmoother.m
+++ b/lib/solvers/smoothers/P1JacobiSmoother.m
@@ -7,10 +7,8 @@
 
 classdef P1JacobiSmoother < MultilevelSmoother
     properties (Access=protected)
-        fes
-        blf
-        inverseDiagonal
         P
+        inverseDiagonal
         intergridMatrix
         freeVertices
         freeVerticesOld
@@ -26,30 +24,28 @@ classdef P1JacobiSmoother < MultilevelSmoother
     methods (Access=public)
         function obj = P1JacobiSmoother(fes, blf, P)
             arguments
-                fes FeSpace
-                blf BilinearForm
+                fes
+                blf
                 P Prolongation
             end
-            obj = obj@MultilevelSmoother();
+            obj = obj@MultilevelSmoother(fes, blf);
             
             assert(fes.finiteElement.order == 1, ...
                 'P1JacobiSmoother only works for lowest order finite elements.')
             
-            obj.fes = fes;
             obj.P = P;
-            obj.blf = blf;
-            
             mesh = fes.mesh;
             obj.listenerHandle = mesh.listener('IsAboutToRefine', @obj.getChangedPatches);
         end
 
-        function setup(obj, A)
+        function nonInvertedSmootherMatrix = setup(obj, A)
             setup@MultilevelSmoother(obj, A);
             
             L = obj.nLevels;
             obj.freeVerticesOld = obj.freeVertices;
             obj.freeVertices = getFreeDofs(obj.fes);
             
+            nonInvertedSmootherMatrix = A;
             obj.inverseDiagonal{L} = full(diag(A)).^(-1);
             
             if L >= 2