diff --git a/src/foam/matrices/blockLduMatrix/BlockAmg/BlockAmgCycle.C b/src/foam/matrices/blockLduMatrix/BlockAmg/BlockAmgCycle.C
index 1607b4670..d325fa3a3 100644
--- a/src/foam/matrices/blockLduMatrix/BlockAmg/BlockAmgCycle.C
+++ b/src/foam/matrices/blockLduMatrix/BlockAmg/BlockAmgCycle.C
@@ -133,10 +133,7 @@ void Foam::BlockAmgCycle<Type>::fixedCycle
         Field<Type>& bCoarse = coarseLevelPtr_->levelPtr_->b();
 
         // Zero out coarse x
-        forAll (xCoarse,i)
-        {
-            xCoarse[i] *= 0.0;
-        }
+        xCoarse = pTraits<Type>::zero;
 
         // Restrict residual: optimisation on number of pre-sweeps
         levelPtr_->restrictResidual
diff --git a/src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.C b/src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.C
index 1d3b08de5..728c87c0e 100644
--- a/src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.C
+++ b/src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.C
@@ -25,7 +25,7 @@ Class
     BlockMatrixAgglomeration
 
 Description
-    Agglomerative block matrix AMG corsening, adjusted for BlockLduMatrix
+    Agglomerative block matrix AMG corsening
 
 Author
     Klas Jareteg, 2012-12-13
@@ -33,6 +33,7 @@ Author
 \*---------------------------------------------------------------------------*/
 
 #include "BlockMatrixAgglomeration.H"
+#include "boolList.H"
 #include "coeffFields.H"
 #include "addToRunTimeSelectionTable.H"
 #include "BlockGAMGInterfaceField.H"
@@ -40,6 +41,12 @@ Author
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
+template<class Type>
+const Foam::scalar Foam::BlockMatrixAgglomeration<Type>::diagFactor_
+(
+    debug::tolerances("aamgDiagFactor", 1e-8)
+);
+
 template<class Type>
 const Foam::scalar Foam::BlockMatrixAgglomeration<Type>::weightFactor_ = 0.65;
 
@@ -47,7 +54,7 @@ const Foam::scalar Foam::BlockMatrixAgglomeration<Type>::weightFactor_ = 0.65;
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class Type>
-void Foam::BlockMatrixAgglomeration<Type>::calcChild()
+void Foam::BlockMatrixAgglomeration<Type>::calcAgglomeration()
 {
     // Algorithm:
     // 1) Create temporary equation addressing using a double-pass algorithm.
@@ -96,9 +103,9 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
 
     // Create column and coefficient index array
     {
-        // Use child array to count number of entries per row.
+        // Use agglomIndex to count number of entries per row.
         // Reset the list for counting
-        labelList& nPerRow = child_;
+        labelList& nPerRow = agglomIndex_;
         nPerRow = 0;
 
         forAll (upperAddr, coeffI)
@@ -123,8 +130,8 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
             nPerRow[lowerAddr[coeffI]]++;
         }
 
-        // Reset child array
-        child_ = -1;
+        // Reset agglomeration index array
+        agglomIndex_ = -1;
     }
 
 
@@ -132,43 +139,122 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
 
     // Get matrix coefficients
     const CoeffField<Type>& diag = matrix_.diag();
-    const CoeffField<Type>& upper = matrix_.upper();
+
+    // Coefficient magnitudes are pre-calculated
+    scalarField magDiag(diag.size());
+    scalarField magOffDiag(upperAddr.size());
+
+    normPtr_->coeffMag(diag, magDiag);
+
+    if (matrix_.asymmetric())
+    {
+        scalarField magUpper(upperAddr.size());
+        scalarField magLower(upperAddr.size());
+
+        normPtr_->coeffMag(matrix_.upper(), magUpper);
+        normPtr_->coeffMag(matrix_.lower(), magLower);
+
+        magOffDiag = Foam::max(magUpper, magLower);
+    }
+    else if (matrix_.symmetric())
+    {
+        normPtr_->coeffMag(matrix_.upper(), magOffDiag);
+    }
+    else
+    {
+        // Diag only matrix.  Reset and return
+        agglomIndex_ = 0;
+        nCoarseEqns_ = 1;
+
+        return;
+    }
 
     labelList sizeOfGroups(nRows, 0);
 
     nCoarseEqns_ = 0;
 
-    label indexU, indexG, colI, curEqn, nextEqn, groupPassI;
+    // Gather disconnected and weakly connected equations into cluster zero
+    // Weak connection is assumed to be the one where the off-diagonal
+    // coefficient is smaller than diagFactor_*diag
+    {
+        // Algorithm
+        // Mark all cells to belong to zero cluster
+        // Go through all upper and lower coefficients and for the ones
+        // larger than threshold mark the equations out of cluster zero
 
-    scalar magRowDiag, magColDiag, weight, weightU, weightG;
+        scalarField magScaledDiag = diagFactor_*magDiag;
 
-    // Magnitudes pre-calculated
-    Field<scalar> magDiag(diag.size());
-    Field<scalar> magUpper(upper.size());
+        boolList zeroCluster(diag.size(), true);
 
-    normPtr_->coeffMag(diag,magDiag);
-    normPtr_->coeffMag(upper,magUpper);
+        forAll (magOffDiag, coeffI)
+        {
+            if (magOffDiag[coeffI] > magScaledDiag[upperAddr[coeffI]])
+            {
+                zeroCluster[upperAddr[coeffI]] = false;
+            }
+
+            if (magOffDiag[coeffI] > magScaledDiag[lowerAddr[coeffI]])
+            {
+                zeroCluster[lowerAddr[coeffI]] = false;
+            }
+        }
+
+        // Collect solo equations
+        label nSolo = 0;
+
+        forAll (zeroCluster, eqnI)
+        {
+            if (zeroCluster[eqnI])
+            {
+                // Found solo equation
+                nSolo++;
+
+                agglomIndex_[eqnI] = nCoarseEqns_;
+            }
+        }
+
+        reduce(nSolo, sumOp<label>());
+
+        if (nSolo > 0)
+        {
+            // Found solo equations
+            nCoarseEqns_++;
+
+            if (BlockLduMatrix<Type>::debug >= 2)
+            {
+                Info<< "Found " << nSolo << " weakly connected equations."
+                    << endl;
+            }
+        }
+    }
+
+    // Go through the off-diagonal and create clusters, marking the child array
+    label indexUngrouped, indexGrouped;
+    label colI, curEqn, nextEqn, groupPassI;
+
+    scalar magRowDiag, magColDiag;
+    scalar weight, weightUngrouped, weightGrouped;
 
     for (label eqnI = 0; eqnI < nRows; eqnI++)
     {
-        if (child_[eqnI] == -1)
+        if (agglomIndex_[eqnI] == -1)
         {
             curEqn = eqnI;
 
-            indexU = -1;
-            indexG = -1;
+            indexUngrouped = -1;
+            indexGrouped = -1;
 
-            child_[curEqn] = nCoarseEqns_;
+            agglomIndex_[curEqn] = nCoarseEqns_;
 
             magRowDiag = magDiag[curEqn];
 
             for (groupPassI = 1; groupPassI < groupSize_; groupPassI++)
             {
-                weightU = 0;
-                weightG = 0;
+                weightUngrouped = 0;
+                weightGrouped = 0;
 
-                indexU = -1;
-                indexG = -1;
+                indexUngrouped = -1;
+                indexGrouped = -1;
 
                 for
                 (
@@ -181,40 +267,42 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
 
                     magColDiag = magDiag[colI];
 
-                    weight = (magUpper[cIndex[rowCoeffI]]
-                        / max(magRowDiag, magColDiag));
+                    weight = magOffDiag[cIndex[rowCoeffI]]/max(magRowDiag, magColDiag);
 
-                    if (child_[colI] == -1)
+                    if (agglomIndex_[colI] == -1)
                     {
-                        if (indexU == -1 || weight > weightU)
+                        if (indexUngrouped == -1 || weight > weightUngrouped)
                         {
-                            indexU = rowCoeffI;
-                            weightU = weight;
+                            indexUngrouped = rowCoeffI;
+                            weightUngrouped = weight;
                         }
                     }
-                    else if (child_[curEqn] != child_[colI])
+                    else if (agglomIndex_[curEqn] != agglomIndex_[colI])
                     {
-                        if (indexG == -1 || weight > weightG)
+                        if (indexGrouped == -1 || weight > weightGrouped)
                         {
-                            indexG = rowCoeffI;
-                            weightG = weight;
+                            indexGrouped = rowCoeffI;
+                            weightGrouped = weight;
                         }
                     }
                 }
 
                 if
                 (
-                    indexU != -1
-                 && (indexG == -1 || weightU >= weightFactor_*weightG)
+                    indexUngrouped != -1
+                 && (
+                        indexGrouped == -1
+                     || weightUngrouped >= weightFactor_*weightGrouped
+                    )
                 )
                 {
                     // Found new element of group.  Add it and use as
                     // start of next search
 
-                    nextEqn = cols[indexU];
+                    nextEqn = cols[indexUngrouped];
 
-                    child_[nextEqn] = child_[curEqn];
-                    sizeOfGroups[child_[curEqn]];
+                    agglomIndex_[nextEqn] = agglomIndex_[curEqn];
+                    sizeOfGroups[agglomIndex_[curEqn]];
 
                     curEqn = nextEqn;
                 }
@@ -228,17 +316,20 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
             if
             (
                 groupPassI > 1
-             || indexG == -1
-             || sizeOfGroups[child_[cols[indexG]]] > (groupSize_ + 2)
+             || indexGrouped == -1
+             || (
+                    sizeOfGroups[agglomIndex_[cols[indexGrouped]]]
+                  > (groupSize_ + 2)
+                )
             )
             {
-                sizeOfGroups[child_[eqnI]]++;
+                sizeOfGroups[agglomIndex_[eqnI]]++;
                 nCoarseEqns_++;
             }
             else
             {
-                child_[eqnI] = child_[cols[indexG]];
-                sizeOfGroups[child_[cols[indexG]]]++;
+                agglomIndex_[eqnI] = agglomIndex_[cols[indexGrouped]];
+                sizeOfGroups[agglomIndex_[cols[indexGrouped]]]++;
             }
         }
     }
@@ -258,6 +349,19 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
 
     reduce(coarsen_, andOp<bool>());
 
+    if (BlockLduMatrix<Type>::debug >= 2)
+    {
+        Pout << "Coarse level size: " << nCoarseEqns_;
+
+        if (coarsen_)
+        {
+            Pout << ".  Accepted" << endl;
+        }
+        else
+        {
+            Pout << ".  Rejected" << endl;
+        }
+    }
 }
 
 
@@ -274,19 +378,13 @@ Foam::BlockMatrixAgglomeration<Type>::BlockMatrixAgglomeration
 :
     BlockMatrixCoarsening<Type>(matrix, dict, groupSize, minCoarseEqns),
     matrix_(matrix),
-    normPtr_
-    (
-        BlockCoeffNorm<Type>::New
-        (
-            dict
-        )
-    ),
-    child_(matrix_.lduAddr().size()),
+    normPtr_(BlockCoeffNorm<Type>::New(dict)),
+    agglomIndex_(matrix_.lduAddr().size()),
     groupSize_(groupSize),
     nCoarseEqns_(0),
     coarsen_(false)
 {
-    calcChild();
+    calcAgglomeration();
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -314,7 +412,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
 
     // Construct the coarse matrix and ldu addressing for the next level
     // Algorithm:
-    // 1) Loop through all fine coeffs. If the child labels on two sides are
+    // 1) Loop through all fine coeffs.  If the agglom labels on two sides are
     //    different, this creates a coarse coeff. Define owner and neighbour
     //    for this coeff based on cluster IDs.
     // 2) Check if the coeff has been seen before. If yes, add the coefficient
@@ -332,14 +430,14 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
     const label nFineCoeffs = upperAddr.size();
 
 #   ifdef FULLDEBUG
-    if (child_.size() != matrix_.lduAddr().size())
+    if (agglomIndex_.size() != matrix_.lduAddr().size())
     {
         FatalErrorIn
         (
             "autoPtr<BlockLduMatrix<Type> > BlockMatrixAgglomeration<Type>::"
             "restrictMatrix() const"
-        )   << "Child array does not correspond to fine level. " << endl
-            << " Child size: " << child_.size()
+        )   << "agglomIndex array does not correspond to fine level. " << endl
+            << " Size: " << agglomIndex_.size()
             << " number of equations: " << matrix_.lduAddr().size()
             << abort(FatalError);
     }
@@ -369,8 +467,8 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
     // Loop through all fine coeffs
     forAll (upperAddr, fineCoeffi)
     {
-        label rmUpperAddr = child_[upperAddr[fineCoeffi]];
-        label rmLowerAddr = child_[lowerAddr[fineCoeffi]];
+        label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]];
+        label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]];
 
         if (rmUpperAddr == rmLowerAddr)
         {
@@ -480,6 +578,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
     initCoarseNeighb.setSize(0);
     coarseCoeffMap.setSize(0);
 
+
     // Create coarse-level coupled interfaces
 
     // Create coarse interfaces, addressing and coefficients
@@ -517,7 +616,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
             interfaceFields[intI].interface().initInternalFieldTransfer
             (
                 Pstream::blocking,
-                child_
+                agglomIndex_
             );
         }
     }
@@ -541,7 +640,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
                     fineInterface.internalFieldTransfer
                     (
                         Pstream::blocking,
-                        child_
+                        agglomIndex_
                     )
                 )
             );
@@ -563,7 +662,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
                 (
                     *coarseAddrPtr,
                     fineInterface,
-                    fineInterface.interfaceInternalField(child_),
+                    fineInterface.interfaceInternalField(agglomIndex_),
                     fineInterfaceAddr[intI]
                 ).ptr()
             );
@@ -664,7 +763,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
                 squareTypeField& activeCoarseLower = coarseLower.asSquare();
                 const squareTypeField& activeFineLower = fineLower.asSquare();
 
-                restrictResidual(fineDiag, coarseDiag);
+                restrictDiag(fineDiag, coarseDiag);
 
                 forAll(coeffRestrictAddr, fineCoeffI)
                 {
@@ -688,27 +787,41 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
             }
             else if (fineUpper.activeType() == blockCoeffBase::LINEAR)
             {
-                FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
-                    << "Matrix diagonal of square type and upper of linear type is not implemented"
+                FatalErrorIn
+                (
+                    "autoPtr<amgMatrix>"
+                    "BlockMatrixAgglomeration<Type>::restrictMatrix() const"
+                )   << "Matrix diagonal of square type and upper of "
+                    << "linear type is not implemented"
                     << abort(FatalError);
             }
             else
             {
-                FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
-                    << "Matrix diagonal of square type and upper of scalar type is not implemented"
+                FatalErrorIn
+                (
+                    "autoPtr<amgMatrix>"
+                    "BlockMatrixAgglomeration<Type>::restrictMatrix() const"
+                )   << "Matrix diagonal of square type and upper of "
+                    << "scalar type is not implemented"
                     << abort(FatalError);
             }
         }
         else if (fineDiag.activeType() == blockCoeffBase::LINEAR)
         {
-                FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
-                    << "Matrix diagonal of linear type not implemented"
+                FatalErrorIn
+                (
+                    "autoPtr<amgMatrix>"
+                    "BlockMatrixAgglomeration<Type>::restrictMatrix() const"
+                )   << "Matrix diagonal of linear type not implemented"
                     << abort(FatalError);
         }
         else
         {
-                FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
-                    << "Matrix diagonal of scalar type not implemented"
+                FatalErrorIn
+                (
+                    "autoPtr<amgMatrix>"
+                    "BlockMatrixAgglomeration<Type>::restrictMatrix() const"
+                )   << "Matrix diagonal of scalar type not implemented"
                     << abort(FatalError);
         }
     }
@@ -722,7 +835,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
                 squareTypeField& activeCoarseDiag = coarseDiag.asSquare();
                 const squareTypeField& activeFineUpper = fineUpper.asSquare();
 
-                restrictResidual(fineDiag, coarseDiag);
+                restrictDiag(fineDiag, coarseDiag);
 
                 forAll(coeffRestrictAddr, fineCoeffI)
                 {
@@ -730,27 +843,38 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
 
                     if (cCoeff >= 0)
                     {
-                        activeCoarseUpper[cCoeff]
-                            += activeFineUpper[fineCoeffI];
+                        activeCoarseUpper[cCoeff] +=
+                            activeFineUpper[fineCoeffI];
                     }
                     else
                     {
-                        // Add the fine face coefficient into the diagonal.
-                        activeCoarseDiag[-1 - cCoeff]
-                            += 2*activeFineUpper[fineCoeffI];
+                        // Add the fine face coefficient into the diagonal
+                        // Note: upper and lower coeffs are transpose of
+                        // each other.  HJ, 28/May/2014
+                        activeCoarseDiag[-1 - cCoeff] +=
+                            activeFineUpper[fineCoeffI]
+                          + activeFineUpper[fineCoeffI].T();
                     }
                 }
             }
             else if (fineUpper.activeType() == blockCoeffBase::LINEAR)
             {
-                FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
-                    << "Matrix diagonal of square type and upper of linear type is not implemented"
+                FatalErrorIn
+                (
+                    "autoPtr<amgMatrix>"
+                    "BlockMatrixAgglomeration<Type>::restrictMatrix() const"
+                )   << "Matrix diagonal of square type and upper of "
+                    << "linear type is not implemented"
                     << abort(FatalError);
             }
             else
             {
-                FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const")
-                    << "Matrix diagonal of square type and upper of scalar type is not implemented"
+                FatalErrorIn
+                (
+                    "autoPtr<amgMatrix>"
+                    "BlockMatrixAgglomeration<Type>::restrictMatrix() const"
+                )   << "Matrix diagonal of square type and upper of "
+                    << "scalar type is not implemented"
                     << abort(FatalError);
             }
 
@@ -775,48 +899,47 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
         }
     }
 
-    return autoPtr<BlockLduMatrix<Type> >
-    (
-        new BlockLduMatrix<Type>
-        (
-            coarseMatrix
-        )
-    );
+    return autoPtr<BlockLduMatrix<Type> >(coarseMatrixPtr);
 }
 
 
 template<class Type>
-void Foam::BlockMatrixAgglomeration<Type>::restrictResidual
+void Foam::BlockMatrixAgglomeration<Type>::restrictDiag
 (
-    const CoeffField<Type>& res,
-    CoeffField<Type>& coarseRes
+    const CoeffField<Type>& Coeff,
+    CoeffField<Type>& coarseCoeff
 ) const
 {
     typedef CoeffField<Type> TypeCoeffField;
 
-    if (res.activeType() == blockCoeffBase::SQUARE &&
-            coarseRes.activeType() == blockCoeffBase::SQUARE)
+    if
+    (
+        Coeff.activeType() == blockCoeffBase::SQUARE
+     && coarseCoeff.activeType() == blockCoeffBase::SQUARE
+    )
     {
+        typedef typename TypeCoeffField::squareType squareType;
         typedef typename TypeCoeffField::squareTypeField squareTypeField;
 
-        squareTypeField& activeCoarseRes = coarseRes.asSquare();
-        const squareTypeField& activeRes = res.asSquare();
+        squareTypeField& activeCoarseCoeff = coarseCoeff.asSquare();
+        const squareTypeField& activeCoeff = Coeff.asSquare();
 
-        // KRJ: 2013-02-05: To be done by correct pTraits?
-        forAll (coarseRes, i)
+        forAll (coarseCoeff, i)
         {
-            activeCoarseRes[i] *= 0.0;
+            activeCoarseCoeff[i] = pTraits<squareType>::zero;
         }
 
-        forAll (res, i)
+        forAll (Coeff, i)
         {
-            activeCoarseRes[child_[i]] += activeRes[i];
+            activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i];
         }
     }
     else
     {
-        FatalErrorIn("void  BlockMatrixAgglomeration<Type>::restrictResidual() const")
-            << "Only present for square type coeff type"
+        FatalErrorIn
+        (
+            "void  BlockMatrixAgglomeration<Type>::restrictDiag() const"
+        )   << "Only present for square type coeff type"
             << abort(FatalError);
     }
 }
@@ -833,7 +956,7 @@ void Foam::BlockMatrixAgglomeration<Type>::restrictResidual
 
     forAll (res, i)
     {
-        coarseRes[child_[i]] += res[i];
+        coarseRes[agglomIndex_[i]] += res[i];
     }
 }
 
@@ -847,7 +970,7 @@ void Foam::BlockMatrixAgglomeration<Type>::prolongateCorrection
 {
     forAll (x, i)
     {
-        x[i] += coarseX[child_[i]];
+        x[i] += coarseX[agglomIndex_[i]];
     }
 }
 
diff --git a/src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.H b/src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.H
index 1c2355e56..c7fdf9cf7 100644
--- a/src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.H
+++ b/src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.H
@@ -59,9 +59,6 @@ class BlockMatrixAgglomeration
 :
     public BlockMatrixCoarsening<Type>
 {
-
-public:
-
     // Private Data
 
         //- Reference to matrix
@@ -70,8 +67,9 @@ public:
         //- Reference to a templated norm calculator
         autoPtr<BlockCoeffNorm<Type> > normPtr_;
 
-        //- Child array
-        labelField child_;
+        //- Child array: for each fine equation give agglomeration cluster
+        //  index
+        labelField agglomIndex_;
 
         //- Group size
         label groupSize_;
@@ -82,6 +80,7 @@ public:
         //- Can a coarse level be constructed?
         bool coarsen_;
 
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -90,20 +89,23 @@ public:
         // Disallow default bitwise assignment
         void operator=(const BlockMatrixAgglomeration<Type>&);
 
-        //- Calculate child
-        void calcChild();
+        //- Calculate agglomeration index (child)
+        void calcAgglomeration();
 
 
     // Private Static Data
 
+        //- Diagonal scaling factor
+        static const scalar diagFactor_;
+
         //- Weighting factor
         static const scalar weightFactor_;
 
 
 public:
 
-        //- Runtime type information
-        TypeName("AAMG");
+    //- Runtime type information
+    TypeName("AAMG");
 
 
     // Constructors
@@ -117,6 +119,7 @@ public:
             const label minCoarseEqns
         );
 
+
     // Destructor
 
         virtual ~BlockMatrixAgglomeration();
@@ -131,8 +134,7 @@ public:
         }
 
         //- Restrict matrix
-        virtual autoPtr<BlockLduMatrix<Type> > restrictMatrix
-        ()  const;
+        virtual autoPtr<BlockLduMatrix<Type> > restrictMatrix() const;
 
         //- Restrict residual
         virtual void restrictResidual
@@ -141,11 +143,11 @@ public:
             Field<Type>& coarseRes
         ) const;
 
-        //- Restrict residual using CoeffField
-        void restrictResidual
+        //- Restrict CoeffField.  Used for diag coefficient
+        void restrictDiag
         (
-            const CoeffField<Type>& res,
-            CoeffField<Type>& coarseRes
+            const CoeffField<Type>& Coeff,
+            CoeffField<Type>& coarseCoeff
         ) const;
 
         //- Prolongate correction
@@ -154,7 +156,6 @@ public:
             Field<Type>& x,
             const Field<Type>& coarseX
         ) const;
-
 };
 
 
diff --git a/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C b/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C
index 5486864f7..b7a4f18f0 100644
--- a/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C
+++ b/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C
@@ -198,6 +198,7 @@ void Foam::coarseBlockAmgLevel<Type>::solve
 
     label maxIter = Foam::min(2*coarseningPtr_->minCoarseEqns(), 1000);
 
+    // Create artificial dictionary for top-level solution
     dictionary topLevelDict;
     topLevelDict.add("nDirections", "5");
     topLevelDict.add("minIter", 1);
@@ -233,15 +234,10 @@ void Foam::coarseBlockAmgLevel<Type>::solve
             matrixPtr_,
             topLevelDict
         ).solve(x, b);
-
-        if (BlockLduMatrix<Type>::debug >= 2)
-        {
-            coarseSolverPerf.print();
-        }
     }
     else
     {
-        topLevelDict.add("preconditioner", "ILU0");
+        topLevelDict.add("preconditioner", "Cholesky");
 
         coarseSolverPerf =
         //BlockGaussSeidelSolver<Type>
@@ -252,13 +248,30 @@ void Foam::coarseBlockAmgLevel<Type>::solve
             matrixPtr_,
             topLevelDict
         ).solve(x, b);
-
-        if (BlockLduMatrix<Type>::debug >= 2)
-        {
-            coarseSolverPerf.print();
-        }
     }
 
+    // Escape cases of top-level solver divergence
+    if
+    (
+        coarseSolverPerf.nIterations() == maxIter
+     && (
+            coarseSolverPerf.finalResidual()
+         >= coarseSolverPerf.initialResidual()
+        )
+    )
+    {
+        // Top-level solution failed.  Attempt rescue
+        // HJ, 27/Jul/2013
+        multiply(x, invDiag, b);
+
+        // Print top level correction failure as info for user
+        coarseSolverPerf.print();
+    }
+
+    if (BlockLduMatrix<Type>::debug >= 2)
+    {
+        coarseSolverPerf.print();
+    }
 }
 
 
diff --git a/src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C b/src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C
index 2f755fec1..a0a761e50 100644
--- a/src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C
+++ b/src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C
@@ -180,17 +180,26 @@ void Foam::fineBlockAmgLevel<Type>::solve
     const scalar relTol
 ) const
 {
+    Info<< "Fine level solver"<<endl;
 
-    Info<<"Fine level solver"<<endl;
+    // Create artificial dictionary for finest solution
+    dictionary finestDict;
+//     finestDict.add("nDirections", "5");
+    finestDict.add("minIter", 1);
+    finestDict.add("maxIter", 1000);
+    finestDict.add("tolerance", tolerance);
+    finestDict.add("relTol", relTol);
 
     if (matrix_.symmetric())
     {
+        finestDict.add("preconditioner", "Cholesky");
+
         BlockSolverPerformance<Type> coarseSolverPerf =
             BlockCGSolver<Type>
             (
                 "topLevelCorr",
                 matrix_,
-                dict_
+                finestDict
             ).solve(x, b);
 
         if (BlockLduMatrix<Type>::debug >= 2)
@@ -200,12 +209,14 @@ void Foam::fineBlockAmgLevel<Type>::solve
     }
     else
     {
+        finestDict.add("preconditioner", "Cholesky");
+
         BlockSolverPerformance<Type> coarseSolverPerf =
             BlockBiCGStabSolver<Type>
             (
                 "topLevelCorr",
                 matrix_,
-                dict_
+                finestDict
             ).solve(x, b);
 
         if (BlockLduMatrix<Type>::debug >= 2)
@@ -234,8 +245,8 @@ void Foam::fineBlockAmgLevel<Type>::scaleX
         x
     );
 
-    scalar scalingFactorNum = sumProd(x,b);
-    scalar scalingFactorDenom = sumProd(x,Ax);
+    scalar scalingFactorNum = sumProd(x, b);
+    scalar scalingFactorDenom = sumProd(x, Ax);
 
     vector scalingVector(scalingFactorNum, scalingFactorDenom, 0);
     reduce(scalingVector, sumOp<vector>());
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAmgSolver/BlockAmgSolver.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAmgSolver/BlockAmgSolver.C
index fd2551c3f..c3c6d0c09 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAmgSolver/BlockAmgSolver.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAmgSolver/BlockAmgSolver.C
@@ -31,12 +31,8 @@ Author
 
 #include "BlockAmgSolver.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-//- Construct from matrix and solver data stream
 template<class Type>
 Foam::BlockAmgSolver<Type>::BlockAmgSolver
 (
@@ -76,34 +72,22 @@ Foam::BlockAmgSolver<Type>::solve
         this->fieldName()
     );
 
-    // Create local references to avoid the spread this-> ugliness
-    const BlockLduMatrix<Type>& matrix = this->matrix_;
-
     scalar norm = this->normFactor(x, b);
 
-    Field<Type> wA(x.size());
-
-    // Calculate residual.  Note: sign of residual swapped for efficiency
-    matrix.Amul(wA, x);
-    wA -= b;
-
-    solverPerf.initialResidual() = gSum(cmptMag(wA))/norm;
+    // Calculate initial residual
+    solverPerf.initialResidual() = gSum(cmptMag(amg_.residual(x, b)))/norm;
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
-    if (!solverPerf.checkConvergence(this->tolerance(), this->relTolerance()))
+    if (!this->stop(solverPerf))
     {
         do
         {
             amg_.cycle(x, b);
 
-            // Re-calculate residual.  Note: sign of residual swapped
-            // for efficiency
-            matrix.Amul(wA, x);
-            wA -= b;
+            solverPerf.finalResidual() =
+                gSum(cmptMag(amg_.residual(x, b)))/norm;
 
-            solverPerf.finalResidual() = gSum(cmptMag(wA))/norm;
             solverPerf.nIterations()++;
-
         } while (!this->stop(solverPerf));
     }
 
diff --git a/src/lduSolvers/amg/amgPolicy/pamgPolicy.C b/src/lduSolvers/amg/amgPolicy/pamgPolicy.C
index df4771f56..80cde8871 100644
--- a/src/lduSolvers/amg/amgPolicy/pamgPolicy.C
+++ b/src/lduSolvers/amg/amgPolicy/pamgPolicy.C
@@ -29,6 +29,7 @@ Description
 
 Author
     Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
 \*---------------------------------------------------------------------------*/
 
 #include "pamgPolicy.H"
diff --git a/src/lduSolvers/amg/coarseAmgLevel.C b/src/lduSolvers/amg/coarseAmgLevel.C
index 8c01aa767..b91c384ca 100644
--- a/src/lduSolvers/amg/coarseAmgLevel.C
+++ b/src/lduSolvers/amg/coarseAmgLevel.C
@@ -237,7 +237,7 @@ void Foam::coarseAmgLevel::solve
         ).solve(x, b, cmpt);
     }
 
-    // Escape cases of solver divergence
+    // Escape cases of top-level solver divergence
     if
     (
         coarseSolverPerf.nIterations() == maxIter