diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAMGSolver/BlockAMGSolver.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAMGSolver/BlockAMGSolver.C
index cae83658b..926988fbd 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAMGSolver/BlockAMGSolver.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAMGSolver/BlockAMGSolver.C
@@ -72,10 +72,10 @@ Foam::BlockAMGSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     // Calculate initial residual
-    solverPerf.initialResidual() = gSum(cmptMag(amg_.residual(x, b)))/norm;
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(amg_.residual(x, b))),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Stop solver on divergence
@@ -88,8 +88,8 @@ Foam::BlockAMGSolver<Type>::solve
         {
             amg_.cycle(x, b);
 
-            solverPerf.finalResidual() =
-                gSum(cmptMag(amg_.residual(x, b)))/norm;
+              solverPerf.finalResidual() =
+                  cmptDivide(gSum(cmptMag(amg_.residual(x, b))), norm);
 
             solverPerf.nIterations()++;
 
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C
index d7776e867..1a9516976 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C
@@ -76,7 +76,7 @@ Foam::BlockBiCGStabSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     // Multiplication helper
     typename BlockCoeff<Type>::multiply mult;
@@ -87,7 +87,8 @@ Foam::BlockBiCGStabSolver<Type>::solve
     matrix.Amul(p, x);
     Field<Type> r(b - p);
 
-    solverPerf.initialResidual() = gSum(cmptMag(r))/norm;
+    // NOTE: Normalisation of residual per component! TU, Feb 2019
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(r)),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Check convergence, solve if not converged
@@ -160,7 +161,7 @@ Foam::BlockBiCGStabSolver<Type>::solve
                 r[i] = s[i] - omega*t[i];
             }
 
-            solverPerf.finalResidual() = gSum(cmptMag(r))/norm;
+            solverPerf.finalResidual() = cmptDivide(gSum(cmptMag(r)), norm);
             solverPerf.nIterations()++;
         } while (!this->stop(solverPerf));
     }
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C
index 7f6401a08..375fc9a53 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C
@@ -75,7 +75,7 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     Field<Type> wA(x.size());
 
@@ -83,7 +83,8 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
     matrix.Amul(wA, x);
     Field<Type> rA(b - wA);
 
-    solverPerf.initialResidual() = gSum(cmptMag(rA))/norm;
+    // NOTE: Normalisation of residual per component! TU, Feb 2019
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(rA)),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Check convergence, solve if not converged
@@ -120,7 +121,7 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
             wApA = gSumProd(wA, pA);
 
             // Check for singularity
-            if (solverPerf.checkSingularity(mag(wApA)/norm))
+            if (solverPerf.checkSingularity(mag(wApA)/mag(norm)))
             {
                 break;
             }
@@ -138,7 +139,7 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
                 rA[i] -= alpha*wA[i];
             }
 
-            solverPerf.finalResidual() = gSum(cmptMag(rA))/norm;
+            solverPerf.finalResidual() = cmptDivide(gSum(cmptMag(rA)), norm);
             solverPerf.nIterations()++;
         } while (!this->stop(solverPerf));
     }
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C
index 93486ec26..e1030ed03 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C
@@ -111,7 +111,7 @@ Foam::BlockGMRESSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     Field<Type> wA(x.size());
 
@@ -119,7 +119,7 @@ Foam::BlockGMRESSolver<Type>::solve
     matrix.Amul(wA, x);
     Field<Type> rA(b - wA);
 
-    solverPerf.initialResidual() = gSum(cmptMag(rA))/norm;
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(rA)),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Check convergence, solve if not converged
@@ -234,7 +234,7 @@ Foam::BlockGMRESSolver<Type>::solve
                 rA[raI] = b[raI] - wA[raI];
             }
 
-            solverPerf.finalResidual() = gSum(cmptMag(rA))/norm;
+            solverPerf.finalResidual() = cmptDivide(gSum(cmptMag(rA)), norm);
             solverPerf.nIterations()++;
         } while (!this->stop(solverPerf));
     }
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C
index c86005c4b..22ae14e8f 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C
@@ -72,7 +72,7 @@ Foam::BlockGaussSeidelSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     Field<Type> wA(x.size());
 
@@ -80,7 +80,7 @@ Foam::BlockGaussSeidelSolver<Type>::solve
     matrix.Amul(wA, x);
     wA -= b;
 
-    solverPerf.initialResidual() = gSum(cmptMag(wA))/norm;
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(wA)),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Check convergence, solve if not converged
@@ -103,7 +103,7 @@ Foam::BlockGaussSeidelSolver<Type>::solve
             matrix.Amul(wA, x);
             wA -= b;
 
-            solverPerf.finalResidual() = gSum(cmptMag(wA))/norm;
+            solverPerf.finalResidual() = cmptDivide(gSum(cmptMag(wA)), norm);
         } while (!this->stop(solverPerf));
     }
 
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockILU/BlockILUSolver.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockILU/BlockILUSolver.C
index c8af28ec3..2375f6524 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockILU/BlockILUSolver.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockILU/BlockILUSolver.C
@@ -72,7 +72,7 @@ Foam::BlockILUSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     Field<Type> residual(x.size());
 
@@ -84,7 +84,7 @@ Foam::BlockILUSolver<Type>::solve
         residual[i] = b[i] - residual[i];
     }
 
-    solverPerf.initialResidual() = gSum(cmptMag(residual))/norm;
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(residual)),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Check convergence, solve if not converged
@@ -114,7 +114,7 @@ Foam::BlockILUSolver<Type>::solve
                 residual[i] = b[i] - residual[i];
             }
 
-            solverPerf.finalResidual() = gSum(cmptMag(residual))/norm;
+            solverPerf.finalResidual() = cmptDivide(gSum(cmptMag(residual)), norm);
         } while (!this->stop(solverPerf));
     }
 
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C
index e322fa83c..7d86e07eb 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C
@@ -50,7 +50,7 @@ Foam::BlockIterativeSolver<Type>::BlockIterativeSolver
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
 template<class Type>
-Foam::scalar Foam::BlockIterativeSolver<Type>::normFactor
+Type Foam::BlockIterativeSolver<Type>::normFactor
 (
     Field<Type>& x,
     const Field<Type>& b
@@ -77,7 +77,13 @@ Foam::scalar Foam::BlockIterativeSolver<Type>::normFactor
         Field<Type>(nRows, xRef)
     );
 
-    scalar normFactor = gSum(mag(wA - pA) + mag(b - pA)) + this->small_;
+    // The components of the normalisation factor cannot be equal to 0 since it
+    // appears in the denominator: replace zero with small_
+    Type small = pTraits<Type>::one*this->small_;
+
+    // Norm factor for the block matrix is not a scalar since each component
+    // should be normalised with a value corresponding to a particular variable
+    Type normFactor = gSum(cmptMag(wA - pA) + cmptMag(b - pA)) + small;
 
     if (blockLduMatrix::debug >= 2)
     {
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.H
index e226ad9d4..f92497a41 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.H
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.H
@@ -84,7 +84,7 @@ protected:
     // Protected Member Functions
 
         //- Return normalisation factor
-        scalar normFactor
+        Type normFactor
         (
             Field<Type>& x,
             const Field<Type>& b
@@ -107,10 +107,9 @@ public:
         );
 
 
-    // Destructor
-
-        virtual ~BlockIterativeSolver()
-        {}
+    //- Destructor
+    virtual ~BlockIterativeSolver()
+    {}
 
 
     // Member Functions
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C
index 0874b1307..c2875c1bf 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C
@@ -39,14 +39,9 @@ bool Foam::BlockSolverPerformance<Type>::checkConvergence
     {
         Info<< solverName_
             << ":  Iteration " << nIterations_
-            << " residual = " << finalResidual_;
-
-        if (pTraits<Type>::rank > 0)
-        {
-            Info<< " mag = " << mag(finalResidual_);
-        }
-
-        Info<< " tol = "
+            << " residual = " << finalResidual_
+            << " mag = " << mag(finalResidual_)
+            << " tol = "
             << Foam::max(Tolerance, RelTolerance*mag(initialResidual_))
             << endl;
     }
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C
index 490868c62..7382be091 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C
@@ -41,6 +41,7 @@ defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceVector, 1);
 defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceSphericalTensor, 1);
 defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceSymmTensor, 1);
 defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceTensor, 1);
+
 };