diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 26bb4865feae54c4f470be573cfaad3b20fb313b..3a49df4e218d64d926bb8f3554a994e13e3007fb 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -256,11 +256,6 @@ $(pairGAMGAgglomeration)/pairGAMGAgglomerationCombineLevels.C
 algebraicPairGAMGAgglomeration = $(GAMGAgglomerations)/algebraicPairGAMGAgglomeration
 $(algebraicPairGAMGAgglomeration)/algebraicPairGAMGAgglomeration.C
 
-matrices/LduMatrix/LduMatrix/lduMatrices.C
-matrices/LduMatrix/Preconditioners/lduPreconditioners.C
-matrices/LduMatrix/Smoothers/lduSmoothers.C
-matrices/LduMatrix/Solvers/lduSolvers.C
-
 meshes/lduMesh/lduMesh.C
 
 primitiveShapes = meshes/primitiveShapes
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/CyclicLduInterfaceField.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/CyclicLduInterfaceField.C
deleted file mode 100644
index 9fa8f795821a7581c65ca09b9862dfa64fdee217..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/CyclicLduInterfaceField.C
+++ /dev/null
@@ -1,66 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "CyclicLduInterfaceField.H"
-#include "diagTensorField.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    //defineTypeNameAndDebug(CyclicLduInterfaceField, 0);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::CyclicLduInterfaceField<Type>::~CyclicLduInterfaceField()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-void Foam::CyclicLduInterfaceField<Type>::transformCoupleField
-(
-    Field<Type>& f
-) const
-{
-    if (doTransform())
-    {
-        label sizeby2 = f.size()/2;
-
-        for (label facei=0; facei<sizeby2; facei++)
-        {
-            f[facei] = transform(f[facei], forwardT()[0]);
-            f[facei + sizeby2] = transform(f[facei + sizeby2], forwardT()[0]);
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/CyclicLduInterfaceField.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/CyclicLduInterfaceField.H
deleted file mode 100644
index c783d6916653abc3da291c0b239dbbfa9cda5956..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/CyclicLduInterfaceField.H
+++ /dev/null
@@ -1,103 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::CyclicLduInterfaceField
-
-Description
-    Abstract base class for cyclic coupled interfaces.
-
-SourceFiles
-    CyclicLduInterfaceField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef CyclicLduInterfaceField_H
-#define CyclicLduInterfaceField_H
-
-#include "primitiveFieldsFwd.H"
-#include "typeInfo.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                  Class CyclicLduInterfaceField Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type>
-class CyclicLduInterfaceField
-{
-
-public:
-
-    //- Runtime type information
-    TypeName("CyclicLduInterfaceField");
-
-
-    // Constructors
-
-        //- Construct given coupled patch
-        CyclicLduInterfaceField()
-        {}
-
-
-    // Destructor
-
-        virtual ~CyclicLduInterfaceField();
-
-
-    // Member Functions
-
-        // Access
-
-            //- Is the transform required
-            virtual bool doTransform() const = 0;
-
-            //- Return face transformation tensor
-            virtual const tensorField& forwardT() const = 0;
-
-            //- Return neighbour-cell transformation tensor
-            virtual const tensorField& reverseT() const = 0;
-
-            //- Return rank of component for transform
-            virtual int rank() const = 0;
-
-
-        //- Transform given patch internal field
-        void transformCoupleField(Field<Type>& f) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.C
deleted file mode 100644
index 9f3b6adc7284e38f4db69a02a6d0d69b88e44453..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.C
+++ /dev/null
@@ -1,48 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "lduInterfaceField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(lduInterfaceField, 0);
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-lduInterfaceField::~lduInterfaceField()
-{}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.H
deleted file mode 100644
index ac4752a49fbe2fe6f84800fdeb80a0aa18cfc9dc..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.H
+++ /dev/null
@@ -1,147 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::LduInterfaceField
-
-Description
-    An abstract base class for implicitly-coupled interface fields
-    e.g. processor and cyclic patch fields.
-
-SourceFiles
-    LduInterfaceField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef LduInterfaceField_H
-#define LduInterfaceField_H
-
-#include "lduInterface.H"
-#include "Field.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                     Class LduInterfaceField Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type>
-class LduInterfaceField
-{
-    // Private data
-
-        //- Reference to the coupled patch this field is defined for
-        const lduInterface& interface_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        LduInterfaceField(const LduInterfaceField&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const LduInterfaceField&);
-
-
-public:
-
-        class Amultiplier
-        {
-        public:
-            
-            Amultiplier()
-            {}
-
-            virtual void addAmul
-            (
-                Field<Type>& Apsi,
-                const Field<Type>& psi
-            ) const = 0;
-        };
-
-
-    //- Runtime type information
-    TypeName("LduInterfaceField");
-
-
-    // Constructors
-
-        //- Construct given interface
-        LduInterfaceField(const lduInterface& interface)
-        :
-            interface_(interface)
-        {}
-
-
-    // Destructor
-
-        virtual ~LduInterfaceField();
-
-
-    // Member Functions
-
-        // Access
-
-            //- Return the interface
-            const lduInterface& interface() const
-            {
-                return interface_;
-            }
-
-
-        // Coupled interface matrix update
-
-            //- Initialise neighbour matrix update
-            virtual void initInterfaceMatrixUpdate
-            (
-                Field<Type>& Apsi,
-                const Field<Type>& psi,
-                const Amultiplier&,
-                const Pstream::commsTypes commsType
-            ) const
-            {}
-
-            //- Update result field based on interface functionality
-            virtual void updateInterfaceMatrix
-            (
-                Field<Type>& Apsi,
-                const Field<Type>& psi,
-                const Amultiplier&,
-                const Pstream::commsTypes commsType
-            ) const = 0;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceFieldPtrsList.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceFieldPtrsList.H
deleted file mode 100644
index 4dbe723afaaab3eda6a9d2e1917829df1b69c387..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceFieldPtrsList.H
+++ /dev/null
@@ -1,71 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Type
-    lduInterfaceFieldPtrsList
-
-Description
-    List of coupled interface fields to be used in coupling.
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef LduInterfaceFieldPtrsList_H
-#define LduInterfaceFieldPtrsList_H
-
-#include "LduInterfaceField.H"
-#include "UPtrList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                   Class LduInterfaceFieldPtrsList Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type>
-class LduInterfaceFieldPtrsList
-:
-    public UPtrList<const LduInterfaceField<Type> >
-{
-public:
-
-    LduInterfaceFieldPtrsList(label size)
-    :
-        UPtrList<const LduInterfaceField<Type> >(size)
-    {}
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
-
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/ProcessorLduInterfaceField.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/ProcessorLduInterfaceField.C
deleted file mode 100644
index c0a3aed5a061f04cd256467c38273f3875aea385..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/ProcessorLduInterfaceField.C
+++ /dev/null
@@ -1,67 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "processorLduInterfaceField.H"
-#include "diagTensorField.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    //defineTypeNameAndDebug(processorLduInterfaceField, 0);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::processorLduInterfaceField<Type>::~processorLduInterfaceField()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-void Foam::processorLduInterfaceField<Type>::transformCoupleField
-(
-    Field<Type>& f
-) const
-{
-    if (doTransform())
-    {
-        if (forwardT().size() == 1)
-        {
-            transform(f, forwardT()[0], f);
-        }
-        else
-        {
-            transform(f, forwardT(), f);
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/ProcessorLduInterfaceField.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/ProcessorLduInterfaceField.H
deleted file mode 100644
index 885ba210e4c96cfd2c918c581e8119ed2c544c63..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/ProcessorLduInterfaceField.H
+++ /dev/null
@@ -1,106 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::ProcessorLduInterfaceField
-
-Description
-    Abstract base class for processor coupled interfaces.
-
-SourceFiles
-    ProcessorLduInterfaceField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef ProcessorLduInterfaceField_H
-#define ProcessorLduInterfaceField_H
-
-#include "primitiveFieldsFwd.H"
-#include "typeInfo.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                  Class ProcessorLduInterfaceField Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type>
-class ProcessorLduInterfaceField
-{
-
-public:
-
-    //- Runtime type information
-    TypeName("ProcessorLduInterfaceField");
-
-
-    // Constructors
-
-        //- Construct given coupled patch
-        ProcessorLduInterfaceField()
-        {}
-
-
-    // Destructor
-
-        virtual ~ProcessorLduInterfaceField();
-
-
-    // Member Functions
-
-        // Access
-
-            //- Return processor number
-            virtual int myProcNo() const = 0;
-
-            //- Return neigbour processor number
-            virtual int neighbProcNo() const = 0;
-
-            //- Is the transform required
-            virtual bool doTransform() const = 0;
-
-            //- Return face transformation tensor
-            virtual const tensorField& forwardT() const = 0;
-
-            //- Return rank of component for transform
-            virtual int rank() const = 0;
-
-
-        //- Transform given patch component field
-        void transformCoupleField(Field<Type>& f) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
deleted file mode 100644
index caf49e595ca4d00e67d2330918d4575a77cc78af..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C
+++ /dev/null
@@ -1,391 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "lduMatrix.H"
-#include "IOstreams.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::LduMatrix(const lduMesh& mesh)
-:
-    lduMesh_(mesh),
-    diagPtr_(NULL),
-    upperPtr_(NULL),
-    lowerPtr_(NULL),
-    sourcePtr_(NULL),
-    interfaces_(0),
-    interfacesUpper_(0),
-    interfacesLower_(0)
-{}
-
-
-template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::LduMatrix(const LduMatrix& A)
-:
-    lduMesh_(A.lduMesh_),
-    diagPtr_(NULL),
-    upperPtr_(NULL),
-    lowerPtr_(NULL),
-    sourcePtr_(NULL),
-    interfaces_(0),
-    interfacesUpper_(0),
-    interfacesLower_(0)
-{
-    if (A.diagPtr_)
-    {
-        diagPtr_ = new Field<DType>(*(A.diagPtr_));
-    }
-
-    if (A.upperPtr_)
-    {
-        upperPtr_ = new Field<LUType>(*(A.upperPtr_));
-    }
-
-    if (A.lowerPtr_)
-    {
-        lowerPtr_ = new Field<LUType>(*(A.lowerPtr_));
-    }
-
-    if (A.sourcePtr_)
-    {
-        sourcePtr_ = new Field<Type>(*(A.sourcePtr_));
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::LduMatrix(LduMatrix& A, bool reUse)
-:
-    lduMesh_(A.lduMesh_),
-    diagPtr_(NULL),
-    upperPtr_(NULL),
-    lowerPtr_(NULL),
-    sourcePtr_(NULL),
-    interfaces_(0),
-    interfacesUpper_(0),
-    interfacesLower_(0)
-{
-    if (reUse)
-    {
-        if (A.diagPtr_)
-        {
-            diagPtr_ = A.diagPtr_;
-            A.diagPtr_ = NULL;
-        }
-
-        if (A.upperPtr_)
-        {
-            upperPtr_ = A.upperPtr_;
-            A.upperPtr_ = NULL;
-        }
-
-        if (A.lowerPtr_)
-        {
-            lowerPtr_ = A.lowerPtr_;
-            A.lowerPtr_ = NULL;
-        }
-
-        if (A.sourcePtr_)
-        {
-            sourcePtr_ = A.sourcePtr_;
-            A.sourcePtr_ = NULL;
-        }
-    }
-    else
-    {
-        if (A.diagPtr_)
-        {
-            diagPtr_ = new Field<DType>(*(A.diagPtr_));
-        }
-
-        if (A.upperPtr_)
-        {
-            upperPtr_ = new Field<LUType>(*(A.upperPtr_));
-        }
-
-        if (A.lowerPtr_)
-        {
-            lowerPtr_ = new Field<LUType>(*(A.lowerPtr_));
-        }
-
-        if (A.sourcePtr_)
-        {
-            sourcePtr_ = new Field<Type>(*(A.sourcePtr_));
-        }
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::LduMatrix
-(
-    const lduMesh& mesh,
-    Istream& is
-)
-:
-    lduMesh_(mesh),
-    diagPtr_(new Field<DType>(is)),
-    upperPtr_(new Field<LUType>(is)),
-    lowerPtr_(new Field<LUType>(is)),
-    sourcePtr_(new Field<Type>(is)),
-    interfaces_(0),
-    interfacesUpper_(0),
-    interfacesLower_(0)
-{}
-
-
-// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::~LduMatrix()
-{
-    if (diagPtr_)
-    {
-        delete diagPtr_;
-    }
-
-    if (upperPtr_)
-    {
-        delete upperPtr_;
-    }
-
-    if (lowerPtr_)
-    {
-        delete lowerPtr_;
-    }
-
-    if (sourcePtr_)
-    {
-        delete sourcePtr_;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag()
-{
-    if (!diagPtr_)
-    {
-        diagPtr_ = new Field<DType>(lduAddr().size(), pTraits<DType>::zero);
-    }
-
-    return *diagPtr_;
-}
-
-
-template<class Type, class DType, class LUType>
-Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper()
-{
-    if (!upperPtr_)
-    {
-        if (lowerPtr_)
-        {
-            upperPtr_ = new Field<LUType>(*lowerPtr_);
-        }
-        else
-        {
-            upperPtr_ = new Field<LUType>
-            (
-                lduAddr().lowerAddr().size(),
-                pTraits<LUType>::zero
-            );
-        }
-    }
-
-    return *upperPtr_;
-}
-
-
-template<class Type, class DType, class LUType>
-Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower()
-{
-    if (!lowerPtr_)
-    {
-        if (upperPtr_)
-        {
-            lowerPtr_ = new Field<LUType>(*upperPtr_);
-        }
-        else
-        {
-            lowerPtr_ = new Field<LUType>
-            (
-                lduAddr().lowerAddr().size(),
-                pTraits<LUType>::zero
-            );
-        }
-    }
-
-    return *lowerPtr_;
-}
-
-
-template<class Type, class DType, class LUType>
-Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source()
-{
-    if (!sourcePtr_)
-    {
-        sourcePtr_ = new Field<Type>(lduAddr().size(), pTraits<Type>::zero);
-    }
-
-    return *sourcePtr_;
-}
-
-
-template<class Type, class DType, class LUType>
-const Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag() const
-{
-    if (!diagPtr_)
-    {
-        FatalErrorIn
-        (
-            "const Field<DType>& LduMatrix<Type, DType, LUType>::diag() const"
-        )   << "diagPtr_ unallocated"
-            << abort(FatalError);
-    }
-
-    return *diagPtr_;
-}
-
-
-template<class Type, class DType, class LUType>
-const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper() const
-{
-    if (!lowerPtr_ && !upperPtr_)
-    {
-        FatalErrorIn
-        (
-            "const Field<LUType>& LduMatrix<Type, DType, LUType>::upper() const"
-        )   << "lowerPtr_ or upperPtr_ unallocated"
-            << abort(FatalError);
-    }
-
-    if (upperPtr_)
-    {
-        return *upperPtr_;
-    }
-    else
-    {
-        return *lowerPtr_;
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower() const
-{
-    if (!lowerPtr_ && !upperPtr_)
-    {
-        FatalErrorIn
-        (
-            "const Field<LUType>& LduMatrix<Type, DType, LUType>::lower() const"
-        )   << "lowerPtr_ or upperPtr_ unallocated"
-            << abort(FatalError);
-    }
-
-    if (lowerPtr_)
-    {
-        return *lowerPtr_;
-    }
-    else
-    {
-        return *upperPtr_;
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
-{
-    if (!sourcePtr_)
-    {
-        FatalErrorIn
-        (
-            "const Field<Type>& LduMatrix<Type, DType, LUType>::source() const"
-        )   << "sourcePtr_ unallocated"
-            << abort(FatalError);
-    }
-
-    return *sourcePtr_;
-}
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::Ostream& Foam::operator<<
-(
-    Ostream& os,
-    const LduMatrix<Type, DType, LUType>& ldum
-)
-{
-    if (ldum.diagPtr_)
-    {
-        os  << "Diagonal = "
-            << *ldum.diagPtr_
-            << endl << endl;
-    }
-
-    if (ldum.upperPtr_)
-    {
-        os  << "Upper triangle = "
-            << *ldum.upperPtr_
-            << endl << endl;
-    }
-
-    if (ldum.lowerPtr_)
-    {
-        os  << "Lower triangle = "
-            << *ldum.lowerPtr_
-            << endl << endl;
-    }
-
-    if (ldum.sourcePtr_)
-    {
-        os  << "Source = "
-            << *ldum.sourcePtr_
-            << endl << endl;
-    }
-
-    os.check("Ostream& operator<<(Ostream&, const LduMatrix&");
-
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "LduMatrixOperations.C"
-#include "LduMatrixATmul.C"
-#include "LduMatrixUpdateMatrixInterfaces.C"
-#include "LduMatrixTests.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
deleted file mode 100644
index e134c093641e6b4c67bb406ad7bf7762d956d810..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
+++ /dev/null
@@ -1,935 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::LduMatrix
-
-Description
-    LduMatrix 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.
-
-    Addressing arrays must be supplied for the upper and lower triangles.
-
-Note
-    It might be better if this class were organised as a hierachy starting
-    from an empty matrix, then deriving diagonal, symmetric and asymmetric
-    matrices.
-
-SourceFiles
-    LduMatrixATmul.C
-    LduMatrix.C
-    LduMatrixOperations.C
-    LduMatrixSolver.C
-    LduMatrixPreconditioner.C
-    LduMatrixTests.C
-    LduMatrixUpdateMatrixInterfaces.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef LduMatrix_H
-#define LduMatrix_H
-
-#include "lduMesh.H"
-#include "Field.H"
-#include "FieldField.H"
-#include "LduInterfaceFieldPtrsList.H"
-#include "typeInfo.H"
-#include "autoPtr.H"
-#include "runTimeSelectionTables.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of friend functions and operators
-
-template<class Type, class DType, class LUType>
-class LduMatrix;
-
-template<class Type, class DType, class LUType>
-Ostream& operator<<
-(
-    Ostream&,
-    const LduMatrix<Type, DType, LUType>&
-);
-
-
-/*---------------------------------------------------------------------------*\
-                           Class LduMatrix Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type, class DType, class LUType>
-class LduMatrix
-{
-    // private data
-
-        //- LDU mesh reference
-        const lduMesh& lduMesh_;
-
-        //- Diagonal coefficients
-        Field<DType> *diagPtr_;
-
-        //- Off-diagonal coefficients
-        Field<LUType> *upperPtr_, *lowerPtr_;
-
-        //- Source
-        Field<Type> *sourcePtr_;
-
-        //- Field interfaces (processor patches etc.)
-        LduInterfaceFieldPtrsList<Type> interfaces_;
-
-        //- Off-diagonal coefficients for interfaces
-        FieldField<Field, LUType> interfacesUpper_, interfacesLower_;
-
-
-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 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_;
-            }
-
-            //- Check, store and return singularity
-            bool singular(const Type& wApA);
-
-            //- Is the matrix singular?
-            bool singular() const;
-
-            //- Check, store and return convergence
-            bool converged
-            (
-                const Type& tolerance,
-                const Type& relTolerance
-            );
-
-            //- Has the solver converged?
-            bool converged() const
-            {
-                return converged_;
-            }
-
-            //- Print summary of solver performance to the given stream
-            void print(Ostream& os) const;
-    };
-
-
-    //- Abstract base-class for LduMatrix solvers
-    class solver
-    {
-    protected:
-
-        // Protected data
-
-            word fieldName_;
-            const LduMatrix<Type, DType, LUType>& matrix_;
-
-            //- dictionary of controls
-            dictionary controlDict_;
-
-            //- Maximum number of iterations in the solver
-            label maxIter_;
-
-            //- Final convergence tolerance
-            Type tolerance_;
-
-            //- Convergence tolerance relative to the initial
-            Type relTol_;
-
-
-        // Protected Member Functions
-
-            //- Read a control parameter from controlDict
-            template<class T>
-            inline void readControl
-            (
-                const dictionary& controlDict,
-                T& control,
-                const word& controlName
-            );
-
-
-            //- Read the control parameters from the controlDict_
-            virtual void readControls();
-
-
-    public:
-
-        //- Runtime type information
-        virtual const word& type() const = 0;
-
-
-        // Declare run-time constructor selection tables
-
-            declareRunTimeSelectionTable
-            (
-                autoPtr,
-                solver,
-                symMatrix,
-                (
-                    const word& fieldName,
-                    const LduMatrix<Type, DType, LUType>& matrix,
-                    const dictionary& solverDict
-                ),
-                (
-                    fieldName,
-                    matrix,
-                    solverDict
-                )
-            );
-
-            declareRunTimeSelectionTable
-            (
-                autoPtr,
-                solver,
-                asymMatrix,
-                (
-                    const word& fieldName,
-                    const LduMatrix<Type, DType, LUType>& matrix,
-                    const dictionary& solverDict
-                ),
-                (
-                    fieldName,
-                    matrix,
-                    solverDict
-                )
-            );
-
-
-        // Constructors
-
-            solver
-            (
-                const word& fieldName,
-                const LduMatrix<Type, DType, LUType>& matrix,
-                const dictionary& solverDict
-            );
-
-
-        // Selectors
-
-            //- Return a new solver
-            static autoPtr<solver> New
-            (
-                const word& fieldName,
-                const LduMatrix<Type, DType, LUType>& matrix,
-                const dictionary& solverDict
-            );
-
-
-        // Destructor
-
-            virtual ~solver()
-            {}
-
-
-        // Member functions
-
-            // Access
-
-                const word& fieldName() const
-                {
-                    return fieldName_;
-                }
-
-                const LduMatrix<Type, DType, LUType>& matrix() const
-                {
-                    return matrix_;
-                }
-
-
-            //- Read and reset the solver parameters from the given dictionary
-            virtual void read(const dictionary& solverDict);
-
-            virtual solverPerformance solve
-            (
-                Field<Type>& psi
-            ) const = 0;
-
-            //- Return the matrix norm used to normalise the residual for the
-            //  stopping criterion
-            Type normFactor
-            (
-                const Field<Type>& psi,
-                const Field<Type>& Apsi,
-                Field<Type>& tmpField
-            ) const;
-    };
-
-
-    //- Abstract base-class for LduMatrix smoothers
-    class smoother
-    {
-    protected:
-
-        // Protected data
-
-            word fieldName_;
-            const LduMatrix<Type, DType, LUType>& matrix_;
-
-
-    public:
-
-        //- Runtime type information
-        virtual const word& type() const = 0;
-
-
-        // Declare run-time constructor selection tables
-
-            declareRunTimeSelectionTable
-            (
-                autoPtr,
-                smoother,
-                symMatrix,
-                (
-                    const word& fieldName,
-                    const LduMatrix<Type, DType, LUType>& matrix
-                ),
-                (
-                    fieldName,
-                    matrix
-                )
-            );
-
-            declareRunTimeSelectionTable
-            (
-                autoPtr,
-                smoother,
-                asymMatrix,
-                (
-                    const word& fieldName,
-                    const LduMatrix<Type, DType, LUType>& matrix
-                ),
-                (
-                    fieldName,
-                    matrix
-                )
-            );
-
-
-        // Constructors
-
-            smoother
-            (
-                const word& fieldName,
-                const LduMatrix<Type, DType, LUType>& matrix
-            );
-
-
-        // Selectors
-
-            //- Return a new smoother
-            static autoPtr<smoother> New
-            (
-                const word& fieldName,
-                const LduMatrix<Type, DType, LUType>& matrix,
-                const dictionary& smootherDict
-            );
-
-
-        // Destructor
-
-            virtual ~smoother()
-            {}
-
-
-        // Member functions
-
-            // Access
-
-                const word& fieldName() const
-                {
-                    return fieldName_;
-                }
-
-                const LduMatrix<Type, DType, LUType>& matrix() const
-                {
-                    return matrix_;
-                }
-
-
-            //- Smooth the solution for a given number of sweeps
-            virtual void smooth
-            (
-                Field<Type>& psi,
-                const label nSweeps
-            ) const = 0;
-    };
-
-
-    //- Abstract base-class for LduMatrix preconditioners
-    class preconditioner
-    {
-    protected:
-
-        // Protected data
-
-            //- Reference to the base-solver this preconditioner is used with
-            const solver& solver_;
-
-
-    public:
-
-        //- Runtime type information
-        virtual const word& type() const = 0;
-
-
-        // Declare run-time constructor selection tables
-
-            declareRunTimeSelectionTable
-            (
-                autoPtr,
-                preconditioner,
-                symMatrix,
-                (
-                    const solver& sol,
-                    const dictionary& preconditionerDict
-                ),
-                (sol, preconditionerDict)
-            );
-
-            declareRunTimeSelectionTable
-            (
-                autoPtr,
-                preconditioner,
-                asymMatrix,
-                (
-                    const solver& sol,
-                    const dictionary& preconditionerDict
-                ),
-                (sol, preconditionerDict)
-            );
-
-
-        // Constructors
-
-            preconditioner
-            (
-                const solver& sol
-            )
-            :
-                solver_(sol)
-            {}
-
-
-        // Selectors
-
-            //- Return a new preconditioner
-            static autoPtr<preconditioner> New
-            (
-                const solver& sol,
-                const dictionary& preconditionerDict
-            );
-
-
-        // Destructor
-
-            virtual ~preconditioner()
-            {}
-
-
-        // Member functions
-
-            //- Read and reset the preconditioner parameters
-            //  from the given dictionary
-            virtual void read(const dictionary& preconditionerDict)
-            {}
-
-            //- Return wA the preconditioned form of residual rA
-            virtual void precondition
-            (
-                Field<Type>& wA,
-                const Field<Type>& rA
-            ) const = 0;
-
-            //- Return wT the transpose-matrix preconditioned form of
-            //  residual rT.
-            //  This is only required for preconditioning asymmetric matrices.
-            virtual void preconditionT
-            (
-                Field<Type>& wT,
-                const Field<Type>& rT
-            ) const
-            {
-                notImplemented
-                (
-                    type() +"::preconditionT"
-                    "(Field<Type>& wT, const Field<Type>& rT)"
-                );
-            }
-    };
-
-
-    // Static data
-
-        // Declare name of the class and its debug switch
-        ClassName("LduMatrix");
-
-        //- Large Type for the use in solvers
-        static const scalar great_;
-
-        //- Small Type for the use in solvers
-        static const scalar small_;
-
-        //- Very small Type for the use in solvers
-        static const scalar vsmall_;
-
-
-    // Constructors
-
-        //- Construct given an LDU addressed mesh.
-        //  The coefficients are initially empty for subsequent setting.
-        LduMatrix(const lduMesh&);
-
-        //- Construct as copy
-        LduMatrix(const LduMatrix<Type, DType, LUType>&);
-
-        //- Construct as copy or re-use as specified.
-        LduMatrix(LduMatrix<Type, DType, LUType>&, bool reUse);
-
-        //- Construct given an LDU addressed mesh and an Istream
-        //  from which the coefficients are read
-        LduMatrix(const lduMesh&, Istream&);
-
-
-    // Destructor
-
-        ~LduMatrix();
-
-
-    // Member functions
-
-        // Access to addressing
-
-            //- Return the LDU mesh from which the addressing is obtained
-            const lduMesh& mesh() const
-            {
-                return lduMesh_;
-            }
-
-            //- Return the LDU addressing
-            const lduAddressing& lduAddr() const
-            {
-                return lduMesh_.lduAddr();
-            }
-
-            //- Return the patch evaluation schedule
-            const lduSchedule& patchSchedule() const
-            {
-                return lduAddr().patchSchedule();
-            }
-
-            //- Return interfaces
-            const LduInterfaceFieldPtrsList<Type>& interfaces() const
-            {
-                return interfaces_;
-            }
-
-
-        // Access to coefficients
-
-            Field<DType>& diag();
-            Field<LUType>& upper();
-            Field<LUType>& lower();
-            Field<Type>& source();
-
-            FieldField<Field, LUType>& interfacesUpper()
-            {
-                return interfacesUpper_;
-            }
-
-            FieldField<Field, LUType>& interfacesLower()
-            {
-                return interfacesLower_;
-            }
-
-
-            const Field<DType>& diag() const;
-            const Field<LUType>& upper() const;
-            const Field<LUType>& lower() const;
-            const Field<Type>& source() const;
-
-            const FieldField<Field, LUType>& interfacesUpper() const
-            {
-                return interfacesUpper_;
-            }
-
-            const FieldField<Field, LUType>& interfacesLower() const
-            {
-                return interfacesLower_;
-            }
-
-
-            bool hasDiag() const
-            {
-                return (diagPtr_);
-            }
-
-            bool hasUpper() const
-            {
-                return (upperPtr_);
-            }
-
-            bool hasLower() const
-            {
-                return (lowerPtr_);
-            }
-
-            bool hasSource() const
-            {
-                return (sourcePtr_);
-            }
-
-            bool diagonal() const
-            {
-                return (diagPtr_ && !lowerPtr_ && !upperPtr_);
-            }
-
-            bool symmetric() const
-            {
-                return (diagPtr_ && (!lowerPtr_ && upperPtr_));
-            }
-
-            bool asymmetric() const
-            {
-                return (diagPtr_ && lowerPtr_ && upperPtr_);
-            }
-
-
-        // operations
-
-            void sumDiag();
-            void negSumDiag();
-
-            void sumMagOffDiag(Field<LUType>& sumOff) const;
-
-            //- Matrix multiplication
-            void Amul(Field<Type>&, const tmp<Field<Type> >&) const;
-    
-            //- Matrix transpose multiplication
-            void Tmul(Field<Type>&, const tmp<Field<Type> >&) const;
-
-
-            //- Sum the coefficients on each row of the matrix
-            void sumA(Field<Type>&) const;
-
-
-            void residual(Field<Type>& rA, const Field<Type>& psi) const;
-
-            tmp<Field<Type> > residual(const Field<Type>& psi) const;
-
-
-            //- Initialise the update of interfaced interfaces
-            //  for matrix operations
-            void initMatrixInterfaces
-            (
-                const FieldField<Field, LUType>& interfaceCoeffs,
-                const Field<Type>& psiif,
-                Field<Type>& result
-            ) const;
-    
-            //- Update interfaced interfaces for matrix operations
-            void updateMatrixInterfaces
-            (
-                const FieldField<Field, LUType>& interfaceCoeffs,
-                const Field<Type>& psiif,
-                Field<Type>& result
-            ) const;
-    
-
-            tmp<Field<Type> > H(const Field<Type>&) const;
-            tmp<Field<Type> > H(const tmp<Field<Type> >&) const;
-
-            tmp<Field<Type> > faceH(const Field<Type>&) const;
-            tmp<Field<Type> > faceH(const tmp<Field<Type> >&) const;
-
-
-    // Member operators
-
-        void operator=(const LduMatrix<Type, DType, LUType>&);
-
-        void negate();
-
-        void operator+=(const LduMatrix<Type, DType, LUType>&);
-        void operator-=(const LduMatrix<Type, DType, LUType>&);
-
-        void operator*=(const scalarField&);
-        void operator*=(scalar);
-
-
-    // Ostream operator
-
-        friend Ostream& operator<< <Type, DType, LUType>
-        (
-            Ostream&,
-            const LduMatrix<Type, DType, LUType>&
-        );
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#define makeLduMatrix(Type, DType, LUType)                                    \
-                                                                              \
-typedef Foam::LduMatrix<Type, DType, LUType>                                  \
-    ldu##Type##DType##LUType##Matrix;                                         \
-                                                                              \
-defineNamedTemplateTypeNameAndDebug(ldu##Type##DType##LUType##Matrix, 0);     \
-                                                                              \
-                                                                              \
-template<>                                                                    \
-const scalar ldu##Type##DType##LUType##Matrix::great_(1e15);                  \
-                                                                              \
-template<>                                                                    \
-const scalar ldu##Type##DType##LUType##Matrix::small_(1e-15);                 \
-                                                                              \
-template<>                                                                    \
-const scalar ldu##Type##DType##LUType##Matrix::vsmall_(VSMALL);               \
-                                                                              \
-                                                                              \
-typedef LduMatrix<Type, DType, LUType>::smoother                              \
-    ldu##Type##DType##LUType##Smoother;                                       \
-                                                                              \
-defineTemplateRunTimeSelectionTable                                           \
-(                                                                             \
-    ldu##Type##DType##LUType##Smoother,                                       \
-    symMatrix                                                                 \
-);                                                                            \
-                                                                              \
-defineTemplateRunTimeSelectionTable                                           \
-(                                                                             \
-    ldu##Type##DType##LUType##Smoother,                                       \
-    asymMatrix                                                                \
-);                                                                            \
-                                                                              \
-                                                                              \
-typedef LduMatrix<Type, DType, LUType>::preconditioner                        \
-    ldu##Type##DType##LUType##Preconditioner;                                 \
-                                                                              \
-defineTemplateRunTimeSelectionTable                                           \
-(                                                                             \
-    ldu##Type##DType##LUType##Preconditioner,                                 \
-    symMatrix                                                                 \
-);                                                                            \
-                                                                              \
-defineTemplateRunTimeSelectionTable                                           \
-(                                                                             \
-    ldu##Type##DType##LUType##Preconditioner,                                 \
-    asymMatrix                                                                \
-);                                                                            \
-                                                                              \
-                                                                              \
-typedef LduMatrix<Type, DType, LUType>::solver                                \
-    ldu##Type##DType##LUType##Solver;                                         \
-                                                                              \
-defineTemplateRunTimeSelectionTable                                           \
-(                                                                             \
-    ldu##Type##DType##LUType##Solver,                                         \
-    symMatrix                                                                 \
-);                                                                            \
-                                                                              \
-defineTemplateRunTimeSelectionTable                                           \
-(                                                                             \
-    ldu##Type##DType##LUType##Solver,                                         \
-    asymMatrix                                                                \
-);
-
-
-#define makeLduPreconditioner(Precon, Type, DType, LUType)                    \
-                                                                              \
-typedef Precon<Type, DType, LUType>                                           \
-    Precon##Type##DType##LUType##Preconditioner;                              \
-defineNamedTemplateTypeNameAndDebug                                           \
-(                                                                             \
-    Precon##Type##DType##LUType##Preconditioner,                              \
-    0                                                                         \
-);
-
-#define makeLduSymPreconditioner(Precon, Type, DType, LUType)                 \
-                                                                              \
-LduMatrix<Type, DType, LUType>::preconditioner::                              \
-addsymMatrixConstructorToTable<Precon##Type##DType##LUType##Preconditioner>   \
-add##Precon##Type##DType##LUType##PreconditionerSymMatrixConstructorToTable_;
-
-#define makeLduAsymPreconditioner(Precon, Type, DType, LUType)                \
-                                                                              \
-LduMatrix<Type, DType, LUType>::preconditioner::                              \
-addasymMatrixConstructorToTable<Precon##Type##DType##LUType##Preconditioner>  \
-add##Precon##Type##DType##LUType##PreconditionerAsymMatrixConstructorToTable_;
-
-
-#define makeLduSmoother(Smoother, Type, DType, LUType)                        \
-                                                                              \
-typedef Smoother<Type, DType, LUType>                                         \
-    Smoother##Type##DType##LUType##Smoother;                                  \
-                                                                              \
-defineNamedTemplateTypeNameAndDebug                                           \
-(                                                                             \
-    Smoother##Type##DType##LUType##Smoother,                                  \
-    0                                                                         \
-);
-
-#define makeLduSymSmoother(Smoother, Type, DType, LUType)                     \
-                                                                              \
-LduMatrix<Type, DType, LUType>::smoother::                                    \
-    addsymMatrixConstructorToTable<Smoother##Type##DType##LUType##Smoother>   \
-    add##Smoother##Type##DType##LUType##SymMatrixConstructorToTable_;
-
-#define makeLduAsymSmoother(Smoother, Type, DType, LUType)                    \
-                                                                              \
-LduMatrix<Type, DType, LUType>::smoother::                                    \
-    addasymMatrixConstructorToTable<Smoother##Type##DType##LUType##Smoother>  \
-    add##Smoother##Type##DType##LUType##AsymMatrixConstructorToTable_;
-
-
-#define makeLduSolver(Solver, Type, DType, LUType)                            \
-                                                                              \
-typedef Solver<Type, DType, LUType>                                           \
-    Solver##Type##DType##LUType##Solver;                                      \
-                                                                              \
-defineNamedTemplateTypeNameAndDebug                                           \
-(                                                                             \
-    Solver##Type##DType##LUType##Solver,                                      \
-    0                                                                         \
-);
-
-#define makeLduSymSolver(Solver, Type, DType, LUType)                         \
-                                                                              \
-LduMatrix<Type, DType, LUType>::solver::                                      \
-    addsymMatrixConstructorToTable<Solver##Type##DType##LUType##Solver>       \
-    add##Solver##Type##DType##LUType##SymMatrixConstructorToTable_;
-
-#define makeLduAsymSolver(Solver, Type, DType, LUType)                        \
-                                                                              \
-LduMatrix<Type, DType, LUType>::solver::                                      \
-    addasymMatrixConstructorToTable<Solver##Type##DType##LUType##Solver>      \
-    add##Solver##Type##DType##LUType##AsymMatrixConstructorToTable_;
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "LduMatrixI.H"
-#   include "LduMatrix.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixATmul.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixATmul.C
deleted file mode 100644
index 27832d64a14d96825be74da5a49ad1f7d633f9d0..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixATmul.C
+++ /dev/null
@@ -1,291 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-template<class Type, class LUType>
-class Amultiplier
-:
-    public LduInterfaceField<Type>::Amultiplier
-{
-    const Field<LUType>& A_;
-
-public:
-
-    Amultiplier(const Field<LUType>& A)
-    :
-        A_(A)
-    {}
-
-    virtual void addAmul(Field<Type>& Apsi, const Field<Type>& psi) const
-    {
-        Apsi += A_*psi;
-    }
-};
-
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::Amul
-(
-    Field<Type>& Apsi,
-    const tmp<Field<Type> >& tpsi
-) const
-{
-    Type* __restrict__ ApsiPtr = Apsi.begin();
-
-    const Field<Type>& psi = tpsi();
-    const Type* const __restrict__ psiPtr = psi.begin();
-
-    const DType* const __restrict__ diagPtr = diag().begin();
-
-    const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
-    const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
-
-    const LUType* const __restrict__ upperPtr = upper().begin();
-    const LUType* const __restrict__ lowerPtr = lower().begin();
-
-    // Initialise the update of interfaced interfaces
-    initMatrixInterfaces
-    (
-        interfacesUpper_,
-        psi,
-        Apsi
-    );
-
-    register const label nCells = diag().size();
-    for (register label cell=0; cell<nCells; cell++)
-    {
-        ApsiPtr[cell] = dot(diagPtr[cell], psiPtr[cell]);
-    }
-
-
-    register const label nFaces = upper().size();
-    for (register label face=0; face<nFaces; face++)
-    {
-        ApsiPtr[uPtr[face]] += dot(lowerPtr[face], psiPtr[lPtr[face]]);
-        ApsiPtr[lPtr[face]] += dot(upperPtr[face], psiPtr[uPtr[face]]);
-    }
-
-    // Update interface interfaces
-    updateMatrixInterfaces
-    (
-        interfacesUpper_,
-        psi,
-        Apsi
-    );
-
-    tpsi.clear();
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::Tmul
-(
-    Field<Type>& Tpsi,
-    const tmp<Field<Type> >& tpsi
-) const
-{
-    Type* __restrict__ TpsiPtr = Tpsi.begin();
-
-    const Field<Type>& psi = tpsi();
-    const Type* const __restrict__ psiPtr = psi.begin();
-
-    const DType* const __restrict__ diagPtr = diag().begin();
-
-    const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
-    const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
-
-    const LUType* const __restrict__ lowerPtr = lower().begin();
-    const LUType* const __restrict__ upperPtr = upper().begin();
-
-    // Initialise the update of interfaced interfaces
-    initMatrixInterfaces
-    (
-        interfacesLower_,
-        psi,
-        Tpsi
-    );
-
-    register const label nCells = diag().size();
-    for (register label cell=0; cell<nCells; cell++)
-    {
-        TpsiPtr[cell] = dot(diagPtr[cell], psiPtr[cell]);
-    }
-
-    register const label nFaces = upper().size();
-    for (register label face=0; face<nFaces; face++)
-    {
-        TpsiPtr[uPtr[face]] += dot(upperPtr[face], psiPtr[lPtr[face]]);
-        TpsiPtr[lPtr[face]] += dot(lowerPtr[face], psiPtr[uPtr[face]]);
-    }
-
-    // Update interface interfaces
-    updateMatrixInterfaces
-    (
-        interfacesLower_,
-        psi,
-        Tpsi
-    );
-
-    tpsi.clear();
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::sumA
-(
-    Field<Type>& sumA
-) const
-{
-    Type* __restrict__ sumAPtr = sumA.begin();
-
-    const DType* __restrict__ diagPtr = diag().begin();
-
-    const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
-    const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
-
-    const LUType* __restrict__ lowerPtr = lower().begin();
-    const LUType* __restrict__ upperPtr = upper().begin();
-
-    register const label nCells = diag().size();
-    register const label nFaces = upper().size();
-
-    for (register label cell=0; cell<nCells; cell++)
-    {
-        sumAPtr[cell] = dot(diagPtr[cell], pTraits<Type>::one);
-    }
-
-    for (register label face=0; face<nFaces; face++)
-    {
-        sumAPtr[uPtr[face]] += dot(lowerPtr[face], pTraits<Type>::one);
-        sumAPtr[lPtr[face]] += dot(upperPtr[face], pTraits<Type>::one);
-    }
-
-    // Add the interface internal coefficients to diagonal
-    // and the interface boundary coefficients to the sum-off-diagonal
-    forAll(interfaces_, patchI)
-    {
-        if (interfaces_.set(patchI))
-        {
-            const unallocLabelList& pa = lduAddr().patchAddr(patchI);
-            const Field<LUType>& pCoeffs = interfacesUpper_[patchI];
-
-            forAll(pa, face)
-            {
-                sumAPtr[pa[face]] -= dot(pCoeffs[face], pTraits<Type>::one);
-            }
-        }
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::residual
-(
-    Field<Type>& rA,
-    const Field<Type>& psi
-) const
-{
-    Type* __restrict__ rAPtr = rA.begin();
-
-    const Type* const __restrict__ psiPtr = psi.begin();
-    const DType* const __restrict__ diagPtr = diag().begin();
-    const Type* const __restrict__ sourcePtr = source().begin();
-
-    const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
-    const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
-
-    const LUType* const __restrict__ upperPtr = upper().begin();
-    const LUType* const __restrict__ lowerPtr = lower().begin();
-
-    // Parallel boundary initialisation.
-    // Note: there is a change of sign in the coupled
-    // interface update to add the contibution to the r.h.s.
-
-    FieldField<Field, LUType> mBouCoeffs(interfacesUpper_.size());
-
-    forAll(mBouCoeffs, patchi)
-    {
-        if (interfaces_.set(patchi))
-        {
-            mBouCoeffs.set(patchi, -interfacesUpper_[patchi]);
-        }
-    }
-
-    // Initialise the update of interfaced interfaces
-    initMatrixInterfaces
-    (
-        mBouCoeffs,
-        psi,
-        rA
-    );
-
-    register const label nCells = diag().size();
-    for (register label cell=0; cell<nCells; cell++)
-    {
-        rAPtr[cell] = sourcePtr[cell] - dot(diagPtr[cell], psiPtr[cell]);
-    }
-
-
-    register const label nFaces = upper().size();
-    for (register label face=0; face<nFaces; face++)
-    {
-        rAPtr[uPtr[face]] -= dot(lowerPtr[face], psiPtr[lPtr[face]]);
-        rAPtr[lPtr[face]] -= dot(upperPtr[face], psiPtr[uPtr[face]]);
-    }
-
-    // Update interface interfaces
-    updateMatrixInterfaces
-    (
-        mBouCoeffs,
-        psi,
-        rA
-    );
-}
-
-
-template<class Type, class DType, class LUType>
-Foam::tmp<Foam::Field<Type> > Foam::LduMatrix<Type, DType, LUType>::residual
-(
-    const Field<Type>& psi
-) const
-{
-    tmp<Field<Type> > trA(new Field<Type>(psi.size()));
-    residual(trA(), psi);
-    return trA;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixI.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixI.H
deleted file mode 100644
index af010040025baa6fa2910a4b166d846c0156b38a..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixI.H
+++ /dev/null
@@ -1,45 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-template<class T>
-inline void Foam::LduMatrix<Type, DType, LUType>::solver::readControl
-(
-    const dictionary& controlDict,
-    T& control,
-    const word& controlName
-)
-{
-    if (controlDict.found(controlName))
-    {
-        controlDict.lookup(controlName) >> control;
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C
deleted file mode 100644
index 75ab0ecd097805e1794886b9f55569efe13a438c..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C
+++ /dev/null
@@ -1,481 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "lduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::sumDiag()
-{
-    const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
-    const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
-    Field<DType>& Diag = diag();
-
-    const unallocLabelList& l = lduAddr().lowerAddr();
-    const unallocLabelList& u = lduAddr().upperAddr();
-
-    for (register label face=0; face<l.size(); face++)
-    {
-        Diag[l[face]] += Lower[face];
-        Diag[u[face]] += Upper[face];
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::negSumDiag()
-{
-    const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
-    const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
-    Field<DType>& Diag = diag();
-
-    const unallocLabelList& l = lduAddr().lowerAddr();
-    const unallocLabelList& u = lduAddr().upperAddr();
-
-    for (register label face=0; face<l.size(); face++)
-    {
-        Diag[l[face]] -= Lower[face];
-        Diag[u[face]] -= Upper[face];
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::sumMagOffDiag
-(
-    Field<LUType>& sumOff
-) const
-{
-    const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
-    const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
-
-    const unallocLabelList& l = lduAddr().lowerAddr();
-    const unallocLabelList& u = lduAddr().upperAddr();
-
-    for (register label face = 0; face < l.size(); face++)
-    {
-        sumOff[u[face]] += cmptMag(Lower[face]);
-        sumOff[l[face]] += cmptMag(Upper[face]);
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-Foam::tmp<Foam::Field<Type> >
-Foam::LduMatrix<Type, DType, LUType>::H(const Field<Type>& psi) const
-{
-    tmp<Field<Type> > tHpsi
-    (
-        new Field<Type>(lduAddr().size(), pTraits<Type>::zero)
-    );
-
-    if (lowerPtr_ || upperPtr_)
-    {
-        Field<Type> & Hpsi = tHpsi();
-
-        Type* __restrict__ HpsiPtr = Hpsi.begin();
-
-        const Type* __restrict__ psiPtr = psi.begin();
-
-        const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
-        const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
-
-        const LUType* __restrict__ lowerPtr = lower().begin();
-        const LUType* __restrict__ upperPtr = upper().begin();
-
-        register const label nFaces = upper().size();
-
-        for (register label face=0; face<nFaces; face++)
-        {
-            HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
-            HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
-        }
-    }
-
-    return tHpsi;
-}
-
-template<class Type, class DType, class LUType>
-Foam::tmp<Foam::Field<Type> >
-Foam::LduMatrix<Type, DType, LUType>::H(const tmp<Field<Type> >& tpsi) const
-{
-    tmp<Field<Type> > tHpsi(H(tpsi()));
-    tpsi.clear();
-    return tHpsi;
-}
-
-
-template<class Type, class DType, class LUType>
-Foam::tmp<Foam::Field<Type> >
-Foam::LduMatrix<Type, DType, LUType>::faceH(const Field<Type>& psi) const
-{
-    const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
-    const Field<LUType>& 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++)
-    {
-        faceHpsi[face] = Upper[face]*psi[u[face]] - Lower[face]*psi[l[face]];
-    }
-
-    return tfaceHpsi;
-}
-
-
-template<class Type, class DType, class LUType>
-Foam::tmp<Foam::Field<Type> >
-Foam::LduMatrix<Type, DType, LUType>::faceH(const tmp<Field<Type> >& tpsi) const
-{
-    tmp<Field<Type> > tfaceHpsi(faceH(tpsi()));
-    tpsi.clear();
-    return tfaceHpsi;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::operator=(const LduMatrix& A)
-{
-    if (this == &A)
-    {
-        FatalErrorIn
-        (
-            "LduMatrix<Type, DType, LUType>::operator=(const LduMatrix&)"
-        )   << "attempted assignment to self"
-            << abort(FatalError);
-    }
-
-    if (A.diagPtr_)
-    {
-        diag() = A.diag();
-    }
-
-    if (A.upperPtr_)
-    {
-        upper() = A.upper();
-    }
-    else if (upperPtr_)
-    {
-        delete upperPtr_;
-        upperPtr_ = NULL;
-    }
-
-    if (A.lowerPtr_)
-    {
-        lower() = A.lower();
-    }
-    else if (lowerPtr_)
-    {
-        delete lowerPtr_;
-        lowerPtr_ = NULL;
-    }
-
-    if (A.sourcePtr_)
-    {
-        source() = A.source();
-    }
-
-    interfacesUpper_ = A.interfacesUpper_;
-    interfacesLower_ = A.interfacesLower_;
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::negate()
-{
-    if (diagPtr_)
-    {
-        diagPtr_->negate();
-    }
-
-    if (upperPtr_)
-    {
-        upperPtr_->negate();
-    }
-
-    if (lowerPtr_)
-    {
-        lowerPtr_->negate();
-    }
-
-    if (sourcePtr_)
-    {
-        sourcePtr_->negate();
-    }
-
-    negate(interfacesUpper_);
-    negate(interfacesLower_);
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::operator+=(const LduMatrix& A)
-{
-    if (A.diagPtr_)
-    {
-        diag() += A.diag();
-    }
-
-    if (A.sourcePtr_)
-    {
-        source() += A.source();
-    }
-
-    if (symmetric() && A.symmetric())
-    {
-        upper() += A.upper();
-    }
-    else if (symmetric() && A.asymmetric())
-    {
-        if (upperPtr_)
-        {
-            lower();
-        }
-        else
-        {
-            upper();
-        }
-
-        upper() += A.upper();
-        lower() += A.lower();
-    }
-    else if (asymmetric() && A.symmetric())
-    {
-        if (A.upperPtr_)
-        {
-            lower() += A.upper();
-            upper() += A.upper();
-        }
-        else
-        {
-            lower() += A.lower();
-            upper() += A.lower();
-        }
-
-    }
-    else if (asymmetric() && A.asymmetric())
-    {
-        lower() += A.lower();
-        upper() += A.upper();
-    }
-    else if (diagonal())
-    {
-        if (A.upperPtr_)
-        {
-            upper() = A.upper();
-        }
-
-        if (A.lowerPtr_)
-        {
-            lower() = A.lower();
-        }
-    }
-    else if (A.diagonal())
-    {
-    }
-    else
-    {
-        FatalErrorIn
-        (
-            "LduMatrix<Type, DType, LUType>::operator+=(const LduMatrix& A)"
-        )   << "Unknown matrix type combination"
-            << abort(FatalError);
-    }
-
-    interfacesUpper_ += A.interfacesUpper_;
-    interfacesLower_ += A.interfacesLower_;
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::operator-=(const LduMatrix& A)
-{
-    if (A.diagPtr_)
-    {
-        diag() -= A.diag();
-    }
-
-    if (A.sourcePtr_)
-    {
-        source() -= A.source();
-    }
-
-    if (symmetric() && A.symmetric())
-    {
-        upper() -= A.upper();
-    }
-    else if (symmetric() && A.asymmetric())
-    {
-        if (upperPtr_)
-        {
-            lower();
-        }
-        else
-        {
-            upper();
-        }
-
-        upper() -= A.upper();
-        lower() -= A.lower();
-    }
-    else if (asymmetric() && A.symmetric())
-    {
-        if (A.upperPtr_)
-        {
-            lower() -= A.upper();
-            upper() -= A.upper();
-        }
-        else
-        {
-            lower() -= A.lower();
-            upper() -= A.lower();
-        }
-
-    }
-    else if (asymmetric() && A.asymmetric())
-    {
-        lower() -= A.lower();
-        upper() -= A.upper();
-    }
-    else if (diagonal())
-    {
-        if (A.upperPtr_)
-        {
-            upper() = -A.upper();
-        }
-
-        if (A.lowerPtr_)
-        {
-            lower() = -A.lower();
-        }
-    }
-    else if (A.diagonal())
-    {
-    }
-    else
-    {
-        FatalErrorIn
-        (
-            "LduMatrix<Type, DType, LUType>::operator-=(const LduMatrix& A)"
-        )   << "Unknown matrix type combination"
-            << abort(FatalError);
-    }
-
-    interfacesUpper_ -= A.interfacesUpper_;
-    interfacesLower_ -= A.interfacesLower_;
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::operator*=
-(
-    const scalarField& sf
-)
-{
-    if (diagPtr_)
-    {
-        *diagPtr_ *= sf;
-    }
-
-    if (sourcePtr_)
-    {
-        *sourcePtr_ *= sf;
-    }
-
-    if (upperPtr_)
-    {
-        Field<LUType>& upper = *upperPtr_;
-
-        const unallocLabelList& l = lduAddr().lowerAddr();
-
-        for (register label face=0; face<upper.size(); face++)
-        {
-            upper[face] *= sf[l[face]];
-        }
-    }
-
-    if (lowerPtr_)
-    {
-        Field<LUType>& lower = *lowerPtr_;
-
-        const unallocLabelList& u = lduAddr().upperAddr();
-
-        for (register label face=0; face<lower.size(); face++)
-        {
-            lower[face] *= sf[u[face]];
-        }
-    }
-
-    FatalErrorIn
-    (
-        "LduMatrix<Type, DType, LUType>::operator*=(const scalarField& sf)"
-    )   << "Scaling a matrix by scalarField is not currently supported\n"
-           "because scaling interfacesUpper_ and interfacesLower_ "
-           "require special transfers"
-        << abort(FatalError);
-
-    //interfacesUpper_ *= ;
-    //interfacesLower_ *= sf;
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::operator*=(scalar s)
-{
-    if (diagPtr_)
-    {
-        *diagPtr_ *= s;
-    }
-
-    if (sourcePtr_)
-    {
-        *sourcePtr_ *= s;
-    }
-
-    if (upperPtr_)
-    {
-        *upperPtr_ *= s;
-    }
-
-    if (lowerPtr_)
-    {
-        *lowerPtr_ *= s;
-    }
-
-    interfacesUpper_ *= s;
-    interfacesLower_ *= s;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixPreconditioner.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixPreconditioner.C
deleted file mode 100644
index 50e4315b151ab1c676b7926824572a322f183396..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixPreconditioner.C
+++ /dev/null
@@ -1,116 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::autoPtr<typename Foam::LduMatrix<Type, DType, LUType>::preconditioner>
-Foam::LduMatrix<Type, DType, LUType>::preconditioner::New
-(
-    const solver& sol,
-    const dictionary& preconditionerDict
-)
-{
-    word preconditionerName = preconditionerDict.lookup("preconditioner");
-
-    if (sol.matrix().symmetric())
-    {
-        typename symMatrixConstructorTable::iterator constructorIter =
-            symMatrixConstructorTablePtr_->find(preconditionerName);
-
-        if (constructorIter == symMatrixConstructorTablePtr_->end())
-        {
-            FatalIOErrorIn
-            (
-                "LduMatrix<Type, DType, LUType>::preconditioner::New"
-                "(const solver&, Istream&)",
-                preconditionerDict
-            )   << "Unknown symmetric matrix preconditioner "
-                << preconditionerName << endl << endl
-                << "Valid symmetric matrix preconditioners are :" << endl
-                << symMatrixConstructorTablePtr_->toc()
-                << exit(FatalIOError);
-        }
-
-        return autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
-        (
-            constructorIter()
-            (
-                sol,
-                preconditionerDict
-            )
-        );
-    }
-    else if (sol.matrix().asymmetric())
-    {
-        typename asymMatrixConstructorTable::iterator constructorIter =
-            asymMatrixConstructorTablePtr_->find(preconditionerName);
-
-        if (constructorIter == asymMatrixConstructorTablePtr_->end())
-        {
-            FatalIOErrorIn
-            (
-                "LduMatrix<Type, DType, LUType>::preconditioner::New"
-                "(const solver&, Istream&)",
-                preconditionerDict
-            )   << "Unknown asymmetric matrix preconditioner "
-                << preconditionerName << endl << endl
-                << "Valid asymmetric matrix preconditioners are :" << endl
-                << asymMatrixConstructorTablePtr_->toc()
-                << exit(FatalIOError);
-        }
-
-        return autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
-        (
-            constructorIter()
-            (
-                sol,
-                preconditionerDict
-            )
-        );
-    }
-    else
-    {
-        FatalIOErrorIn
-        (
-            "LduMatrix<Type, DType, LUType>::preconditioner::New"
-            "(const solver&, Istream&)",
-            preconditionerDict
-        )   << "cannot preconditione incomplete matrix, "
-               "no diagonal or off-diagonal coefficient"
-            << exit(FatalIOError);
-
-        return autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
-        (
-            NULL
-        );
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSmoother.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSmoother.C
deleted file mode 100644
index 33f1ae28e0698c404f5b28363847055c0ab37ff9..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSmoother.C
+++ /dev/null
@@ -1,121 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::autoPtr<typename Foam::LduMatrix<Type, DType, LUType>::smoother>
-Foam::LduMatrix<Type, DType, LUType>::smoother::New
-(
-    const word& fieldName,
-    const LduMatrix<Type, DType, LUType>& matrix,
-    const dictionary& smootherDict
-)
-{
-    word smootherName = smootherDict.lookup("smoother");
-
-    if (matrix.symmetric())
-    {
-        typename symMatrixConstructorTable::iterator constructorIter =
-            symMatrixConstructorTablePtr_->find(smootherName);
-
-        if (constructorIter == symMatrixConstructorTablePtr_->end())
-        {
-            FatalIOErrorIn
-            (
-                "LduMatrix<Type, DType, LUType>::smoother::New", smootherDict
-            )   << "Unknown symmetric matrix smoother " << smootherName
-                << endl << endl
-                << "Valid symmetric matrix smoothers are :" << endl
-                << symMatrixConstructorTablePtr_->toc()
-                << exit(FatalIOError);
-        }
-
-        return autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
-        (
-            constructorIter()
-            (
-                fieldName,
-                matrix
-            )
-        );
-    }
-    else if (matrix.asymmetric())
-    {
-        typename asymMatrixConstructorTable::iterator constructorIter =
-            asymMatrixConstructorTablePtr_->find(smootherName);
-
-        if (constructorIter == asymMatrixConstructorTablePtr_->end())
-        {
-            FatalIOErrorIn
-            (
-                "LduMatrix<Type, DType, LUType>::smoother::New", smootherDict
-            )   << "Unknown asymmetric matrix smoother " << smootherName
-                << endl << endl
-                << "Valid asymmetric matrix smoothers are :" << endl
-                << asymMatrixConstructorTablePtr_->toc()
-                << exit(FatalIOError);
-        }
-
-        return autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
-        (
-            constructorIter()
-            (
-                fieldName,
-                matrix
-            )
-        );
-    }
-    else
-    {
-        FatalIOErrorIn
-        (
-            "LduMatrix<Type, DType, LUType>::smoother::New", smootherDict
-        )   << "cannot solve incomplete matrix, no off-diagonal coefficients"
-            << exit(FatalIOError);
-
-        return autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>(NULL);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::smoother::smoother
-(
-    const word& fieldName,
-    const LduMatrix<Type, DType, LUType>& matrix
-)
-:
-    fieldName_(fieldName),
-    matrix_(matrix)
-{}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C
deleted file mode 100644
index 2fbf05e5cde0ec1a5f8ee983f5222f58708805f3..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C
+++ /dev/null
@@ -1,191 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "LduMatrix.H"
-#include "DiagonalSolver.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::autoPtr<typename Foam::LduMatrix<Type, DType, LUType>::solver>
-Foam::LduMatrix<Type, DType, LUType>::solver::New
-(
-    const word& fieldName,
-    const LduMatrix<Type, DType, LUType>& matrix,
-    const dictionary& solverDict
-)
-{
-    word solverName = solverDict.lookup("solver");
-
-    if (matrix.diagonal())
-    {
-        return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>
-        (
-            new DiagonalSolver<Type, DType, LUType>
-            (
-                fieldName,
-                matrix,
-                solverDict
-            )
-        );
-    }
-    else if (matrix.symmetric())
-    {
-        typename symMatrixConstructorTable::iterator constructorIter =
-            symMatrixConstructorTablePtr_->find(solverName);
-
-        if (constructorIter == symMatrixConstructorTablePtr_->end())
-        {
-            FatalIOErrorIn
-            (
-                "LduMatrix<Type, DType, LUType>::solver::New", solverDict
-            )   << "Unknown symmetric matrix solver " << solverName
-                << endl << endl
-                << "Valid symmetric matrix solvers are :" << endl
-                << symMatrixConstructorTablePtr_->toc()
-                << exit(FatalIOError);
-        }
-
-        return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>
-        (
-            constructorIter()
-            (
-                fieldName,
-                matrix,
-                solverDict
-            )
-        );
-    }
-    else if (matrix.asymmetric())
-    {
-        typename asymMatrixConstructorTable::iterator constructorIter =
-            asymMatrixConstructorTablePtr_->find(solverName);
-
-        if (constructorIter == asymMatrixConstructorTablePtr_->end())
-        {
-            FatalIOErrorIn
-            (
-                "LduMatrix<Type, DType, LUType>::solver::New", solverDict
-            )   << "Unknown asymmetric matrix solver " << solverName
-                << endl << endl
-                << "Valid asymmetric matrix solvers are :" << endl
-                << asymMatrixConstructorTablePtr_->toc()
-                << exit(FatalIOError);
-        }
-
-        return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>
-        (
-            constructorIter()
-            (
-                fieldName,
-                matrix,
-                solverDict
-            )
-        );
-    }
-    else
-    {
-        FatalIOErrorIn
-        (
-            "LduMatrix<Type, DType, LUType>::solver::New", solverDict
-        )   << "cannot solve incomplete matrix, "
-               "no diagonal or off-diagonal coefficient"
-            << exit(FatalIOError);
-
-        return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>(NULL);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::LduMatrix<Type, DType, LUType>::solver::solver
-(
-    const word& fieldName,
-    const LduMatrix<Type, DType, LUType>& matrix,
-    const dictionary& solverDict
-)
-:
-    fieldName_(fieldName),
-    matrix_(matrix),
-
-    controlDict_(solverDict),
-
-    maxIter_(1000),
-    tolerance_(1e-6*pTraits<Type>::one),
-    relTol_(pTraits<Type>::zero)
-{
-    readControls();
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::solver::readControls()
-{
-    readControl(controlDict_, maxIter_, "maxIter");
-    readControl(controlDict_, tolerance_, "tolerance");
-    readControl(controlDict_, relTol_, "relTol");
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::solver::read
-(
-    const dictionary& solverDict
-)
-{
-    controlDict_ = solverDict;
-    readControls();
-}
-
-
-template<class Type, class DType, class LUType>
-Type Foam::LduMatrix<Type, DType, LUType>::solver::normFactor
-(
-    const Field<Type>& psi,
-    const Field<Type>& Apsi,
-    Field<Type>& tmpField
-) const
-{
-    // --- Calculate A dot reference value of psi
-    matrix_.sumA(tmpField);
-    cmptMultiply(tmpField, tmpField, gAverage(psi));
-
-    return stabilise
-    (
-        gSum(cmptMag(Apsi - tmpField) + cmptMag(matrix_.source() - tmpField)),
-        matrix_.small_
-    );
-
-    // At convergence this simpler method is equivalent to the above
-    // return stabilise(2*gSumCmptMag(matrix_.source()), matrix_.small_);
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C
deleted file mode 100644
index cb769fd0c891f3bd2d76fa0b8d38d6ed33b305f0..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C
+++ /dev/null
@@ -1,119 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular
-(
-    const Type& wApA
-)
-{
-    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-    {
-        singular_[cmpt] = component(wApA, cmpt) < vsmall_;
-    }
-
-    return singular();
-}
-
-
-template<class Type, class DType, class LUType>
-bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular() const
-{
-    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-    {
-        if (!singular_[cmpt]) return false;
-    }
-
-    return true;
-}
-
-
-template<class Type, class DType, class LUType>
-bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::converged
-(
-    const Type& Tolerance,
-    const Type& RelTolerance
-)
-{
-    if (debug >= 2)
-    {
-        Info<< solverName_
-            << ":  Iteration " << noIterations_
-            << " residual = " << finalResidual_
-            << endl;
-    }
-
-    if
-    (
-        finalResidual_ < Tolerance
-     || (
-            RelTolerance > small_*pTraits<Type>::one
-         && finalResidual_ < cmptMultiply(RelTolerance, initialResidual_)
-        )
-    )
-    {
-        converged_ = true;
-    }
-    else
-    {
-        converged_ = false;
-    }
-
-    return converged_;
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::solverPerformance::print
-(
-    Ostream& os
-) const
-{
-    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
-    {
-        os  << solverName_ << ":  Solving for "
-            << word(fieldName_ + pTraits<Type>::componentNames[cmpt]);
-
-        if (singular_[cmpt])
-        {
-            os  << ":  solution singularity" << endl;
-        }
-        else
-        {
-            os  << ", Initial residual = " << component(initialResidual_, cmpt)
-                << ", Final residual = " << component(finalResidual_, cmpt)
-                << ", No Iterations " << noIterations_
-                << endl;
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixUpdateMatrixInterfaces.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixUpdateMatrixInterfaces.C
deleted file mode 100644
index fbf608912ac33209e93546ef9b8ee6b69fa95653..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixUpdateMatrixInterfaces.C
+++ /dev/null
@@ -1,195 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "LduMatrix.H"
-#include "LduInterfaceField.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::initMatrixInterfaces
-(
-    const FieldField<Field, LUType>& interfaceCoeffs,
-    const Field<Type>& psiif,
-    Field<Type>& result
-) const
-{
-    if
-    (
-        Pstream::defaultCommsType == Pstream::blocking
-     || Pstream::defaultCommsType == Pstream::nonBlocking
-    )
-    {
-        forAll (interfaces_, interfaceI)
-        {
-            if (interfaces_.set(interfaceI))
-            {
-                interfaces_[interfaceI].initInterfaceMatrixUpdate
-                (
-                    result,
-                    psiif,
-                    Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
-                    Pstream::defaultCommsType
-                );
-            }
-        }
-    }
-    else if (Pstream::defaultCommsType == Pstream::scheduled)
-    {
-        const lduSchedule& patchSchedule = this->patchSchedule();
-
-        // Loop over the "global" patches are on the list of interfaces but
-        // beyond the end of the schedule which only handles "normal" patches
-        for
-        (
-            label interfaceI=patchSchedule.size()/2;
-            interfaceI<interfaces_.size();
-            interfaceI++
-        )
-        {
-            if (interfaces_.set(interfaceI))
-            {
-                interfaces_[interfaceI].initInterfaceMatrixUpdate
-                (
-                    result,
-                    psiif,
-                    Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
-                    Pstream::blocking
-                );
-            }
-        }
-    }
-    else
-    {
-        FatalErrorIn("LduMatrix<Type, DType, LUType>::initMatrixInterfaces")
-            << "Unsuported communications type "
-            << Pstream::commsTypeNames[Pstream::defaultCommsType]
-            << exit(FatalError);
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
-(
-    const FieldField<Field, LUType>& interfaceCoeffs,
-    const Field<Type>& psiif,
-    Field<Type>& result
-) const
-{
-    if
-    (
-        Pstream::defaultCommsType == Pstream::blocking
-     || Pstream::defaultCommsType == Pstream::nonBlocking
-    )
-    {
-        // Block until all sends/receives have been finished
-        if (Pstream::defaultCommsType == Pstream::nonBlocking)
-        {
-            IPstream::waitRequests();
-            OPstream::waitRequests();
-        }
-
-        forAll (interfaces_, interfaceI)
-        {
-            if (interfaces_.set(interfaceI))
-            {
-                interfaces_[interfaceI].updateInterfaceMatrix
-                (
-                    result,
-                    psiif,
-                    Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
-                    Pstream::defaultCommsType
-                );
-            }
-        }
-    }
-    else if (Pstream::defaultCommsType == Pstream::scheduled)
-    {
-        const lduSchedule& patchSchedule = this->patchSchedule();
-
-        // Loop over all the "normal" interfaces relating to standard patches
-        forAll (patchSchedule, i)
-        {
-            label interfaceI = patchSchedule[i].patch;
-
-            if (interfaces_.set(interfaceI))
-            {
-                if (patchSchedule[i].init)
-                {
-                    interfaces_[interfaceI].initInterfaceMatrixUpdate
-                    (
-                        result,
-                        psiif,
-                        Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
-                        Pstream::scheduled
-                    );
-                }
-                else
-                {
-                    interfaces_[interfaceI].updateInterfaceMatrix
-                    (
-                        result,
-                        psiif,
-                        Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
-                        Pstream::scheduled
-                    );
-                }
-            }
-        }
-
-        // Loop over the "global" patches are on the list of interfaces but
-        // beyond the end of the schedule which only handles "normal" patches
-        for
-        (
-            label interfaceI=patchSchedule.size()/2;
-            interfaceI<interfaces_.size();
-            interfaceI++
-        )
-        {
-            if (interfaces_.set(interfaceI))
-            {
-                interfaces_[interfaceI].updateInterfaceMatrix
-                (
-                    result,
-                    psiif,
-                    Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
-                    Pstream::blocking
-                );
-            }
-        }
-    }
-    else
-    {
-        FatalErrorIn("LduMatrix<Type, DType, LUType>::updateMatrixInterfaces")
-            << "Unsuported communications type "
-            << Pstream::commsTypeNames[Pstream::defaultCommsType]
-            << exit(FatalError);
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/lduMatrices.C b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/lduMatrices.C
deleted file mode 100644
index 7155c8d993c5cfe70790ef1c637384353501a8de..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/lduMatrices.C
+++ /dev/null
@@ -1,11 +0,0 @@
-#include "LduMatrix.H"
-#include "fieldTypes.H"
-
-namespace Foam
-{
-    makeLduMatrix(scalar, scalar, scalar);
-    makeLduMatrix(vector, scalar, scalar);
-    makeLduMatrix(sphericalTensor, scalar, scalar);
-    makeLduMatrix(symmTensor, scalar, scalar);
-    makeLduMatrix(tensor, scalar, scalar);
-};
diff --git a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DILUPreconditioner/TDILUPreconditioner.C b/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DILUPreconditioner/TDILUPreconditioner.C
deleted file mode 100644
index c7c5c40743eff8a4d8f5f5968317261367fa7ec3..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DILUPreconditioner/TDILUPreconditioner.C
+++ /dev/null
@@ -1,180 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "TDILUPreconditioner.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::TDILUPreconditioner<Type, DType, LUType>::TDILUPreconditioner
-(
-    const typename LduMatrix<Type, DType, LUType>::solver& sol,
-    const dictionary&
-)
-:
-    LduMatrix<Type, DType, LUType>::preconditioner(sol),
-    rD_(sol.matrix().diag())
-{
-    calcInvD(rD_, sol.matrix());
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::TDILUPreconditioner<Type, DType, LUType>::calcInvD
-(
-    Field<DType>& rD,
-    const LduMatrix<Type, DType, LUType>& matrix
-)
-{
-    DType* __restrict__ rDPtr = rD.begin();
-
-    const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
-    const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();
-
-    const LUType* const __restrict__ upperPtr = matrix.upper().begin();
-    const LUType* const __restrict__ lowerPtr = matrix.lower().begin();
-
-    register label nFaces = matrix.upper().size();
-    for (register label face=0; face<nFaces; face++)
-    {
-        rDPtr[uPtr[face]] -= 
-            dot(dot(upperPtr[face], lowerPtr[face]), inv(rDPtr[lPtr[face]]));
-    }
-
-
-    // Calculate the reciprocal of the preconditioned diagonal
-    register label nCells = rD.size();
-
-    for (register label cell=0; cell<nCells; cell++)
-    {
-        rDPtr[cell] = inv(rDPtr[cell]);
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::TDILUPreconditioner<Type, DType, LUType>::precondition
-(
-    Field<Type>& wA,
-    const Field<Type>& rA
-) const
-{
-    Type* __restrict__ wAPtr = wA.begin();
-    const Type* __restrict__ rAPtr = rA.begin();
-    const DType* __restrict__ rDPtr = rD_.begin();
-
-    const label* const __restrict__ uPtr =
-        this->solver_.matrix().lduAddr().upperAddr().begin();
-    const label* const __restrict__ lPtr =
-        this->solver_.matrix().lduAddr().lowerAddr().begin();
-    const label* const __restrict__ losortPtr =
-        this->solver_.matrix().lduAddr().losortAddr().begin();
-
-    const LUType* const __restrict__ upperPtr =
-        this->solver_.matrix().upper().begin();
-    const LUType* const __restrict__ lowerPtr =
-        this->solver_.matrix().lower().begin();
-
-    register label nCells = wA.size();
-    register label nFaces = this->solver_.matrix().upper().size();
-    register label nFacesM1 = nFaces - 1;
-
-    for (register label cell=0; cell<nCells; cell++)
-    {
-        wAPtr[cell] = dot(rDPtr[cell], rAPtr[cell]);
-    }
-
-
-    register label sface;
-
-    for (register label face=0; face<nFaces; face++)
-    {
-        sface = losortPtr[face];
-        wAPtr[uPtr[sface]] -=
-            dot(rDPtr[uPtr[sface]], dot(lowerPtr[sface], wAPtr[lPtr[sface]]));
-    }
-
-    for (register label face=nFacesM1; face>=0; face--)
-    {
-        wAPtr[lPtr[face]] -=
-            dot(rDPtr[lPtr[face]], dot(upperPtr[face], wAPtr[uPtr[face]]));
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::TDILUPreconditioner<Type, DType, LUType>::preconditionT
-(
-    Field<Type>& wT,
-    const Field<Type>& rT
-) const
-{
-    Type* __restrict__ wTPtr = wT.begin();
-    const Type* __restrict__ rTPtr = rT.begin();
-    const DType* __restrict__ rDPtr = rD_.begin();
-
-    const label* const __restrict__ uPtr =
-        this->solver_.matrix().lduAddr().upperAddr().begin();
-    const label* const __restrict__ lPtr =
-        this->solver_.matrix().lduAddr().lowerAddr().begin();
-    const label* const __restrict__ losortPtr =
-        this->solver_.matrix().lduAddr().losortAddr().begin();
-
-    const LUType* const __restrict__ upperPtr =
-        this->solver_.matrix().upper().begin();
-    const LUType* const __restrict__ lowerPtr =
-        this->solver_.matrix().lower().begin();
-
-    register label nCells = wT.size();
-    register label nFaces = this->solver_.matrix().upper().size();
-    register label nFacesM1 = nFaces - 1;
-
-    for (register label cell=0; cell<nCells; cell++)
-    {
-        wTPtr[cell] = dot(rDPtr[cell], rTPtr[cell]);
-    }
-
-    for (register label face=0; face<nFaces; face++)
-    {
-        wTPtr[uPtr[face]] -=
-            dot(rDPtr[uPtr[face]], dot(upperPtr[face], wTPtr[lPtr[face]]));
-    }
-
-
-    register label sface;
-
-    for (register label face=nFacesM1; face>=0; face--)
-    {
-        sface = losortPtr[face];
-        wTPtr[lPtr[sface]] -=
-            dot(rDPtr[lPtr[sface]], dot(lowerPtr[sface], wTPtr[uPtr[sface]]));
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DILUPreconditioner/TDILUPreconditioner.H b/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DILUPreconditioner/TDILUPreconditioner.H
deleted file mode 100644
index 73a8bf11d62bf46e8fb403829a80edebadd7d7c3..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DILUPreconditioner/TDILUPreconditioner.H
+++ /dev/null
@@ -1,127 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::TDILUPreconditioner
-
-Description
-    Simplified diagonal-based incomplete LU preconditioner for asymmetric
-    matrices.
-
-    The inverse (reciprocal for scalar) of the preconditioned diagonal is
-    calculated and stored.
-
-SourceFiles
-    TDILUPreconditioner.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef TDILUPreconditioner_H
-#define TDILUPreconditioner_H
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class TDILUPreconditioner Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type, class DType, class LUType>
-class TDILUPreconditioner
-:
-    public LduMatrix<Type, DType, LUType>::preconditioner
-{
-    // Private data
-
-        //- The inverse (reciprocal for scalar) preconditioned diagonal
-        Field<DType> rD_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("DILU");
-
-
-    // Constructors
-
-        //- Construct from matrix components and preconditioner data dictionary
-        TDILUPreconditioner
-        (
-            const typename LduMatrix<Type, DType, LUType>::solver& sol,
-            const dictionary& preconditionerDict
-        );
-
-
-    // Destructor
-
-        virtual ~TDILUPreconditioner()
-        {}
-
-
-    // Member Functions
-
-        //- Calculate the reciprocal of the preconditioned diagonal
-        static void calcInvD
-        (
-            Field<DType>& rD,
-            const LduMatrix<Type, DType, LUType>& matrix
-        );
-
-        //- Return wA the preconditioned form of residual rA
-        virtual void precondition
-        (
-            Field<Type>& wA,
-            const Field<Type>& rA
-        ) const;
-
-        //- Return wT the transpose-matrix preconditioned form of
-        //  residual rT.
-        virtual void preconditionT
-        (
-            Field<Type>& wT,
-            const Field<Type>& rT
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "TDILUPreconditioner.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DiagonalPreconditioner/DiagonalPreconditioner.C b/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DiagonalPreconditioner/DiagonalPreconditioner.C
deleted file mode 100644
index bb085774062dffb5ae77c5334225831d5cf5e0c5..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DiagonalPreconditioner/DiagonalPreconditioner.C
+++ /dev/null
@@ -1,81 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "DiagonalPreconditioner.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::DiagonalPreconditioner<Type, DType, LUType>::DiagonalPreconditioner
-(
-    const typename LduMatrix<Type, DType, LUType>::solver& sol,
-    const dictionary&
-)
-:
-    LduMatrix<Type, DType, LUType>::preconditioner(sol),
-    rD(sol.matrix().diag().size())
-{
-    DType* __restrict__ rDPtr = rD.begin();
-    const DType* __restrict__ DPtr = this->solver_.matrix().diag().begin();
-
-    register label nCells = rD.size();
-
-    // Generate inverse (reciprocal for scalar) diagonal
-    for (register label cell=0; cell<nCells; cell++)
-    {
-        rDPtr[cell] = inv(DPtr[cell]);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::DiagonalPreconditioner<Type, DType, LUType>::read(const dictionary&)
-{}
-
-
-template<class Type, class DType, class LUType>
-void Foam::DiagonalPreconditioner<Type, DType, LUType>::precondition
-(
-    Field<Type>& wA,
-    const Field<Type>& rA
-) const
-{
-    Type* __restrict__ wAPtr = wA.begin();
-    const Type* __restrict__ rAPtr = rA.begin();
-    const DType* __restrict__ rDPtr = rD.begin();
-
-    register label nCells = wA.size();
-
-    for (register label cell=0; cell<nCells; cell++)
-    {
-        wAPtr[cell] = dot(rDPtr[cell], rAPtr[cell]);
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DiagonalPreconditioner/DiagonalPreconditioner.H b/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DiagonalPreconditioner/DiagonalPreconditioner.H
deleted file mode 100644
index b261a52ab95730156f38d97681a0ab1cd86a2a08..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/DiagonalPreconditioner/DiagonalPreconditioner.H
+++ /dev/null
@@ -1,135 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::DiagonalPreconditioner
-
-Description
-    Diagonal preconditioner for both symmetric and asymmetric matrices.
-
-    The inverse (reciprocal for scalar) of the diagonal is calculated and
-    stored.
-
-SourceFiles
-    DiagonalPreconditioner.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef DiagonalPreconditioner_H
-#define DiagonalPreconditioner_H
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class DiagonalPreconditioner Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type, class DType, class LUType>
-class DiagonalPreconditioner
-:
-    public LduMatrix<Type, DType, LUType>::preconditioner
-{
-    // Private data
-
-        //- The inverse (reciprocal for scalar) diagonal 
-        Field<DType> rD;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        DiagonalPreconditioner(const DiagonalPreconditioner&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const DiagonalPreconditioner&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("diagonal");
-
-
-    // Constructors
-
-        //- Construct from matrix components and preconditioner data dictionary
-        DiagonalPreconditioner
-        (
-            const typename LduMatrix<Type, DType, LUType>::solver& sol,
-            const dictionary& preconditionerDict
-        );
-
-
-    // Destructor
-
-        virtual ~DiagonalPreconditioner()
-        {}
-
-
-    // Member Functions
-
-        //- Read and reset the preconditioner parameters from the given
-        //  dictionary
-        virtual void read(const dictionary& preconditionerDict);
-
-        //- Return wA the preconditioned form of residual rA
-        virtual void precondition
-        (
-            Field<Type>& wA,
-            const Field<Type>& rA
-        ) const;
-
-        //- Return wT the transpose-matrix preconditioned form of
-        //  residual rT.
-        virtual void preconditionT
-        (
-            Field<Type>& wT,
-            const Field<Type>& rT
-        ) const
-        {
-            return(precondition(wT, rT));
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "DiagonalPreconditioner.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/NoPreconditioner/NoPreconditioner.C b/src/OpenFOAM/matrices/LduMatrix/Preconditioners/NoPreconditioner/NoPreconditioner.C
deleted file mode 100644
index d3a910cec821362eabec95f84c9fac61d0013753..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/NoPreconditioner/NoPreconditioner.C
+++ /dev/null
@@ -1,60 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "NoPreconditioner.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::NoPreconditioner<Type, DType, LUType>::NoPreconditioner
-(
-    const typename LduMatrix<Type, DType, LUType>::solver& sol,
-    const dictionary&
-)
-:
-    LduMatrix<Type, DType, LUType>::preconditioner(sol)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::NoPreconditioner<Type, DType, LUType>::read(const dictionary&)
-{}
-
-
-template<class Type, class DType, class LUType>
-void Foam::NoPreconditioner<Type, DType, LUType>::precondition
-(
-    Field<Type>& wA,
-    const Field<Type>& rA
-) const
-{
-    wA = rA;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/NoPreconditioner/NoPreconditioner.H b/src/OpenFOAM/matrices/LduMatrix/Preconditioners/NoPreconditioner/NoPreconditioner.H
deleted file mode 100644
index f139efd311e6bac545b784cc89e5fe4e4930a99e..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/NoPreconditioner/NoPreconditioner.H
+++ /dev/null
@@ -1,126 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::NoPreconditioner
-
-Description
-    Null preconditioner for both symmetric and asymmetric matrices.
-
-SourceFiles
-    NoPreconditioner.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef NoPreconditioner_H
-#define NoPreconditioner_H
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class NoPreconditioner Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type, class DType, class LUType>
-class NoPreconditioner
-:
-    public LduMatrix<Type, DType, LUType>::preconditioner
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        NoPreconditioner(const NoPreconditioner&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const NoPreconditioner&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("none");
-
-
-    // Constructors
-
-        //- Construct from matrix components and preconditioner data dictionary
-        NoPreconditioner
-        (
-            const typename LduMatrix<Type, DType, LUType>::solver& sol,
-            const dictionary& preconditionerDict
-        );
-
-
-    // Destructor
-
-        virtual ~NoPreconditioner()
-        {}
-
-
-    // Member Functions
-
-        //- Read and reset the preconditioner parameters from the given
-        //  dictionary
-        virtual void read(const dictionary& preconditionerDict);
-
-        //- Return wA the preconditioned form of residual rA
-        virtual void precondition
-        (
-            Field<Type>& wA,
-            const Field<Type>& rA
-        ) const;
-
-        //- Return wT the transpose-matrix preconditioned form of
-        //  residual rT.
-        virtual void preconditionT
-        (
-            Field<Type>& wT,
-            const Field<Type>& rT
-        ) const
-        {
-            return(precondition(wT, rT));
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "NoPreconditioner.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/lduPreconditioners.C b/src/OpenFOAM/matrices/LduMatrix/Preconditioners/lduPreconditioners.C
deleted file mode 100644
index de2a717d68c1aaf4c81cffb3f98cec1188d2379d..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Preconditioners/lduPreconditioners.C
+++ /dev/null
@@ -1,26 +0,0 @@
-#include "NoPreconditioner.H"
-#include "DiagonalPreconditioner.H"
-#include "TDILUPreconditioner.H"
-#include "fieldTypes.H"
-
-#define makeLduPreconditioners(Type, DType, LUType)                           \
-                                                                              \
-    makeLduPreconditioner(NoPreconditioner, Type, DType, LUType);             \
-    makeLduSymPreconditioner(NoPreconditioner, Type, DType, LUType);          \
-    makeLduAsymPreconditioner(NoPreconditioner, Type, DType, LUType);         \
-                                                                              \
-    makeLduPreconditioner(DiagonalPreconditioner, Type, DType, LUType);       \
-    makeLduSymPreconditioner(DiagonalPreconditioner, Type, DType, LUType);    \
-    makeLduAsymPreconditioner(DiagonalPreconditioner, Type, DType, LUType);   \
-                                                                              \
-    makeLduPreconditioner(TDILUPreconditioner, Type, DType, LUType);          \
-    makeLduAsymPreconditioner(TDILUPreconditioner, Type, DType, LUType);
-
-namespace Foam
-{
-    makeLduPreconditioners(scalar, scalar, scalar);
-    makeLduPreconditioners(vector, scalar, scalar);
-    makeLduPreconditioners(sphericalTensor, scalar, scalar);
-    makeLduPreconditioners(symmTensor, scalar, scalar);
-    makeLduPreconditioners(tensor, scalar, scalar);
-};
diff --git a/src/OpenFOAM/matrices/LduMatrix/Smoothers/GaussSeidel/TGaussSeidelSmoother.C b/src/OpenFOAM/matrices/LduMatrix/Smoothers/GaussSeidel/TGaussSeidelSmoother.C
deleted file mode 100644
index d23644f0e26c118c1aa870887c6e76f38fd573a5..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Smoothers/GaussSeidel/TGaussSeidelSmoother.C
+++ /dev/null
@@ -1,168 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "TGaussSeidelSmoother.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::TGaussSeidelSmoother<Type, DType, LUType>::TGaussSeidelSmoother
-(
-    const word& fieldName,
-    const LduMatrix<Type, DType, LUType>& matrix
-)
-:
-    LduMatrix<Type, DType, LUType>::smoother
-    (
-        fieldName,
-        matrix
-    ),
-    rD_(matrix.diag().size())
-{
-    register const label nCells = matrix.diag().size();
-    register const DType* const __restrict__ diagPtr = matrix.diag().begin();
-    register DType* __restrict__ rDPtr = rD_.begin();
-
-    for (register label cellI=0; cellI<nCells; cellI++)
-    {
-        rDPtr[cellI] = inv(diagPtr[cellI]);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::TGaussSeidelSmoother<Type, DType, LUType>::smooth
-(
-    const word& fieldName_,
-    Field<Type>& psi,
-    const LduMatrix<Type, DType, LUType>& matrix_,
-    const Field<DType>& rD_,
-    const label nSweeps
-)
-{
-    register Type* __restrict__ psiPtr = psi.begin();
-
-    register const label nCells = psi.size();
-
-    Field<Type> bPrime(nCells);
-    register Type* __restrict__ bPrimePtr = bPrime.begin();
-
-    register const DType* const __restrict__ rDPtr = rD_.begin();
-
-    register const LUType* const __restrict__ upperPtr =
-        matrix_.upper().begin();
-
-    register const LUType* const __restrict__ lowerPtr =
-        matrix_.lower().begin();
-
-    register const label* const __restrict__ uPtr =
-        matrix_.lduAddr().upperAddr().begin();
-
-    register const label* const __restrict__ ownStartPtr =
-        matrix_.lduAddr().ownerStartAddr().begin();
-
-
-    // Parallel boundary initialisation.  The parallel boundary is treated
-    // as an effective jacobi interface in the boundary.
-    // Note: there is a change of sign in the coupled
-    // interface update to add the contibution to the r.h.s.
-
-    FieldField<Field, LUType> mBouCoeffs(matrix_.interfacesUpper().size());
-
-    forAll(mBouCoeffs, patchi)
-    {
-        if (matrix_.interfaces().set(patchi))
-        {
-            mBouCoeffs.set(patchi, -matrix_.interfacesUpper()[patchi]);
-        }
-    }
-
-    for (label sweep=0; sweep<nSweeps; sweep++)
-    {
-        bPrime = matrix_.source();
-
-        matrix_.initMatrixInterfaces
-        (
-            mBouCoeffs,
-            psi,
-            bPrime
-        );
-
-        matrix_.updateMatrixInterfaces
-        (
-            mBouCoeffs,
-            psi,
-            bPrime
-        );
-
-        Type curPsi;
-        register label fStart;
-        register label fEnd = ownStartPtr[0];
-
-        for (register label cellI=0; cellI<nCells; cellI++)
-        {
-            // Start and end of this row
-            fStart = fEnd;
-            fEnd = ownStartPtr[cellI + 1];
-
-            // Get the accumulated neighbour side
-            curPsi = bPrimePtr[cellI];
-
-            // Accumulate the owner product side
-            for (register label curFace=fStart; curFace<fEnd; curFace++)
-            {
-                curPsi -= dot(upperPtr[curFace], psiPtr[uPtr[curFace]]);
-            }
-
-            // Finish current psi
-            curPsi = dot(rDPtr[cellI], curPsi);
-
-            // Distribute the neighbour side using current psi
-            for (register label curFace=fStart; curFace<fEnd; curFace++)
-            {
-                bPrimePtr[uPtr[curFace]] -= dot(lowerPtr[curFace], curPsi);
-            }
-
-            psiPtr[cellI] = curPsi;
-        }
-    }
-}
-
-
-template<class Type, class DType, class LUType>
-void Foam::TGaussSeidelSmoother<Type, DType, LUType>::smooth
-(
-    Field<Type>& psi,
-    const label nSweeps
-) const
-{
-    smooth(this->fieldName_, psi, this->matrix_, rD_, nSweeps);
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Smoothers/GaussSeidel/TGaussSeidelSmoother.H b/src/OpenFOAM/matrices/LduMatrix/Smoothers/GaussSeidel/TGaussSeidelSmoother.H
deleted file mode 100644
index 8fd452b8161f35f1393458f6c4be41e8f2702f6a..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Smoothers/GaussSeidel/TGaussSeidelSmoother.H
+++ /dev/null
@@ -1,113 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::TGaussSeidelSmoother
-
-Description
-    Foam::TGaussSeidelSmoother
-
-SourceFiles
-    TGaussSeidelSmoother.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef TGaussSeidelSmoother_H
-#define TGaussSeidelSmoother_H
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class TGaussSeidelSmoother Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type, class DType, class LUType>
-class TGaussSeidelSmoother
-:
-    public LduMatrix<Type, DType, LUType>::smoother
-{
-    // Private data
-
-        //- The inverse (reciprocal for scalars) diagonal
-        Field<DType> rD_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("GaussSeidel");
-
-
-    // Constructors
-
-        //- Construct from components
-        TGaussSeidelSmoother
-        (
-            const word& fieldName,
-            const LduMatrix<Type, DType, LUType>& matrix
-        );
-
-
-    // Member Functions
-
-        //- Smooth for the given number of sweeps
-        static void smooth
-        (
-            const word& fieldName,
-            Field<Type>& psi,
-            const LduMatrix<Type, DType, LUType>& matrix,
-            const Field<DType>& rD,
-            const label nSweeps
-        );
-
-
-        //- Smooth the solution for a given number of sweeps
-        virtual void smooth
-        (
-            Field<Type>& psi,
-            const label nSweeps
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "TGaussSeidelSmoother.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Smoothers/lduSmoothers.C b/src/OpenFOAM/matrices/LduMatrix/Smoothers/lduSmoothers.C
deleted file mode 100644
index 736963e656f39f97e566b1280f84aab16059d15b..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Smoothers/lduSmoothers.C
+++ /dev/null
@@ -1,17 +0,0 @@
-#include "TGaussSeidelSmoother.H"
-#include "fieldTypes.H"
-
-#define makeLduSmoothers(Type, DType, LUType)                                 \
-                                                                              \
-    makeLduSmoother(TGaussSeidelSmoother, Type, DType, LUType);               \
-    makeLduSymSmoother(TGaussSeidelSmoother, Type, DType, LUType);            \
-    makeLduAsymSmoother(TGaussSeidelSmoother, Type, DType, LUType);
-
-namespace Foam
-{
-    makeLduSmoothers(scalar, scalar, scalar);
-    makeLduSmoothers(vector, scalar, scalar);
-    makeLduSmoothers(sphericalTensor, scalar, scalar);
-    makeLduSmoothers(symmTensor, scalar, scalar);
-    makeLduSmoothers(tensor, scalar, scalar);
-};
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C
deleted file mode 100644
index 5d9daa6aaadc9162c47a92a06770610661158eac..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C
+++ /dev/null
@@ -1,80 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "DiagonalSolver.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::DiagonalSolver<Type, DType, LUType>::DiagonalSolver
-(
-    const word& fieldName,
-    const LduMatrix<Type, DType, LUType>& matrix,
-    const dictionary& solverDict
-)
-:
-    LduMatrix<Type, DType, LUType>::solver
-    (
-        fieldName,
-        matrix,
-        solverDict
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::DiagonalSolver<Type, DType, LUType>::read
-(
-    const dictionary&
-)
-{}
-
-
-template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
-Foam::DiagonalSolver<Type, DType, LUType>::solve
-(
-    Field<Type>& psi
-) const
-{
-    psi = this->matrix_.source()/this->matrix_.diag();
-
-    return typename LduMatrix<Type, DType, LUType>::solverPerformance
-    (
-        typeName,
-        this->fieldName_,
-        pTraits<Type>::zero,
-        pTraits<Type>::zero,
-        0,
-        true,
-        false
-    );
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H
deleted file mode 100644
index b4c3d6ba8ddb7b7b78340e7f0d91c4128f2cc23b..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H
+++ /dev/null
@@ -1,108 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::DiagonalSolver
-
-Description
-    Foam::DiagonalSolver
-
-SourceFiles
-    DiagonalSolver.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef DiagonalSolver_H
-#define DiagonalSolver_H
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class DiagonalSolver Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type, class DType, class LUType>
-class DiagonalSolver
-:
-    public LduMatrix<Type, DType, LUType>::solver
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        DiagonalSolver(const DiagonalSolver&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const DiagonalSolver&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("diagonal");
-
-
-    // Constructors
-
-        //- Construct from matrix
-        DiagonalSolver
-        (
-            const word& fieldName,
-            const LduMatrix<Type, DType, LUType>& matrix,
-            const dictionary& solverDict
-        );
-
-
-    // Member Functions
-
-        //- Read and reset the solver parameters from the given dictionary
-        void read(const dictionary& solverDict);
-
-        //- Solve the matrix with this solver
-        typename LduMatrix<Type, DType, LUType>::solverPerformance solve
-        (
-            Field<Type>& psi
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "DiagonalSolver.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCG/TPBiCG.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCG/TPBiCG.C
deleted file mode 100644
index 927a8b2bb71e41f8fb1b3b254ff184b491a63947..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCG/TPBiCG.C
+++ /dev/null
@@ -1,195 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "TPBiCG.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::TPBiCG<Type, DType, LUType>::TPBiCG
-(
-    const word& fieldName,
-    const LduMatrix<Type, DType, LUType>& matrix,
-    const dictionary& solverDict
-)
-:
-    LduMatrix<Type, DType, LUType>::solver
-    (
-        fieldName,
-        matrix,
-        solverDict
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
-Foam::TPBiCG<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
-    (
-        preconditionerName + typeName,
-        this->fieldName_
-    );
-
-    register label nCells = psi.size();
-
-    Type* __restrict__ psiPtr = psi.begin();
-
-    Field<Type> pA(nCells);
-    Type* __restrict__ pAPtr = pA.begin();
-
-    Field<Type> pT(nCells, pTraits<Type>::zero);
-    Type* __restrict__ pTPtr = pT.begin();
-
-    Field<Type> wA(nCells);
-    Type* __restrict__ wAPtr = wA.begin();
-
-    Field<Type> wT(nCells);
-    Type* __restrict__ wTPtr = wT.begin();
-
-    Type wArT = this->matrix_.great_*pTraits<Type>::one;
-    Type wArTold = wArT;
-
-    // --- Calculate A.psi and T.psi
-    this->matrix_.Amul(wA, psi);
-    this->matrix_.Tmul(wT, psi);
-
-    // --- Calculate initial residual and transpose residual fields
-    Field<Type> rA(this->matrix_.source() - wA);
-    Field<Type> rT(this->matrix_.source() - wT);
-    Type* __restrict__ rAPtr = rA.begin();
-    Type* __restrict__ rTPtr = rT.begin();
-
-    // --- Calculate normalisation factor
-    Type normFactor = this->normFactor(psi, wA, pA);
-
-    if (LduMatrix<Type, DType, LUType>::debug >= 2)
-    {
-        Info<< "   Normalisation factor = " << normFactor << endl;
-    }
-
-    // --- Calculate normalised residual norm
-    solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
-    solverPerf.finalResidual() = solverPerf.initialResidual();
-
-    // --- Check convergence, solve if not converged
-    if (!solverPerf.converged(this->tolerance_, this->relTol_))
-    {
-        // --- Select and construct the preconditioner
-        autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
-        preconPtr = LduMatrix<Type, DType, LUType>::preconditioner::New
-        (
-            *this,
-            this->controlDict_
-        );
-
-        // --- Solver iteration
-        do
-        {
-            // --- Store previous wArT
-            wArTold = wArT;
-
-            // --- Precondition residuals
-            preconPtr->precondition(wA, rA);
-            preconPtr->preconditionT(wT, rT);
-
-            // --- Update search directions:
-            wArT = gSumCmptProd(wA, rT);
-
-            if (solverPerf.nIterations() == 0)
-            {
-                for (register label cell=0; cell<nCells; cell++)
-                {
-                    pAPtr[cell] = wAPtr[cell];
-                    pTPtr[cell] = wTPtr[cell];
-                }
-            }
-            else
-            {
-                Type beta = cmptDivide
-                (
-                    wArT,
-                    stabilise(wArTold, this->matrix_.vsmall_)
-                );
-
-                for (register label cell=0; cell<nCells; cell++)
-                {
-                    pAPtr[cell] = wAPtr[cell] + cmptMultiply(beta, pAPtr[cell]);
-                    pTPtr[cell] = wTPtr[cell] + cmptMultiply(beta, pTPtr[cell]);
-                }
-            }
-
-
-            // --- Update preconditioned residuals
-            this->matrix_.Amul(wA, pA);
-
-            this->matrix_.Tmul(wT, pT);
-
-            Type wApT = gSumCmptProd(wA, pT);
-
-            // --- Test for singularity
-            if (solverPerf.singular(cmptDivide(cmptMag(wApT), normFactor)))
-            {
-                break;
-            }
-
-
-            // --- Update solution and residual:
-
-            Type alpha = cmptDivide
-            (
-                wArT,
-                stabilise(wApT, this->matrix_.vsmall_)
-            );
-
-            for (register label cell=0; cell<nCells; cell++)
-            {
-                psiPtr[cell] += cmptMultiply(alpha, pAPtr[cell]);
-                rAPtr[cell] -= cmptMultiply(alpha, wAPtr[cell]);
-                rTPtr[cell] -= cmptMultiply(alpha, wTPtr[cell]);
-            }
-
-            solverPerf.finalResidual() =
-                cmptDivide(gSumCmptMag(rA), normFactor);
-
-        } while
-        (
-            solverPerf.nIterations()++ < this->maxIter_
-        && !(solverPerf.converged(this->tolerance_, this->relTol_))
-        );
-    }
-
-    return solverPerf;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCG/TPBiCG.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCG/TPBiCG.H
deleted file mode 100644
index d75e6e3b8a4e3108176f841cdf8956d63a58d14a..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCG/TPBiCG.H
+++ /dev/null
@@ -1,112 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::TPBiCG
-
-Description
-    Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices
-    using a run-time selectable preconditiioner.
-
-SourceFiles
-    TPBiCG.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef TPBiCG_H
-#define TPBiCG_H
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class TPBiCG Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type, class DType, class LUType>
-class TPBiCG
-:
-    public LduMatrix<Type, DType, LUType>::solver
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        TPBiCG(const TPBiCG&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const TPBiCG&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("PBiCG");
-
-
-    // Constructors
-
-        //- Construct from matrix components and solver data dictionary
-        TPBiCG
-        (
-            const word& fieldName,
-            const LduMatrix<Type, DType, LUType>& matrix,
-            const dictionary& solverDict
-        );
-
-
-    // Destructor
-
-        virtual ~TPBiCG()
-        {}
-
-
-    // Member Functions
-
-        //- Solve the matrix with this solver
-        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
-        (
-            Field<Type>& psi
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "TPBiCG.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCGScalarAlpha/PBiCGScalarAlpha.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCGScalarAlpha/PBiCGScalarAlpha.C
deleted file mode 100644
index fd02cf7a739db88c8e5641e2006ada1b2c805669..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCGScalarAlpha/PBiCGScalarAlpha.C
+++ /dev/null
@@ -1,191 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "TPBiCG.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::TPBiCG<Type, DType, LUType>::TPBiCG
-(
-    const word& fieldName,
-    const LduMatrix<Type, DType, LUType>& matrix,
-    const dictionary& solverDict
-)
-:
-    LduMatrix<Type, DType, LUType>::solver
-    (
-        fieldName,
-        matrix,
-        solverDict
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
-Foam::TPBiCG<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
-    (
-        preconditionerName + typeName,
-        this->fieldName_
-    );
-
-    register label nCells = psi.size();
-
-    Type* __restrict__ psiPtr = psi.begin();
-
-    Field<Type> pA(nCells);
-    Type* __restrict__ pAPtr = pA.begin();
-
-    Field<Type> pT(nCells, pTraits<Type>::zero);
-    Type* __restrict__ pTPtr = pT.begin();
-
-    Field<Type> wA(nCells);
-    Type* __restrict__ wAPtr = wA.begin();
-
-    Field<Type> wT(nCells);
-    Type* __restrict__ wTPtr = wT.begin();
-
-    scalar wArT = 1e15; //this->matrix_.great_;
-    scalar wArTold = wArT;
-
-    // --- Calculate A.psi and T.psi
-    this->matrix_.Amul(wA, psi);
-    this->matrix_.Tmul(wT, psi);
-
-    // --- Calculate initial residual and transpose residual fields
-    Field<Type> rA(this->matrix_.source() - wA);
-    Field<Type> rT(this->matrix_.source() - wT);
-    Type* __restrict__ rAPtr = rA.begin();
-    Type* __restrict__ rTPtr = rT.begin();
-
-    // --- Calculate normalisation factor
-    Type normFactor = this->normFactor(psi, wA, pA);
-
-    if (LduMatrix<Type, DType, LUType>::debug >= 2)
-    {
-        Info<< "   Normalisation factor = " << normFactor << endl;
-    }
-
-    // --- Calculate normalised residual norm
-    solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
-    solverPerf.finalResidual() = solverPerf.initialResidual();
-
-    // --- Check convergence, solve if not converged
-    if (!solverPerf.converged(this->tolerance_, this->relTol_))
-    {
-        // --- Select and construct the preconditioner
-        autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
-        preconPtr = LduMatrix<Type, DType, LUType>::preconditioner::New
-        (
-            *this,
-            this->controlDict_
-        );
-
-        // --- Solver iteration
-        do
-        {
-            // --- Store previous wArT
-            wArTold = wArT;
-
-            // --- Precondition residuals
-            preconPtr->precondition(wA, rA);
-            preconPtr->preconditionT(wT, rT);
-
-            // --- Update search directions:
-            //wArT = gSumProd(wA, rT);
-            wArT = sum(wA & rT);
-
-            if (solverPerf.nIterations() == 0)
-            {
-                for (register label cell=0; cell<nCells; cell++)
-                {
-                    pAPtr[cell] = wAPtr[cell];
-                    pTPtr[cell] = wTPtr[cell];
-                }
-            }
-            else
-            {
-                scalar beta = cmptDivide(wArT, wArTold);
-
-                for (register label cell=0; cell<nCells; cell++)
-                {
-                    pAPtr[cell] = wAPtr[cell] + (beta* pAPtr[cell]);
-                    pTPtr[cell] = wTPtr[cell] + (beta* pTPtr[cell]);
-                }
-            }
-
-
-            // --- Update preconditioned residuals
-            this->matrix_.Amul(wA, pA);
-            this->matrix_.Tmul(wT, pT);
-
-            scalar wApT = sum(wA & pT);
-
-
-            // --- Test for singularity
-            //if
-            //(
-            //    solverPerf.checkSingularity(ratio(mag(wApT)normFactor))
-            //) break;
-
-
-            // --- Update solution and residual:
-
-            scalar alpha = cmptDivide(wArT, wApT);
-
-            for (register label cell=0; cell<nCells; cell++)
-            {
-                psiPtr[cell] += (alpha* pAPtr[cell]);
-                rAPtr[cell] -= (alpha* wAPtr[cell]);
-                rTPtr[cell] -= (alpha* wTPtr[cell]);
-            }
-
-            solverPerf.finalResidual() =
-                cmptDivide(gSumCmptMag(rA), normFactor);
-
-        } while
-        (
-            solverPerf.nIterations()++ < this->maxIter_
-        && !(solverPerf.converged(this->tolerance_, this->relTol_))
-        );
-    }
-
-    return solverPerf;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCGScalarAlpha/PBiCGScalarAlpha.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCGScalarAlpha/PBiCGScalarAlpha.H
deleted file mode 100644
index 47b467f900afa0c98cb1cd46008270aa57e81c49..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCGScalarAlpha/PBiCGScalarAlpha.H
+++ /dev/null
@@ -1,112 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::TPBiCG
-
-Description
-    Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices
-    using a run-time selectable preconditiioner.
-
-SourceFiles
-    TPBiCG.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef TPBiCG_H
-#define TPBiCG_H
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class TPBiCG Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type, class DType, class LUType>
-class TPBiCG
-:
-    public LduMatrix<Type, DType, LUType>::solver
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        TPBiCG(const TPBiCG&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const TPBiCG&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("PBiCG");
-
-
-    // Constructors
-
-        //- Construct from matrix components and solver data dictionary
-        TPBiCG
-        (
-            const word& fieldName,
-            const LduMatrix<Type, DType, LUType>& matrix,
-            const dictionary& solverDict
-        );
-
-
-    // Destructor
-
-        virtual ~TPBiCG()
-        {}
-
-
-    // Member Functions
-
-        //- Solve the matrix with this solver
-        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
-        (
-            Field<Type>& psi
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "PBiCGScalarAlpha.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCG/TPCG.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCG/TPCG.C
deleted file mode 100644
index c39aef1fc62e274db4b4a2daf5b5d70856e081d2..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCG/TPCG.C
+++ /dev/null
@@ -1,181 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "TPCG.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::TPCG<Type, DType, LUType>::TPCG
-(
-    const word& fieldName,
-    const LduMatrix<Type, DType, LUType>& matrix,
-    const dictionary& solverDict
-)
-:
-    LduMatrix<Type, DType, LUType>::solver
-    (
-        fieldName,
-        matrix,
-        solverDict
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
-Foam::TPCG<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
-    (
-        preconditionerName + typeName,
-        this->fieldName_
-    );
-
-    register label nCells = psi.size();
-
-    Type* __restrict__ psiPtr = psi.begin();
-
-    Field<Type> pA(nCells);
-    Type* __restrict__ pAPtr = pA.begin();
-
-    Field<Type> wA(nCells);
-    Type* __restrict__ wAPtr = wA.begin();
-
-    Type wArA = this->matrix_.great_*pTraits<Type>::one;
-    Type wArAold = wArA;
-
-    // --- Calculate A.psi
-    this->matrix_.Amul(wA, psi);
-
-    // --- Calculate initial residual field
-    Field<Type> rA(this->matrix_.source() - wA);
-    Type* __restrict__ rAPtr = rA.begin();
-
-    // --- Calculate normalisation factor
-    Type normFactor = this->normFactor(psi, wA, pA);
-
-    if (LduMatrix<Type, DType, LUType>::debug >= 2)
-    {
-        Info<< "   Normalisation factor = " << normFactor << endl;
-    }
-
-    // --- Calculate normalised residual norm
-    solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
-    solverPerf.finalResidual() = solverPerf.initialResidual();
-
-    // --- Check convergence, solve if not converged
-    if (!solverPerf.converged(this->tolerance_, this->relTol_))
-    {
-        // --- Select and construct the preconditioner
-        autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
-        preconPtr = LduMatrix<Type, DType, LUType>::preconditioner::New
-        (
-            *this,
-            this->controlDict_
-        );
-
-        // --- Solver iteration
-        do
-        {
-            // --- Store previous wArA
-            wArAold = wArA;
-
-            // --- Precondition residual
-            preconPtr->precondition(wA, rA);
-
-            // --- Update search directions:
-            wArA = gSumCmptProd(wA, rA);
-
-            if (solverPerf.nIterations() == 0)
-            {
-                for (register label cell=0; cell<nCells; cell++)
-                {
-                    pAPtr[cell] = wAPtr[cell];
-                }
-            }
-            else
-            {
-                Type beta = cmptDivide
-                (
-                    wArA,
-                    stabilise(wArAold, this->matrix_.vsmall_)
-                );
-
-                for (register label cell=0; cell<nCells; cell++)
-                {
-                    pAPtr[cell] = wAPtr[cell] + cmptMultiply(beta, pAPtr[cell]);
-                }
-            }
-
-
-            // --- Update preconditioned residual
-            this->matrix_.Amul(wA, pA);
-
-            Type wApA = gSumCmptProd(wA, pA);
-
-
-            // --- Test for singularity
-            if (solverPerf.singular(cmptDivide(cmptMag(wApA), normFactor)))
-            {
-                break;
-            }
-
-
-            // --- Update solution and residual:
-
-            Type alpha = cmptDivide
-            (
-                wArA,
-                stabilise(wApA, this->matrix_.vsmall_)
-            );
-
-            for (register label cell=0; cell<nCells; cell++)
-            {
-                psiPtr[cell] += cmptMultiply(alpha, pAPtr[cell]);
-                rAPtr[cell] -= cmptMultiply(alpha, wAPtr[cell]);
-            }
-
-            solverPerf.finalResidual() =
-                cmptDivide(gSumCmptMag(rA), normFactor);
-
-        } while
-        (
-            solverPerf.nIterations()++ < this->maxIter_
-        && !(solverPerf.converged(this->tolerance_, this->relTol_))
-        );
-    }
-
-    return solverPerf;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCG/TPCG.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/PCG/TPCG.H
deleted file mode 100644
index f92f1ebe59b71423742a1860130e2467da08a725..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/PCG/TPCG.H
+++ /dev/null
@@ -1,112 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::TPCG
-
-Description
-    Preconditioned conjugate gradient solver for symmetric lduMatrices
-    using a run-time selectable preconditiioner.
-
-SourceFiles
-    TPCG.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef TPCG_H
-#define TPCG_H
-
-#include "LduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class TPCG Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type, class DType, class LUType>
-class TPCG
-:
-    public LduMatrix<Type, DType, LUType>::solver
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        TPCG(const TPCG&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const TPCG&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("PCG");
-
-
-    // Constructors
-
-        //- Construct from matrix components and solver data dictionary
-        TPCG
-        (
-            const word& fieldName,
-            const LduMatrix<Type, DType, LUType>& matrix,
-            const dictionary& solverDict
-        );
-
-
-    // Destructor
-
-        virtual ~TPCG()
-        {}
-
-
-    // Member Functions
-
-        //- Solve the matrix with this solver
-        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
-        (
-            Field<Type>& psi
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "TPCG.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C
deleted file mode 100644
index 5a16706f295e26fcb6b9d2ff256036a01cd4e1fa..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C
+++ /dev/null
@@ -1,154 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "SmoothSolver.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-Foam::SmoothSolver<Type, DType, LUType>::SmoothSolver
-(
-    const word& fieldName,
-    const LduMatrix<Type, DType, LUType>& matrix,
-    const dictionary& solverDict
-)
-:
-    LduMatrix<Type, DType, LUType>::solver
-    (
-        fieldName,
-        matrix,
-        solverDict
-    ),
-    nSweeps_(1)
-{
-    readControls();
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type, class DType, class LUType>
-void Foam::SmoothSolver<Type, DType, LUType>::readControls()
-{
-    LduMatrix<Type, DType, LUType>::solver::readControls();
-    readControl(this->controlDict_, nSweeps_, "nSweeps");
-}
-
-
-template<class Type, class DType, class LUType>
-typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
-Foam::SmoothSolver<Type, DType, LUType>::solve(Field<Type>& psi) const
-{
-    // --- Setup class containing solver performance data
-    typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
-    (
-        typeName,
-        this->fieldName_
-    );
-
-    // If the nSweeps_ is negative do a fixed number of sweeps
-    if (nSweeps_ < 0)
-    {
-        autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
-        smootherPtr = LduMatrix<Type, DType, LUType>::smoother::New
-        (
-            this->fieldName_,
-            this->matrix_,
-            this->controlDict_
-        );
-
-        smootherPtr->smooth(psi, -nSweeps_);
-
-        solverPerf.nIterations() -= nSweeps_;
-    }
-    else
-    {
-        Type normFactor = pTraits<Type>::zero;
-
-        {
-            Field<Type> Apsi(psi.size());
-            Field<Type> temp(psi.size());
-
-            // Calculate A.psi
-            this->matrix_.Amul(Apsi, psi);
-
-            // Calculate normalisation factor
-            normFactor = this->normFactor(psi, Apsi, temp);
-
-            // Calculate residual magnitude
-            solverPerf.initialResidual() = cmptDivide
-            (
-                gSumCmptMag(this->matrix_.source() - Apsi),
-                normFactor
-            );
-            solverPerf.finalResidual() = solverPerf.initialResidual();
-        }
-
-        if (LduMatrix<Type, DType, LUType>::debug >= 2)
-        {
-            Info<< "   Normalisation factor = " << normFactor << endl;
-        }
-
-
-        // Check convergence, solve if not converged
-        if (!solverPerf.converged(this->tolerance_, this->relTol_))
-        {
-            autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
-            smootherPtr = LduMatrix<Type, DType, LUType>::smoother::New
-            (
-                this->fieldName_,
-                this->matrix_,
-                this->controlDict_
-            );
-
-            // Smoothing loop
-            do
-            {
-                smootherPtr->smooth
-                (
-                    psi,
-                    nSweeps_
-                );
-
-                // Calculate the residual to check convergence
-                solverPerf.finalResidual() = cmptDivide
-                (
-                    gSumCmptMag(this->matrix_.residual(psi)),
-                    normFactor
-                );
-            } while
-            (
-                (solverPerf.nIterations() += nSweeps_) < this->maxIter_
-             && !(solverPerf.converged(this->tolerance_, this->relTol_))
-            );
-        }
-    }
-
-    return solverPerf;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H b/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H
deleted file mode 100644
index e7baec3b26398fcfdb4e696d5a556716786ad16c..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H
+++ /dev/null
@@ -1,111 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::SmoothSolver
-
-Description
-    Iterative solver for symmetric and assymetric matrices which uses a
-    run-time selected smoother e.g. GaussSeidel to converge the solution to
-    the required tolerance.  To improve efficiency, the residual is evaluated 
-    after every nSweeps smoothing iterations.
-
-SourceFiles
-    SmoothSolver.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef SmoothSolver_H
-#define SmoothSolver_H
-
-#include "lduMatrix.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class SmoothSolver Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type, class DType, class LUType>
-class SmoothSolver
-:
-    public LduMatrix<Type, DType, LUType>::solver
-{
-
-protected:
-
-    // Protected data
-
-        //- Number of sweeps before the evaluation of residual
-        label nSweeps_;
-
-        //- Read the control parameters from the controlDict_
-        virtual void readControls();
-
-
-public:
-
-    //- Runtime type information
-    TypeName("SmoothSolver");
-
-
-    // Constructors
-
-        //- Construct from matrix components and solver data dictionary
-        SmoothSolver
-        (
-            const word& fieldName,
-            const LduMatrix<Type, DType, LUType>& matrix,
-            const dictionary& solverDict
-        );
-
-
-    // Member Functions
-
-        //- Solve the matrix with this solver
-        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
-        (
-            Field<Type>& psi
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "SmoothSolver.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LduMatrix/Solvers/lduSolvers.C b/src/OpenFOAM/matrices/LduMatrix/Solvers/lduSolvers.C
deleted file mode 100644
index 4b9cbfac4e0dc0a9b92b885684cc31edabdbca35..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/LduMatrix/Solvers/lduSolvers.C
+++ /dev/null
@@ -1,29 +0,0 @@
-#include "TPCG.H"
-#include "TPBiCG.H"
-#include "SmoothSolver.H"
-#include "fieldTypes.H"
-
-#define makeLduSolvers(Type, DType, LUType)                                   \
-                                                                              \
-    makeLduSolver(DiagonalSolver, Type, DType, LUType);                       \
-    makeLduSymSolver(DiagonalSolver, Type, DType, LUType);                    \
-    makeLduAsymSolver(DiagonalSolver, Type, DType, LUType);                   \
-                                                                              \
-    makeLduSolver(TPCG, Type, DType, LUType);                                 \
-    makeLduSymSolver(TPCG, Type, DType, LUType);                              \
-                                                                              \
-    makeLduSolver(TPBiCG, Type, DType, LUType);                               \
-    makeLduAsymSolver(TPBiCG, Type, DType, LUType);                           \
-                                                                              \
-    makeLduSolver(SmoothSolver, Type, DType, LUType);                         \
-    makeLduSymSolver(SmoothSolver, Type, DType, LUType);                      \
-    makeLduAsymSolver(SmoothSolver, Type, DType, LUType);
-
-namespace Foam
-{
-    makeLduSolvers(scalar, scalar, scalar);
-    makeLduSolvers(vector, scalar, scalar);
-    makeLduSolvers(sphericalTensor, scalar, scalar);
-    makeLduSolvers(symmTensor, scalar, scalar);
-    makeLduSolvers(tensor, scalar, scalar);
-};