Added missing dirs

src/OpenFOAM/matrices/Matrix/Matrix.C

@@ -0,0 +1,407 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 "Matrix.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template<class T>
+void Foam::Matrix<T>::allocate()
+{
+    if (n_ && m_)
+    {
+        v_ = new T*[n_];
+        v_[0] = new T[n_*m_];
+
+        for (register label i=1; i<n_; i++)
+        {
+            v_[i] = v_[i-1] + m_;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
+
+template<class T>
+Foam::Matrix<T>::~Matrix()
+{
+    if (v_)
+    {
+        delete[] (v_[0]);
+        delete[] v_;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class T>
+const Foam::Matrix<T>& Foam::Matrix<T>::null()
+{
+    Matrix<T>* nullPtr = reinterpret_cast<Matrix<T>*>(NULL);
+    return *nullPtr;
+}
+
+
+template<class T>
+Foam::Matrix<T>::Matrix(const label n, const label m)
+:
+    v_(NULL),
+    n_(n),
+    m_(m)
+{
+    if (n_ < 0 || m_ < 0)
+    {
+        FatalErrorIn("Matrix<T>::Matrix(const label n, const label m)")
+            << "bad n, m " << n_ << ", " << m_
+            << abort(FatalError);
+    }
+
+    allocate();
+}
+
+
+template<class T>
+Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
+:
+    v_(NULL),
+    n_(n),
+    m_(m)
+{
+    if (n_ < 0 || m_ < 0)
+    {
+        FatalErrorIn
+        (
+            "Matrix<T>::Matrix(const label n, const label m, const T&)"
+        )   << "bad n, m " << n_ << ", " << m_
+            << abort(FatalError);
+    }
+
+    allocate();
+
+    if (v_)
+    {
+        T* v = v_[0];
+
+        label nm = n_*m_;
+
+        for (register label i=0; i<nm; i++)
+        {
+            v[i] = a;
+        }
+    }
+}
+
+
+template<class T>
+Foam::Matrix<T>::Matrix(const Matrix<T>& a)
+:
+    v_(NULL),
+    n_(a.n_),
+    m_(a.m_)
+{
+    if (a.v_)
+    {
+        allocate();
+        T* v = v_[0];
+        const T* av = a.v_[0];
+
+        label nm = n_*m_;
+        for (register label i=0; i<nm; i++)
+        {
+            v[i] = av[i];
+        }
+    }
+}
+
+        
+template<class T>
+void Foam::Matrix<T>::clear()
+{
+    if (v_)
+    {
+        delete[] (v_[0]);
+        delete[] v_;
+    }
+    n_ = 0;
+    m_ = 0;
+    v_ = NULL;
+}
+
+
+template<class T>
+void Foam::Matrix<T>::transfer(Matrix<T>& a)
+{
+    clear();
+
+    n_ = a.n_;
+    a.n_ = 0;
+
+    m_ = a.m_;
+    a.m_ = 0;
+
+    v_ = a.v_;
+    a.v_ = NULL;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class T>
+void Foam::Matrix<T>::operator=(const T& t)
+{
+    if (v_)
+    {
+        T* v = v_[0];
+
+        label nm = n_*m_;
+        for (register label i=0; i<nm; i++)
+        {
+            v[i] = t;
+        }
+    }
+}
+
+
+// Assignment operator. Takes linear time.
+template<class T>
+void Foam::Matrix<T>::operator=(const Matrix<T>& a)
+{
+    if (this == &a)
+    {
+        FatalErrorIn("Matrix<T>::operator=(const Matrix<T>&)")
+            << "attempted assignment to self"
+            << abort(FatalError);
+    }
+
+    if (n_ != a.n_ || m_ != a.m_)
+    {
+        clear();
+        n_ = a.n_;
+        m_ = a.m_;
+        allocate();
+    } 
+
+    if (v_)
+    {
+        T* v = v_[0];
+        const T* av = a.v_[0];
+
+        label nm = n_*m_;
+        for (register label i=0; i<nm; i++)
+        {
+            v[i] = av[i];
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+template<class T>
+const T& Foam::max(const Matrix<T>& a)
+{
+    label nm = a.n_*a.m_;
+
+    if (nm)
+    {
+        label curMaxI = 0;
+        const T* v = a.v_[0];
+
+        for (register label i=1; i<nm; i++)
+        {
+            if (v[i] > v[curMaxI])
+            {
+                curMaxI = i;
+            }
+        }
+
+        return v[curMaxI];
+    }
+    else
+    {
+        FatalErrorIn("max(const Matrix<T>&)")
+            << "matrix is empty"
+            << abort(FatalError);
+
+        // Return in error to keep compiler happy
+        return a[0][0];
+    }
+}
+
+
+template<class T>
+const T& Foam::min(const Matrix<T>& a)
+{
+    label nm = a.n_*a.m_;
+
+    if (nm)
+    {
+        label curMinI = 0;
+        const T* v = a.v_[0];
+
+        for (register label i=1; i<nm; i++)
+        {
+            if (v[i] < v[curMinI])
+            {
+                curMinI = i;
+            }
+        }
+
+        return v[curMinI];
+    }
+    else
+    {
+        FatalErrorIn("min(const Matrix<T>&)")
+            << "matrix is empty"
+            << abort(FatalError);
+
+        // Return in error to keep compiler happy
+        return a[0][0];
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
+
+template<class T>
+Foam::Matrix<T> Foam::operator-(const Matrix<T>& a)
+{
+    Matrix<T> na(a.n_, a.m_);
+
+    if (a.n_ && a.m_)
+    {
+        T* nav = na.v_[0];
+        const T* av = a.v_[0];
+
+        label nm = a.n_*a.m_;
+        for (register label i=0; i<nm; i++)
+        {
+            nav[i] = -av[i];
+        }
+    }
+
+    return na;
+}
+
+
+template<class T>
+Foam::Matrix<T> Foam::operator+(const Matrix<T>& a, const Matrix<T>& b)
+{
+    if (a.n_ != b.n_)
+    {
+        FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
+            << "attempted add matrices with different number of rows: "
+            << a.n_ << ", " << b.n_
+            << abort(FatalError);
+    }
+
+    if (a.m_ != b.m_)
+    {
+        FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
+            << "attempted add matrices with different number of columns: "
+            << a.m_ << ", " << b.m_
+            << abort(FatalError);
+    }
+
+    Matrix<T> ab(a.n_, a.m_);
+
+    T* abv = ab.v_[0];
+    const T* av = a.v_[0];
+    const T* bv = b.v_[0];
+
+    label nm = a.n_*a.m_;;
+    for (register label i=0; i<nm; i++)
+    {
+        abv[i] = av[i] + bv[i];
+    }
+
+    return ab;
+}
+
+
+template<class T>
+Foam::Matrix<T> Foam::operator-(const Matrix<T>& a, const Matrix<T>& b)
+{
+    if (a.n_ != b.n_)
+    {
+        FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
+            << "attempted add matrices with different number of rows: "
+            << a.n_ << ", " << b.n_
+            << abort(FatalError);
+    }
+
+    if (a.m_ != b.m_)
+    {
+        FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
+            << "attempted add matrices with different number of columns: "
+            << a.m_ << ", " << b.m_
+            << abort(FatalError);
+    }
+
+    Matrix<T> ab(a.n_, a.m_);
+
+    T* abv = ab.v_[0];
+    const T* av = a.v_[0];
+    const T* bv = b.v_[0];
+
+    label nm = a.n_*a.m_;;
+    for (register label i=0; i<nm; i++)
+    {
+        abv[i] = av[i] - bv[i];
+    }
+
+    return ab;
+}
+
+
+template<class T>
+Foam::Matrix<T> Foam::operator*(const scalar s, const Matrix<T>& a)
+{
+    Matrix<T> sa(a.n_, a.m_);
+
+    if (a.n_ && a.m_)
+    {
+        T* sav = sa.v_[0];
+        const T* av = a.v_[0];
+
+        label nm = a.n_*a.m_;;
+        for (register label i=0; i<nm; i++)
+        {
+            sav[i] = s*av[i];
+        }
+    }
+
+    return sa;
+}
+
+
+// * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
+
+#include "MatrixIO.C"
+
+// ************************************************************************* //

src/OpenFOAM/matrices/Matrix/Matrix.H

@@ -0,0 +1,211 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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::Matrix
+
+Description
+    A templated 2D matrix of objects of \<T\>, where the n x m matrix
+    dimensions are known and used for subscript bounds checking, etc.
+
+SourceFiles
+    Matrix.C
+    MatrixI.H
+    MatrixIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Matrix_H
+#define Matrix_H
+
+#include "List.H"
+#include "label.H"
+#include "bool.H"
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class T> class Matrix;
+
+template<class T> const T& max(const Matrix<T>&);
+template<class T> const T& min(const Matrix<T>&);
+
+template<class T> Matrix<T> operator-(const Matrix<T>&);
+template<class T> Matrix<T> operator+(const Matrix<T>&, const Matrix<T>&);
+template<class T> Matrix<T> operator-(const Matrix<T>&, const Matrix<T>&);
+template<class T> Matrix<T> operator*(const scalar, const Matrix<T>&);
+
+template<class T> Istream& operator>>(Istream&, Matrix<T>&);
+template<class T> Ostream& operator<<(Ostream&, const Matrix<T>&);
+
+
+/*---------------------------------------------------------------------------*\
+                           Class Matrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class T>
+class Matrix
+{
+    // Private data
+
+        //- Row pointers
+        T** __restrict__ v_;
+
+        //- Number of rows and columns in Matrix.
+        label n_, m_;
+
+        //- Allocate the storage for the row-pointers and the data
+        //  and set the row pointers
+        void allocate();
+
+
+public:
+
+    // Constructors
+
+        //- Null constructor.
+        inline Matrix();
+
+        //- Construct given number of rows and columns.
+        Matrix(const label n, const label m);
+
+        //- Construct with given number of rows and columns
+        //  and value for all elements.
+        Matrix(const label n, const label m, const T&);
+
+        //- Copy constructor.
+        Matrix(const Matrix<T>&);
+
+        //- Construct from Istream.
+        Matrix(Istream&);
+
+        //- Clone
+        inline autoPtr<Matrix<T> > clone() const;
+
+
+    // Destructor
+
+        ~Matrix();
+
+
+    // Member functions
+
+        //- Return a null Matrix
+        static const Matrix<T>& null();
+
+
+        // Access
+
+            //- Return the number of rows
+            inline label n() const;
+
+            //- Return the number of columns
+            inline label m() const;
+
+            //- Return the number of elements in matrix (n*m)
+            inline label size() const;
+
+
+        // Check
+
+            //- Check index i is within valid range (0 ... n-1).
+            inline void checki(const label i) const;
+
+            //- Check index j is within valid range (0 ... m-1).
+            inline void checkj(const label j) const;
+
+
+        // Edit
+
+            //- Clear the Matrix, i.e. set sizes to zero.
+            void clear();
+
+            //- Transfer the contents of the argument Matrix into this Matrix
+            //  and annull the argument Matrix.
+            void transfer(Matrix<T>&);
+
+
+    // Member operators
+
+        //- Return subscript-checked element of Matrix.
+        inline T* operator[](const label);
+
+        //- Return subscript-checked element of constant Matrix.
+        inline const T* operator[](const label) const;
+
+        //- Assignment operator. Takes linear time.
+        void operator=(const Matrix<T>&);
+
+        //- Assignment of all entries to the given value
+        void operator=(const T&);
+
+
+    // Friend functions
+
+        friend const T& max<T>(const Matrix<T>&);
+        friend const T& min<T>(const Matrix<T>&);
+
+
+    // Friend operators
+
+        friend Matrix<T> operator-<T>(const Matrix<T>&);
+        friend Matrix<T> operator+<T>(const Matrix<T>&, const Matrix<T>&);
+        friend Matrix<T> operator-<T>(const Matrix<T>&, const Matrix<T>&);
+        friend Matrix<T> operator*<T>(const scalar, const Matrix<T>&);
+
+
+    // IOstream operators
+
+        //- Read Matrix from Istream, discarding contents of existing Matrix.
+        friend Istream& operator>> <T>(Istream&, Matrix<T>&);
+
+        // Write Matrix to Ostream.
+        friend Ostream& operator<< <T>(Ostream&, const Matrix<T>&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "MatrixI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Matrix.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/Matrix/MatrixI.H

@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class T>
+inline Foam::Matrix<T>::Matrix()
+:
+    v_(NULL),
+    n_(0),
+    m_(0)
+{}
+
+
+template<class T>
+inline Foam::autoPtr<Foam::Matrix<T> > Foam::Matrix<T>::clone() const
+{
+    return autoPtr<Matrix<T> >(new Matrix<T>(*this));
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//- Return the number of rows
+template<class T>
+inline Foam::label Foam::Matrix<T>::n() const
+{
+    return n_;
+}
+
+
+//- Return the number of columns
+template<class T>
+inline Foam::label Foam::Matrix<T>::m() const
+{
+    return m_;
+}
+
+
+//- Return the number of columns
+template<class T>
+inline Foam::label Foam::Matrix<T>::size() const
+{
+    return n_*m_;
+}
+
+
+// Check index i is within valid range (0 ... n-1).
+template<class T>
+inline void Foam::Matrix<T>::checki(const label i) const
+{
+    if (!n_)
+    {
+        FatalErrorIn("Matrix<T>::checki(const label)")
+            << "attempt to access element from zero sized row"
+            << abort(FatalError);
+    }
+    else if (i<0 || i>=n_)
+    {
+        FatalErrorIn("Matrix<T>::checki(const label)")
+            << "index " << i << " out of range 0 ... " << n_-1
+            << abort(FatalError);
+    }
+}
+
+
+// Check index j is within valid range (0 ... n-1).
+template<class T>
+inline void Foam::Matrix<T>::checkj(const label j) const
+{
+    if (!m_)
+    {
+        FatalErrorIn("Matrix<T>::checkj(const label)")
+            << "attempt to access element from zero sized column"
+            << abort(FatalError);
+    }
+    else if (j<0 || j>=m_)
+    {
+        FatalErrorIn("Matrix<T>::checkj(const label)")
+            << "index " << j << " out of range 0 ... " << m_-1
+            << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+// Return subscript-checked element access
+template<class T>
+inline T* Foam::Matrix<T>::operator[](const label i)
+{
+#   ifdef FULLDEBUG
+    checki(i);
+#   endif
+    return v_[i];
+}
+
+
+// Return subscript-checked const element access
+template<class T>
+inline const T* Foam::Matrix<T>::operator[](const label i) const
+{
+#   ifdef FULLDEBUG
+    checki(i);
+#   endif
+    return v_[i];
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/Matrix/MatrixIO.C

@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 "Matrix.H"
+#include "Istream.H"
+#include "Ostream.H"
+#include "token.H"
+#include "contiguous.H"
+
+// * * * * * * * * * * * * * * * Ostream Operator *  * * * * * * * * * * * * //
+
+template<class T>
+Foam::Matrix<T>::Matrix(Istream& is)
+:
+    v_(NULL),
+    n_(0),
+    m_(0)
+{
+    operator>>(is, *this);
+}
+
+
+template<class T>
+Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
+{
+    // Anull matrix
+    M.clear();
+
+    is.fatalCheck("operator>>(Istream&, Matrix<T>&)");
+
+    token firstToken(is);
+
+    is.fatalCheck("operator>>(Istream&, Matrix<T>&) : reading first token");
+
+    if (firstToken.isLabel())
+    {
+        M.n_ = firstToken.labelToken();
+        M.m_ = readLabel(is);
+
+        label nm = M.n_*M.m_;
+
+        // Read list contents depending on data format
+        if (is.format() == IOstream::ASCII || !contiguous<T>())
+        {
+            // Read beginning of contents
+            char listDelimiter = is.readBeginList("Matrix");
+
+            if (nm)
+            {
+                M.allocate();
+                T* v = M.v_[0];
+
+                if (listDelimiter == token::BEGIN_LIST)
+                {
+                    label k = 0;
+
+                    // loop over rows
+                    for (register label i=0; i<M.n(); i++)
+                    {
+                        listDelimiter = is.readBeginList("MatrixRow");
+
+                        for (register label j=0; j<M.m(); j++)
+                        {
+                            is >> v[k++];
+
+                            is.fatalCheck
+                            (
+                                "operator>>(Istream&, Matrix<T>&) : "
+                                "reading entry"
+                            );
+                        }
+
+                        is.readEndList("MatrixRow");
+                    }
+                }
+                else
+                {
+                    T element;
+                    is >> element;
+
+                    is.fatalCheck
+                    (
+                        "operator>>(Istream&, Matrix<T>&) : "
+                        "reading the single entry"
+                    );
+
+                    for (register label i=0; i<nm; i++)
+                    {
+                        v[i] = element;
+                    }
+                }
+            }
+
+            // Read end of contents
+            is.readEndList("Matrix");
+        }
+        else
+        {
+            if (nm)
+            {
+                M.allocate();
+                T* v = M.v_[0];
+
+                is.read(reinterpret_cast<char*>(v), nm*sizeof(T));
+
+                is.fatalCheck
+                (
+                    "operator>>(Istream&, Matrix<T>&) : "
+                    "reading the binary block"
+                );
+            }
+        }
+    }
+    else
+    {
+        FatalIOErrorIn("operator>>(Istream&, Matrix<T>&)", is)
+            << "incorrect first token, expected <int>, found "
+            << firstToken.info()
+            << exit(FatalIOError);
+    }
+
+    return is;
+}
+
+
+template<class T>
+Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<T>& M)
+{
+    label nm = M.n_*M.m_;
+
+    os  << M.n() << token::SPACE << M.m();
+
+    // Write list contents depending on data format
+    if (os.format() == IOstream::ASCII || !contiguous<T>())
+    {
+        if (nm)
+        {
+            bool uniform = false;
+
+            const T* v = M.v_[0];
+
+            if (nm > 1 && contiguous<T>())
+            {
+                uniform = true;
+
+                for (register label i=0; i< nm; i++)
+                {
+                    if (v[i] != v[0])
+                    {
+                        uniform = false;
+                        break;
+                    }
+                }
+            }
+
+            if (uniform)
+            {
+                // Write size of list and start contents delimiter
+                os  << token::BEGIN_BLOCK;
+
+                // Write list contents
+                os << v[0];
+
+                // Write end of contents delimiter
+                os << token::END_BLOCK;
+            }
+            // Fix: matrices smaller than 20x20 will be written square.
+            // HJ, 22/Jan/2009
+            else if (nm < 400 && contiguous<T>())
+            {
+                // Write size of list and start contents delimiter
+                os  << token::BEGIN_LIST;
+
+                label k = 0;
+
+                // loop over rows
+                for (register label i=0; i< M.n(); i++)
+                {
+                    os  << nl << token::BEGIN_LIST;
+
+                    // Write row
+                    for (register label j=0; j< M.m(); j++)
+                    {
+                        if (j > 0) os << token::SPACE;
+                        os << v[k++];
+                    }
+
+                    os << token::END_LIST;
+                }
+
+                // Write end of contents delimiter
+                os << token::END_LIST;
+            }
+            else
+            {
+                // Write size of list and start contents delimiter
+                os  << nl << token::BEGIN_LIST;
+
+                label k = 0;
+
+                // loop over rows
+                for (register label i=0; i< M.n(); i++)
+                {
+                    os  << nl << token::BEGIN_LIST;
+
+                    // Write row
+                    for (register label j=0; j< M.m(); j++)
+                    {
+                        os << nl << v[k++];
+                    }
+
+                    os << nl << token::END_LIST;
+                }
+
+                // Write end of contents delimiter
+                os << nl << token::END_LIST << nl;
+            }
+        }
+        else
+        {
+            os  << token::BEGIN_LIST << token::END_LIST << nl;
+        }
+    }
+    else
+    {
+        if (nm)
+        {
+            os.write(reinterpret_cast<const char*>(M.v_[0]), nm*sizeof(T));
+        }
+    }
+
+    // Check state of IOstream
+    os.check("Ostream& operator<<(Ostream&, const Matrix&)");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/Matrix/tools/DenseMatrixTools.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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
+    DenseMatrixTools
+
+\*----------------------------------------------------------------------------*/
+
+#include "DenseMatrixTools.H"
+#include "Swap.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Solve with pivoting.  Note: Destroys matrix and b
+template<class T>
+static void Foam::DenseMatrixTools::solve
+(
+    Matrix<T>& A,
+    List<T>& x,
+    List<T>& b
+)
+{
+    const label nRows = A.n();
+
+    // Create and initialise permutation array
+    labelList p(nRows);
+
+    forAll (p, i)
+    {
+        p[i] = i;
+    }
+
+
+    // Eliminate equations
+    for (label k = 0; k < nRows - 1; k++)
+    {
+        // Swap rows
+        label m = k;
+
+        for(label i = k + 1; i < nRows; ++i)
+        {
+            if (Foam::mag(A[p[i]][k]) > Foam::mag(A[p[m]][k]))
+            {
+                m = i;
+            }
+
+            Swap(p[k], p[m]);
+        }
+
+
+        for (label i = k + 1; i < nRows; i++)
+        {
+            label pi = p[i];
+            label pk = p[k];
+
+//             T r = A[pi][k]/A[pk][k];
+            T r = cmptDivide(A[pi][k], A[pk][k]);
+
+            for (label j = k + 1; j < nRows; j++)
+            {
+//                 A[pi][j] -= A[pk][j]*r;
+                A[pi][j] -= cmptMultiply(A[pk][j], r);
+            }
+
+//             b[pi] -= b[pk]*r;
+            b[pi] -= cmptMultiply(b[pk], r);
+        }
+    }
+
+    // Back substitute
+    for (label i = nRows - 1; i >= 0; i--)
+    {
+        label pi = p[i];
+
+        T sum = b[pi];
+
+        for (label j = i + 1; j < nRows; j++)
+        {
+//             sum -= A[pi][j]*x[j];
+            sum -= cmptMultiply(A[pi][j], x[j]);
+        }
+
+//         x[i] = sum/A[pi][i];
+        x[i] = cmptDivide(sum, A[pi][i]);
+    }
+}
+
+
+template<class T>
+static void Foam::DenseMatrixTools::qrDecompose
+(
+    const label nCols,
+    FieldField<Field, T>& A,
+    Matrix<T>& R
+)
+{
+    // Note: consider Arnoldi algorithm for speed-up.  HJ, 14/Sep/2006
+
+    for (label j = 0; j < nCols; j++)
+    {
+        R[j][j] = Foam::sqrt(gSumSqr(A[j]));
+
+        if (mag(R[j][j]) > SMALL)
+        {
+            A[j] /= R[j][j];
+        }
+
+        for (label k = j + 1; k < nCols; k++)
+        {
+            const Field<T>& Aj = A[j];
+            Field<T>& Ak = A[k];
+            T& Rjk = R[j][k];
+
+            Rjk = gSumProd(Aj, Ak);
+
+            for (label i = 0; i < Ak.size(); i++)
+            {
+                Ak[i] -= Rjk*Aj[i];
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/Matrix/tools/DenseMatrixTools.H

@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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
+
+Namespace
+    DenseMatrixTools
+
+SourceFiles
+    DenseMatrixTools.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DenseMatrixTools_H
+#define DenseMatrixTools_H
+
+#include "Matrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class DenseMatrixTools Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace DenseMatrixTools
+{
+    //- Solve using Gauss Elimination with pivoting.  Destroys matrix and b
+    template<class T>
+    static void solve(Matrix<T>& A, List<T>& x, List<T>& b);
+
+    //- Q-R decomposition
+    template<class T>
+    static void qrDecompose
+    (
+        const label nCols,
+        FieldField<Field, T>& A,
+        Matrix<T>& R
+    );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "DenseMatrixTools.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/Matrix/tools/EigenSolver/EigenSolver.C

@@ -0,0 +1,297 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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
+    EigenSolver
+
+Description
+
+\*----------------------------------------------------------------------------*/
+
+#include "EigenSolver.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class T>
+void Foam::EigenSolver<T>::checkMatrix(const Matrix<T>& mtx) const
+{
+    // Check shape
+    if (mtx.m() != mtx.n())
+    {
+        FatalErrorIn
+        (
+            "void EigenSolver<T>::checkMatrix(const Matrix<T>& m) const"
+        )   << "Matrix is not square, cannot calculate eigen-values.  "
+            << "M = " << mtx.m() << " N = " << mtx.n()
+            << abort(FatalError);
+    }
+
+    // Check symmetry
+    T assymetry = pTraits<T>::zero;
+
+    for (label i = 0; i < mtx.m(); i++)
+    {
+        for (label j = i; j < mtx.n(); j++)
+        {
+            assymetry += mag(mtx[i][j] - mtx[j][i]);
+        }
+    }
+
+    if (assymetry > mtx.m()*SMALL)
+    {
+        FatalErrorIn
+        (
+            "void EigenSolver<T>::checkMatrix(const Matrix<T>& m) const"
+        )   << "Matrix is not symmetric.  Assymetry = " << assymetry
+            << abort(FatalError);
+    }
+}
+
+
+template<class T>
+void Foam::EigenSolver<T>::factorise(const Matrix<T>& mtx)
+{
+    // Copy the original matrix into scratch space
+    Matrix<T> a = mtx;
+
+    label n = a.m();
+
+    // Create and initialise the matrix to hold eigen-vectors in columns
+    Matrix<T> v(n, n, pTraits<T>::zero);
+
+    for (label ip = 0; ip < n; ip++)
+    {
+        v[ip][ip] = pTraits<T>::one;
+    }
+
+    // Create and initialise scratch vectors
+    List<T> b(n);
+    List<T>& d = values_;
+
+    // Initialise b and d to the diagonal of a
+    for (label ip = 0; ip < n; ip++)
+    {
+        b[ip] = a[ip][ip];
+        d[ip] = a[ip][ip];
+    }
+
+    List<T> z(n, pTraits<T>::zero);
+
+    label nRotations = 0;
+
+    // Main iteration loop
+    for (label i = 0; i <= 50; i++)
+    {
+        T sm = pTraits<T>::zero;
+        for (label ip = 0; ip < n - 1; ip++)
+        {
+            for (label iq = ip + 1; iq < n; iq++)
+            {
+                sm += mag(a[ip][iq]);
+            }
+        }
+
+        if (sm < VSMALL)
+        {
+            // Normal return, relying on quadratic convergence
+
+            // Copy eigen-vectors.  Note that v stores vectors in columns
+            for (label ip = 0; ip < n; ip++)
+            {
+                for (label iq = 0; iq < n; iq++)
+                {
+                    vectors_[iq][ip] = v[ip][iq];
+                }
+            }
+
+            return;
+        }
+
+        // Calculate treshold
+        T treshold = pTraits<T>::zero;
+
+        if (i < 4)
+        {
+            // Operation on first 3 sweeps
+            treshold = 0.2*sm/(n*n);
+        }
+
+        T theta, tau, g, h, t, c, s;
+
+        for (label ip = 0; ip < n - 1; ip++)
+        {
+            for (label iq = ip + 1; iq < n; iq++)
+            {
+                g = 100.0*mag(a[ip][iq]);
+
+                // After four sweeps, skip the rotation if the off-diagonal
+                // element is small
+                if
+                (
+                    i > 4
+                 && mag(d[ip] + g) == mag(d[ip])
+                 && mag(d[iq] + g) == mag(d[iq])
+                )
+                {
+                    a[ip][iq] = pTraits<T>::zero;
+                }
+                else if (mag(a[ip][iq]) > treshold)
+                {
+                    h = d[iq] - d[ip];
+
+                    if (mag(h) + g == mag(h))
+                    {
+                        t = a[ip][iq]/h;
+                    }
+                    else
+                    {
+                        theta = 0.5*h/(a[ip][iq]);
+
+                        t = 1.0/(mag(theta) + Foam::sqrt(1.0 + sqr(theta)));
+
+                        if (theta < 0)
+                        {
+                            t = -t;
+                        }
+                    }
+
+                    c = 1.0/sqrt(1.0 + sqr(t));
+                    s = t*c;
+                    tau = s/(1.0 + c);
+                    h = t*a[ip][iq];
+
+                    z[ip] -= h;
+                    z[iq] += h;
+                    d[ip] -= h;
+                    d[iq] += h;
+
+                    a[ip][iq] = pTraits<T>::zero;
+
+                    // Rotations 0 <= j < p
+                    for (label j = 0; j < ip; j++)
+                    {
+                        rotate(a, s, tau, j, ip, j, iq);
+                    }
+
+                    // Rotations p < j < q
+                    for (label j = ip + 1; j < iq; j++)
+                    {
+                        rotate(a, s, tau, ip, j, j, iq);
+                    }
+
+                    // Rotations q < j < n
+                    for (label j = iq + 1; j < n; j++)
+                    {
+                        rotate(a, s, tau, ip, j, iq, j);
+                    }
+
+                    for (label j = 0; j < n; j++)
+                    {
+                        rotate(v, s, tau, j, ip, j, iq);
+                    }
+
+                    // Increment number of rotations
+                    nRotations++;
+                }
+            }
+        }
+
+        // Update d with the sum of t a_pq and clear z
+        for (label ip = 0; ip < n; ip++)
+        {
+            b[ip] += z[ip];
+            d[ip] = b[ip];
+            z[ip] = pTraits<T>::zero;
+        }
+
+    } // End of main iteration loop
+
+    FatalErrorIn("void Foam::EigenSolver<T>::factorise(const Matrix<T>& mtx)")
+        << "Maximum number of iterations exceeded"
+        << abort(FatalError);
+}
+
+
+template<class T>
+inline void Foam::EigenSolver<T>::rotate
+(
+    Matrix<T>& a,
+    const T s,
+    const T tau,
+    const label i,
+    const label j, 
+    const label k, 
+    const label l
+) const
+{
+    T g = a[i][j];
+    T h = a[k][l];
+
+    a[i][j] = g - s*(h + g*tau);
+    a[k][l] = h + s*(g - h*tau);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+template<class T>
+Foam::EigenSolver<T>::EigenSolver(const Matrix<T>& mtx)
+:
+    values_(mtx.m(), pTraits<T>::zero),
+    vectors_(mtx.m())
+{
+    // Size the eigen vectors
+    forAll (vectors_, rowI)
+    {
+        vectors_[rowI].setSize(mtx.m());
+        vectors_[rowI] = pTraits<T>::zero;
+    }
+
+    this->checkMatrix(mtx);
+
+    factorise(mtx);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class T>
+T Foam::EigenSolver<T>::eigenValue(const label n) const
+{
+    return values_[n];
+}
+
+
+template<class T>
+const Foam::List<T>& Foam::EigenSolver<T>::eigenVector(const label n) const
+{
+    return vectors_[n];
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/Matrix/tools/EigenSolver/EigenSolver.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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
+    EigenSolver
+
+Description
+    Calculate eigen-values and eigen-vectors of a symmetric dense matrix
+
+SourceFiles
+    EigenSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef EigenSolver_H
+#define EigenSolver_H
+
+#include "Matrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Class forward declarations
+
+/*---------------------------------------------------------------------------*\
+                           Class EigenSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class T>
+class EigenSolver
+{
+    // Private data
+
+        //- Eigenvalues
+        List<T> values_;
+
+        //- Eigenvectors
+        List<List<T> > vectors_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        EigenSolver(const EigenSolver&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const EigenSolver&);
+
+
+        //- Check matrix for shape and symmetry
+        void checkMatrix(const Matrix<T>& mtx) const;
+
+        //- Factorise into eigen-values and eigen-vectors
+        void factorise(const Matrix<T>& mtx);
+
+        //- Rotate the matrix
+        inline void rotate
+        (
+            Matrix<T>& a,
+            const T s,
+            const T tau,
+            const label i,
+            const label j,
+            const label k,
+            const label l
+        ) const;
+
+
+public:
+
+    // Static data members
+
+
+    // Constructors
+
+        //- Construct from matrix
+        EigenSolver(const Matrix<T>& mtx);
+
+
+    // Destructor - default
+
+
+    // Member Functions
+
+        //- Return nth eigen value
+        T eigenValue(const label n) const;
+
+        // Return nth eigen vector
+        const List<T>& eigenVector(const label n) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "EigenSolver.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/check/Make/files

@@ -0,0 +1,4 @@
+ggiCheck/ggiCheckFunctionObject.C
+meshCheck/meshCheckFunctionObject.C
+
+LIB = $(FOAM_LIBBIN)/libcheckFunctionObjects

src/postProcessing/functionObjects/check/Make/options

@@ -0,0 +1,4 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+LIB_LIBS = \
+    -lfiniteVolume

wmake/rules/linuxGcc/mplibSYSTEMOPENMPI

@@ -0,0 +1 @@
+PFLAGS     = -DOMPI_SKIP_MPICXX