diff --git a/applications/test/PisoFoam/PisoFoam.C b/applications/test/PisoFoam/PisoFoam.C
index 67f0084e7192f7024446e7d8994b9a3247b5060a..1713e9747c77b09aa357ff547f9bf34743f40ea4 100644
--- a/applications/test/PisoFoam/PisoFoam.C
+++ b/applications/test/PisoFoam/PisoFoam.C
@@ -102,7 +102,7 @@ int main(int argc, char *argv[])
                         "{"
                         "    /*solver          SmoothSolver;*/"
                         "    smoother        GaussSeidel;"
-                        "    solver           PBiCICG;"
+                        "    solver           PBiCCCG;"
                         "    preconditioner   none;"
                         "    tolerance        (1e-7 1e-7 1);"
                         "    relTol           (0 0 0);"
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
index 97b2b2f40e9e34dca76899963571b3858f63eca6..00782ac75b19b49225873331122a6aa0159da50f 100644
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
+++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
@@ -382,7 +382,7 @@ Foam::Ostream& Foam::operator<<
 #include "LduMatrixOperations.C"
 #include "LduMatrixATmul.C"
 #include "LduMatrixUpdateMatrixInterfaces.C"
-#include "LduMatrixTests.C"
+#include "SolverPerformance.C"
 #include "LduMatrixPreconditioner.C"
 #include "LduMatrixSmoother.C"
 #include "LduMatrixSolver.C"
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
index e0ce6d28ba45f25bafe12b92bfaa1ae89ef11c63..507fd40e94f6fe988325699f92e6a4bdbfde44c9 100644
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
+++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
@@ -54,6 +54,7 @@ SourceFiles
 #include "Field.H"
 #include "FieldField.H"
 #include "LduInterfaceFieldPtrsList.H"
+#include "SolverPerformance.H"
 #include "typeInfo.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
@@ -75,27 +76,6 @@ Ostream& operator<<
     const LduMatrix<Type, DType, LUType>&
 );
 
-template<class Type, class DType, class LUType>
-typename LduMatrix<Type, DType, LUType>::solverPerformance max
-(
-    const typename LduMatrix<Type, DType, LUType>::solverPerformance&,
-    const typename LduMatrix<Type, DType, LUType>::solverPerformance&
-);
-
-template<class Type, class DType, class LUType>
-Istream& operator>>
-(
-    Istream&,
-    typename LduMatrix<Type, DType, LUType>::solverPerformance&
-);
-
-template<class Type, class DType, class LUType>
-Ostream& operator<<
-(
-    Ostream&,
-    const typename LduMatrix<Type, DType, LUType>::solverPerformance&
-);
-
 
 /*---------------------------------------------------------------------------*\
                            Class LduMatrix Declaration
@@ -127,168 +107,7 @@ class LduMatrix
 
 public:
 
-    //- Class returned by the solver
-    //  containing performance statistics
-    class solverPerformance
-    {
-        word   solverName_;
-        word   fieldName_;
-        Type   initialResidual_;
-        Type   finalResidual_;
-        label  noIterations_;
-        bool   converged_;
-        FixedList<bool, pTraits<Type>::nComponents> singular_;
-
-
-    public:
-
-        // Constructors
-
-            solverPerformance()
-            :
-                initialResidual_(pTraits<Type>::zero),
-                finalResidual_(pTraits<Type>::zero),
-                noIterations_(0),
-                converged_(false),
-                singular_(false)
-            {}
-
-
-            solverPerformance
-            (
-                const word&  solverName,
-                const word&  fieldName,
-                const Type&  iRes = pTraits<Type>::zero,
-                const Type&  fRes = pTraits<Type>::zero,
-                const label  nIter = 0,
-                const bool   converged = false,
-                const bool   singular = false
-            )
-            :
-                solverName_(solverName),
-                fieldName_(fieldName),
-                initialResidual_(iRes),
-                finalResidual_(fRes),
-                noIterations_(nIter),
-                converged_(converged),
-                singular_(singular)
-            {}
-
-
-        // Member functions
-
-            //- Return solver name
-            const word& solverName() const
-            {
-                return solverName_;
-            }
-
-            //- Return solver name
-            word& solverName()
-            {
-                return solverName_;
-            }
-
-
-            //- Return field name
-            const word& fieldName() const
-            {
-                return fieldName_;
-            }
-
-
-            //- Return initial residual
-            const Type& initialResidual() const
-            {
-                return initialResidual_;
-            }
-
-            //- Return initial residual
-            Type& initialResidual()
-            {
-                return initialResidual_;
-            }
-
-
-            //- Return final residual
-            const Type& finalResidual() const
-            {
-                return finalResidual_;
-            }
-
-            //- Return final residual
-            Type& finalResidual()
-            {
-                return finalResidual_;
-            }
-
-
-            //- Return number of iterations
-            label nIterations() const
-            {
-                return noIterations_;
-            }
-
-            //- Return number of iterations
-            label& nIterations()
-            {
-                return noIterations_;
-            }
-
-
-            //- Has the solver converged?
-            bool converged() const
-            {
-                return converged_;
-            }
-
-            //- Is the matrix singular?
-            bool singular() const;
-
-            //- Check, store and return convergence
-            bool converged
-            (
-                const Type& tolerance,
-                const Type& relTolerance
-            );
-
-            //- Singularity test
-            bool checkSingularity(const Type& residual);
-
-            //- Print summary of solver performance to the given stream
-            void print(Ostream& os) const;
-
-
-        // Member Operators
-
-            bool operator!=(const solverPerformance&) const;
-
-
-        // Friend functions
-
-            //- Return the element-wise maximum of two solverPerformances
-            friend solverPerformance max <Type, DType, LUType>
-            (
-                const solverPerformance&,
-                const solverPerformance&
-            );
-
-
-        // Ostream Operator
-
-            friend Istream& operator>> <Type, DType, LUType>
-            (
-                Istream&,
-                solverPerformance&
-            );
-
-            friend Ostream& operator<< <Type, DType, LUType>
-            (
-                Ostream&,
-                const solverPerformance&
-            );
-    };
-
+    friend class SolverPerformance<Type>;
 
     //- Abstract base-class for LduMatrix solvers
     class solver
@@ -417,7 +236,7 @@ public:
             //- Read and reset the solver parameters from the given dictionary
             virtual void read(const dictionary& solverDict);
 
-            virtual solverPerformance solve
+            virtual SolverPerformance<Type> solve
             (
                 Field<Type>& psi
             ) const = 0;
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.C
similarity index 72%
rename from src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C
rename to src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.C
index f5a675f0440eae2d0e70fd6b456a8ac6fe2c9d2f..a50a6e537ff237bc905354d81bc1ca744ae31ced 100644
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C
+++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.C
@@ -27,23 +27,24 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class Type, class DType, class LUType>
-bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::checkSingularity
+template<class Type>
+bool Foam::SolverPerformance<Type>::checkSingularity
 (
     const Type& wApA
 )
 {
     for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
     {
-        singular_[cmpt] = component(wApA, cmpt) < vsmall_;
+        singular_[cmpt] =
+            component(wApA, cmpt) < LduMatrix<scalar, scalar, scalar>::vsmall_;
     }
 
     return singular();
 }
 
 
-template<class Type, class DType, class LUType>
-bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular() const
+template<class Type>
+bool Foam::SolverPerformance<Type>::singular() const
 {
     for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
     {
@@ -54,14 +55,14 @@ bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular() const
 }
 
 
-template<class Type, class DType, class LUType>
-bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::converged
+template<class Type>
+bool Foam::SolverPerformance<Type>::converged
 (
     const Type& Tolerance,
     const Type& RelTolerance
 )
 {
-    if (debug >= 2)
+    if (LduMatrix<scalar, scalar, scalar>::debug >= 2)
     {
         Info<< solverName_
             << ":  Iteration " << noIterations_
@@ -73,7 +74,8 @@ bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::converged
     (
         finalResidual_ < Tolerance
      || (
-            RelTolerance > small_*pTraits<Type>::one
+            RelTolerance
+          > LduMatrix<scalar, scalar, scalar>::small_*pTraits<Type>::one
          && finalResidual_ < cmptMultiply(RelTolerance, initialResidual_)
         )
     )
@@ -89,8 +91,8 @@ bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::converged
 }
 
 
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::solverPerformance::print
+template<class Type>
+void Foam::SolverPerformance<Type>::print
 (
     Ostream& os
 ) const
@@ -115,10 +117,10 @@ void Foam::LduMatrix<Type, DType, LUType>::solverPerformance::print
 }
 
 
-template<class Type, class DType, class LUType>
-bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::operator!=
+template<class Type>
+bool Foam::SolverPerformance<Type>::operator!=
 (
-    const LduMatrix<Type, DType, LUType>::solverPerformance& sp
+    const SolverPerformance<Type>& sp
 ) const
 {
     return
@@ -134,14 +136,14 @@ bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::operator!=
 }
 
 
-template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance Foam::max
+template<class Type>
+typename Foam::SolverPerformance<Type> Foam::max
 (
-    const typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp1,
-    const typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp2
+    const typename Foam::SolverPerformance<Type>& sp1,
+    const typename Foam::SolverPerformance<Type>& sp2
 )
 {
-    return lduMatrix::solverPerformance
+    return SolverPerformance<Type>
     (
         sp1.solverName(),
         sp1.fieldName_,
@@ -154,14 +156,14 @@ typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance Foam::max
 }
 
 
-template<class Type, class DType, class LUType>
+template<class Type>
 Foam::Istream& Foam::operator>>
 (
     Istream& is,
-    typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp
+    typename Foam::SolverPerformance<Type>& sp
 )
 {
-    is.readBeginList("LduMatrix<Type, DType, LUType>::solverPerformance");
+    is.readBeginList("SolverPerformance<Type>");
     is  >> sp.solverName_
         >> sp.fieldName_
         >> sp.initialResidual_
@@ -169,17 +171,17 @@ Foam::Istream& Foam::operator>>
         >> sp.noIterations_
         >> sp.converged_
         >> sp.singular_;
-    is.readEndList("LduMatrix<Type, DType, LUType>::solverPerformance");
+    is.readEndList("SolverPerformance<Type>");
 
     return is;
 }
 
 
-template<class Type, class DType, class LUType>
+template<class Type>
 Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
-    const typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp
+    const typename Foam::SolverPerformance<Type>& sp
 )
 {
     os  << token::BEGIN_LIST
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.H
new file mode 100644
index 0000000000000000000000000000000000000000..fdfd269e17f49ace816532a2b9830ceb3fc14c5c
--- /dev/null
+++ b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.H
@@ -0,0 +1,258 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::SolverPerformance
+
+Description
+    SolverPerformance is a general matrix class in which the coefficients are
+    stored as three arrays, one for the upper triangle, one for the
+    lower triangle and a third for the diagonal.
+
+SourceFiles
+    SolverPerformance.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SolverPerformance_H
+#define SolverPerformance_H
+
+#include "word.H"
+#include "FixedList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class Type>
+class SolverPerformance;
+
+template<class Type>
+SolverPerformance<Type> max
+(
+    const SolverPerformance<Type>&,
+    const SolverPerformance<Type>&
+);
+
+template<class Type>
+Istream& operator>>
+(
+    Istream&,
+    SolverPerformance<Type>&
+);
+
+template<class Type>
+Ostream& operator<<
+(
+    Ostream&,
+    const SolverPerformance<Type>&
+);
+
+
+
+/*---------------------------------------------------------------------------*\
+                       Class SolverPerformance Declaration
+\*---------------------------------------------------------------------------*/
+
+//- Class returned by the solver
+//  containing performance statistics
+template<class Type>
+class SolverPerformance
+{
+    word   solverName_;
+    word   fieldName_;
+    Type   initialResidual_;
+    Type   finalResidual_;
+    label  noIterations_;
+    bool   converged_;
+    FixedList<bool, pTraits<Type>::nComponents> singular_;
+
+
+public:
+
+    // Constructors
+
+        SolverPerformance()
+        :
+            initialResidual_(pTraits<Type>::zero),
+            finalResidual_(pTraits<Type>::zero),
+            noIterations_(0),
+            converged_(false),
+            singular_(false)
+        {}
+
+
+        SolverPerformance
+        (
+            const word&  solverName,
+            const word&  fieldName,
+            const Type&  iRes = pTraits<Type>::zero,
+            const Type&  fRes = pTraits<Type>::zero,
+            const label  nIter = 0,
+            const bool   converged = false,
+            const bool   singular = false
+        )
+        :
+            solverName_(solverName),
+            fieldName_(fieldName),
+            initialResidual_(iRes),
+            finalResidual_(fRes),
+            noIterations_(nIter),
+            converged_(converged),
+            singular_(singular)
+        {}
+
+
+    // Member functions
+
+        //- Return solver name
+        const word& solverName() const
+        {
+            return solverName_;
+        }
+
+        //- Return solver name
+        word& solverName()
+        {
+            return solverName_;
+        }
+
+
+        //- Return field name
+        const word& fieldName() const
+        {
+            return fieldName_;
+        }
+
+
+        //- Return initial residual
+        const Type& initialResidual() const
+        {
+            return initialResidual_;
+        }
+
+        //- Return initial residual
+        Type& initialResidual()
+        {
+            return initialResidual_;
+        }
+
+
+        //- Return final residual
+        const Type& finalResidual() const
+        {
+            return finalResidual_;
+        }
+
+        //- Return final residual
+        Type& finalResidual()
+        {
+            return finalResidual_;
+        }
+
+
+        //- Return number of iterations
+        label nIterations() const
+        {
+            return noIterations_;
+        }
+
+        //- Return number of iterations
+        label& nIterations()
+        {
+            return noIterations_;
+        }
+
+
+        //- Has the solver converged?
+        bool converged() const
+        {
+            return converged_;
+        }
+
+        //- Is the matrix singular?
+        bool singular() const;
+
+        //- Check, store and return convergence
+        bool converged
+        (
+            const Type& tolerance,
+            const Type& relTolerance
+        );
+
+        //- Singularity test
+        bool checkSingularity(const Type& residual);
+
+        //- Print summary of solver performance to the given stream
+        void print(Ostream& os) const;
+
+
+    // Member Operators
+
+        bool operator!=(const SolverPerformance<Type>&) const;
+
+
+    // Friend functions
+
+        //- Return the element-wise maximum of two SolverPerformance<Type>s
+        friend SolverPerformance<Type> max <Type>
+        (
+            const SolverPerformance<Type>&,
+            const SolverPerformance<Type>&
+        );
+
+
+    // Ostream Operator
+
+        friend Istream& operator>> <Type>
+        (
+            Istream&,
+            SolverPerformance<Type>&
+        );
+
+        friend Ostream& operator<< <Type>
+        (
+            Ostream&,
+            const SolverPerformance<Type>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+//#   include "SolverPerformance.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C
index af128d249a8fab2e8617dfa8ee71c52ab60beee9..adbe9113a271c750af850a9443baf3190840d232 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C
@@ -55,7 +55,7 @@ void Foam::DiagonalSolver<Type, DType, LUType>::read
 
 
 template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
+Foam::SolverPerformance<Type>
 Foam::DiagonalSolver<Type, DType, LUType>::solve
 (
     Field<Type>& psi
@@ -63,7 +63,7 @@ Foam::DiagonalSolver<Type, DType, LUType>::solve
 {
     psi = this->matrix_.source()/this->matrix_.diag();
 
-    return typename LduMatrix<Type, DType, LUType>::solverPerformance
+    return SolverPerformance<Type>
     (
         typeName,
         this->fieldName_,
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H
index 400891b97776681ec9849d007605785213d9b79f..b9dc537c08afaf1ea2a771046842f84dc776aa56 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H
@@ -83,10 +83,7 @@ public:
         void read(const dictionary& solverDict);
 
         //- Solve the matrix with this solver
-        typename LduMatrix<Type, DType, LUType>::solverPerformance solve
-        (
-            Field<Type>& psi
-        ) const;
+        virtual SolverPerformance<Type> solve(Field<Type>& psi) const;
 };
 
 
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C
index 45453ff1d062095d2c3e08f23e467d6887517555..ab0c1f4fc95f3f28a663d2f3470d6b15836085c6 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C
@@ -47,7 +47,7 @@ Foam::PBiCCCG<Type, DType, LUType>::PBiCCCG
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
+Foam::SolverPerformance<Type>
 Foam::PBiCCCG<Type, DType, LUType>::solve
 (
     Field<Type>& psi
@@ -56,7 +56,7 @@ Foam::PBiCCCG<Type, DType, LUType>::solve
     word preconditionerName(this->controlDict_.lookup("preconditioner"));
 
     // --- Setup class containing solver performance data
-    typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
+    SolverPerformance<Type> solverPerf
     (
         preconditionerName + typeName,
         this->fieldName_
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.H
index e5363dc6510a33e8463d09173a7e56d286405143..39c0e56d44768d12e13879c872c6a97035db9fe9 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.H
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.H
@@ -87,10 +87,7 @@ public:
     // Member Functions
 
         //- Solve the matrix with this solver
-        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
-        (
-            Field<Type>& psi
-        ) const;
+        virtual SolverPerformance<Type> solve(Field<Type>& psi) const;
 };
 
 
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.C
index b3665f0617565e4b743c18b326ca1dfabee9fa74..c92499c9fc2d2cd6b424447e026c9862ff611cc6 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.C
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.C
@@ -47,13 +47,13 @@ Foam::PBiCICG<Type, DType, LUType>::PBiCICG
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
+Foam::SolverPerformance<Type>
 Foam::PBiCICG<Type, DType, LUType>::solve(Field<Type>& psi) const
 {
     word preconditionerName(this->controlDict_.lookup("preconditioner"));
 
     // --- Setup class containing solver performance data
-    typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
+    SolverPerformance<Type> solverPerf
     (
         preconditionerName + typeName,
         this->fieldName_
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.H
index 340ae486678db5e6f1291734753c36b4945200e0..1265af2aac4b6197e9c7e2d7d0038ccb846a4a5a 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.H
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.H
@@ -87,10 +87,7 @@ public:
     // Member Functions
 
         //- Solve the matrix with this solver
-        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
-        (
-            Field<Type>& psi
-        ) const;
+        virtual SolverPerformance<Type> solve(Field<Type>& psi) const;
 };
 
 
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C
index e6206e61163c47dc88cf1739bba98e4c03f267ba..83abd1f24996916dec2b3d9a7b5a45484d4866d4 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C
@@ -47,13 +47,13 @@ Foam::PCICG<Type, DType, LUType>::PCICG
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
+typename Foam::SolverPerformance<Type>
 Foam::PCICG<Type, DType, LUType>::solve(Field<Type>& psi) const
 {
     word preconditionerName(this->controlDict_.lookup("preconditioner"));
 
     // --- Setup class containing solver performance data
-    typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
+    SolverPerformance<Type> solverPerf
     (
         preconditionerName + typeName,
         this->fieldName_
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.H
index 250cc41b9a363dea56420fedfc723bd7ad86b019..cafeb7f61d3a79ae55861e4dea79b68942f83c38 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.H
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.H
@@ -87,10 +87,7 @@ public:
     // Member Functions
 
         //- Solve the matrix with this solver
-        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
-        (
-            Field<Type>& psi
-        ) const;
+        virtual SolverPerformance<Type> solve(Field<Type>& psi) const;
 };
 
 
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C
index e6f0241308e5c2377c6e38b55f58693f370941bd..56ff3675e5b9ecfb09a2eff45344462679516c3b 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C
@@ -58,11 +58,11 @@ void Foam::SmoothSolver<Type, DType, LUType>::readControls()
 
 
 template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
+Foam::SolverPerformance<Type>
 Foam::SmoothSolver<Type, DType, LUType>::solve(Field<Type>& psi) const
 {
     // --- Setup class containing solver performance data
-    typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
+    SolverPerformance<Type> solverPerf
     (
         typeName,
         this->fieldName_
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H
index 3d789420b390e3b67710f01fc9d0fba7ed6585fe..54db11be8c569bc80530e486f9b7cada29806ff8 100644
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H
+++ b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H
@@ -86,10 +86,7 @@ public:
     // Member Functions
 
         //- Solve the matrix with this solver
-        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
-        (
-            Field<Type>& psi
-        ) const;
+        virtual SolverPerformance<Type> solve(Field<Type>& psi) const;
 };
 
 
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
index f567254f8655c1b5762fc6450297a93555d5a851..8b72617d0469549901e041790af07ac5f2733270 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
@@ -387,10 +387,18 @@ public:
             //  Solver controls read from fvSolution
             autoPtr<fvSolver> solver();
 
-            //- Solve returning the solution statistics.
+            //- Solve segregated or coupled returning the solution statistics.
             //  Use the given solver controls
             lduMatrix::solverPerformance solve(const dictionary&);
 
+            //- Solve segregated returning the solution statistics.
+            //  Use the given solver controls
+            lduMatrix::solverPerformance solveSegregated(const dictionary&);
+
+            //- Solve coupled returning the solution statistics.
+            //  Use the given solver controls
+            lduMatrix::solverPerformance solveCoupled(const dictionary&);
+
             //- Solve returning the solution statistics.
             //  Solver controls read from fvSolution
             lduMatrix::solverPerformance solve();
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
index bc3fc66d018ee99dfd90359723478396b76a9514..b37bfa2d5b7d117d8ca5838dff129d794396388f 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
@@ -23,6 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "LduMatrix.H"
+#include "diagTensorField.H"
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
@@ -62,12 +65,51 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
             << endl;
     }
 
+    word type(solverControls.lookupOrDefault<word>("type", "segregated"));
+
+    if (type == "segregated")
+    {
+        return solveSegregated(solverControls);
+    }
+    else if (type == "coupled")
+    {
+        return solveCoupled(solverControls);
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "fvMatrix<Type>::solve(const dictionary& solverControls)",
+            solverControls
+        )   << "Unknown type " << type
+            << "; currently supported solver types are segregated and coupled"
+            << exit(FatalIOError);
+
+        return lduMatrix::solverPerformance();
+    }
+}
+
+
+template<class Type>
+Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solveSegregated
+(
+    const dictionary& solverControls
+)
+{
+    if (debug)
+    {
+        Info<< "fvMatrix<Type>::solveSegregated"
+               "(const dictionary& solverControls) : "
+               "solving fvMatrix<Type>"
+            << endl;
+    }
+
     GeometricField<Type, fvPatchField, volMesh>& psi =
        const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);
 
     lduMatrix::solverPerformance solverPerfVec
     (
-        "fvMatrix<Type>::solve",
+        "fvMatrix<Type>::solveSegregated",
         psi.name()
     );
 
@@ -164,6 +206,62 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
 }
 
 
+template<class Type>
+Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solveCoupled
+(
+    const dictionary& solverControls
+)
+{
+    if (debug)
+    {
+        Info<< "fvMatrix<Type>::solveCoupled"
+               "(const dictionary& solverControls) : "
+               "solving fvMatrix<Type>"
+            << endl;
+    }
+
+    GeometricField<Type, fvPatchField, volMesh>& psi =
+       const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);
+
+    LduMatrix<Type, scalar, scalar> coupledMatrix(psi.mesh());
+    coupledMatrix.diag() = diag();
+    coupledMatrix.upper() = upper();
+    coupledMatrix.lower() = lower();
+    coupledMatrix.source() = source();
+
+    addBoundaryDiag(coupledMatrix.diag(), 0);
+    addBoundarySource(coupledMatrix.source(), false);
+
+    coupledMatrix.interfaces() = psi.boundaryField().interfaces();
+    coupledMatrix.interfacesUpper() = boundaryCoeffs().component(0);
+    coupledMatrix.interfacesLower() = internalCoeffs().component(0);
+
+    autoPtr<typename LduMatrix<Type, scalar, scalar>::solver>
+    coupledMatrixSolver
+    (
+        LduMatrix<Type, scalar, scalar>::solver::New
+        (
+            psi.name(),
+            coupledMatrix,
+            solverControls
+        )
+    );
+
+    SolverPerformance<Type> solverPerf
+    (
+        coupledMatrixSolver->solve(psi)
+    );
+
+    solverPerf.print(Info);
+
+    psi.correctBoundaryConditions();
+
+    // psi.mesh().setSolverPerformance(psi.name(), solverPerf);
+
+    return lduMatrix::solverPerformance();
+}
+
+
 template<class Type>
 Foam::autoPtr<typename Foam::fvMatrix<Type>::fvSolver>
 Foam::fvMatrix<Type>::solver()
diff --git a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
index 23a259d574552e48cd73d2214ce8628122ad1182..aed610cd4c2a862f985eb0d2e3be7bf2f7268f27 100644
--- a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
@@ -130,14 +130,15 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::fvSolver::solve
 
 
 template<>
-Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
+Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solveSegregated
 (
     const dictionary& solverControls
 )
 {
     if (debug)
     {
-        Info<< "fvMatrix<scalar>::solve(const dictionary& solverControls) : "
+        Info<< "fvMatrix<scalar>::solveSegregated"
+               "(const dictionary& solverControls) : "
                "solving fvMatrix<scalar>"
             << endl;
     }
diff --git a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H
index 072c9f2cd09e0de014046fb5dfa78b41546bd239..29c51a155b307362f7e8e92faa1ea06af865b8f5 100644
--- a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H
+++ b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H
@@ -68,7 +68,7 @@ lduMatrix::solverPerformance fvMatrix<scalar>::fvSolver::solve
 );
 
 template<>
-lduMatrix::solverPerformance fvMatrix<scalar>::solve
+lduMatrix::solverPerformance fvMatrix<scalar>::solveSegregated
 (
     const dictionary&
 );