diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
index fa5236bd99..11ece2171d 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
@@ -57,22 +57,7 @@ Foam::tmp<Foam::Field<Type> > Foam::lduMatrix::H(const Field<Type>& psi) const
 
         for (register label face=0; face<nFaces; face++)
         {
-            #ifdef ICC_IA64_PREFETCH
-            __builtin_prefetch (&uPtr[face+32],0,0);
-            __builtin_prefetch (&lPtr[face+32],0,0);
-            __builtin_prefetch (&lowerPtr[face+32],0,1);
-            __builtin_prefetch (&psiPtr[lPtr[face+32]],0,1);
-            __builtin_prefetch (&HpsiPtr[uPtr[face+32]],0,1);
-            #endif
-
             HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
-        
-            #ifdef ICC_IA64_PREFETCH
-            __builtin_prefetch (&upperPtr[face+32],0,1);
-            __builtin_prefetch (&psiPtr[uPtr[face+32]],0,1);
-            __builtin_prefetch (&HpsiPtr[lPtr[face+32]],0,1);
-                #endif
-
             HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
         }
     }
@@ -94,22 +79,33 @@ template<class Type>
 Foam::tmp<Foam::Field<Type> >
 Foam::lduMatrix::faceH(const Field<Type>& psi) const
 {
-    const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
-    const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();
-
-    // Take refereces to addressing
-    const unallocLabelList& l = lduAddr().lowerAddr();
-    const unallocLabelList& u = lduAddr().upperAddr();
-
-    tmp<Field<Type> > tfaceHpsi(new Field<Type> (Lower.size()));
-    Field<Type> & faceHpsi = tfaceHpsi();
-
-    for (register label face=0; face<l.size(); face++)
+    if (lowerPtr_ || upperPtr_)
     {
-        faceHpsi[face] = Upper[face]*psi[u[face]] - Lower[face]*psi[l[face]];
-    }
+        const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
+        const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();
 
-    return tfaceHpsi;
+        const unallocLabelList& l = lduAddr().lowerAddr();
+        const unallocLabelList& u = lduAddr().upperAddr();
+
+        tmp<Field<Type> > tfaceHpsi(new Field<Type> (Lower.size()));
+        Field<Type> & faceHpsi = tfaceHpsi();
+
+        for (register label face=0; face<l.size(); face++)
+        {
+            faceHpsi[face] = Upper[face]*psi[u[face]] - Lower[face]*psi[l[face]];
+        }
+
+        return tfaceHpsi;
+    }
+    else
+    {
+        FatalErrorIn("lduMatrix::faceH(const Field<Type>& psi) const")
+            << "Cannot calculate faceH"
+               " the matrix does not have any off-diagonal coefficients."
+            << exit(FatalError);
+
+        return tmp<Field<Type> >(NULL);
+    }
 }
 
 
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 0fcf2e626d..78fede73f8 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -1303,7 +1303,14 @@ Foam::tmp<Foam::fvMatrix<Type> > Foam::correction
     const fvMatrix<Type>& A
 )
 {
-    return A - (A & A.psi());
+    tmp<Foam::fvMatrix<Type> > tAcorr = A - (A & A.psi());
+
+    if ((A.hasUpper() || A.hasLower()) && A.mesh().fluxRequired(A.psi().name()))
+    {
+        tAcorr().faceFluxCorrectionPtr() = (-A.flux()).ptr();
+    }
+
+    return tAcorr;
 }
 
 
@@ -1313,7 +1320,17 @@ Foam::tmp<Foam::fvMatrix<Type> > Foam::correction
     const tmp<fvMatrix<Type> >& tA
 )
 {
-    return tA - (tA() & tA().psi());
+    tmp<Foam::fvMatrix<Type> > tAcorr = tA - (tA() & tA().psi());
+
+    // Note the matrix coefficients are still that of matrix A
+    const fvMatrix<Type>& A = tAcorr();
+
+    if ((A.hasUpper() || A.hasLower()) && A.mesh().fluxRequired(A.psi().name()))
+    {
+        tAcorr().faceFluxCorrectionPtr() = (-A.flux()).ptr();
+    }
+
+    return tAcorr;
 }