Merge branch 'feature-iterative-eigendecomposition' into 'develop'

ENH: Robust Iterative Eigendecomposition and Parallel Low-Memory Dynamic Mode Decomposition

See merge request Development/openfoam!353

applications/test/TestTools/TestTools.H

@@ -0,0 +1,305 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Various tools for test applications.
+
+\*---------------------------------------------------------------------------*/
+
+using namespace Foam;
+
+#include "Random.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Create a non-complex random Matrix.
+template<class MatrixType>
+typename std::enable_if
+<
+    !std::is_same<complex, typename MatrixType::cmptType>:: value,
+    MatrixType
+>::type makeRandomMatrix
+(
+    const labelPair& dims,
+    Random& rndGen
+)
+{
+    MatrixType mat(dims);
+
+    std::generate
+    (
+        mat.begin(),
+        mat.end(),
+        [&]{return rndGen.GaussNormal<scalar>();}
+    );
+
+    return mat;
+}
+
+
+// Create a complex random Matrix.
+template<class MatrixType>
+typename std::enable_if
+<
+    std::is_same<complex, typename MatrixType::cmptType>:: value,
+    MatrixType
+>::type makeRandomMatrix
+(
+    const labelPair& dims,
+    Random& rndGen
+)
+{
+    MatrixType mat(dims);
+
+    for (auto& x : mat)
+    {
+        x = complex(rndGen.GaussNormal<scalar>(), rndGen.GaussNormal<scalar>());
+    }
+
+    return mat;
+}
+
+
+// Copy an initializer list into a DiagonalMatrix
+template<class Type>
+void assignMatrix
+(
+    UList<Type>& A,
+    std::initializer_list<Type> list
+)
+{
+    std::copy(list.begin(), list.end(), A.begin());
+}
+
+
+// Copy an initializer list into a SymmetricSquareMatrix.
+template<class Form, class Type>
+void assignMatrix
+(
+    Matrix<Form, Type>& A,
+    std::initializer_list<typename Matrix<Form, Type>::cmptType> list
+)
+{
+    const label nargs = list.size();
+
+    if (nargs != A.size())
+    {
+        FatalErrorInFunction
+            << "Mismatch in matrix dimension ("
+            << A.m() << ", "
+            << A.n() << ") and number of args (" << nargs << ')' << nl
+            << exit(FatalError);
+     }
+
+    std::copy(list.begin(), list.end(), A.begin());
+}
+
+
+// Return a copy of the Matrix collapsed into one dimension.
+template<class Form, class Type>
+List<Type> flt
+(
+    const Matrix<Form, Type>& A,
+    const bool rowMajorOrder = true
+)
+{
+   List<Type> flatMatrix(A.size());
+
+   if (rowMajorOrder)
+   {
+       std::copy(A.cbegin(), A.cend(), flatMatrix.begin());
+   }
+   else
+   {
+       for (label j = 0; j < A.n(); ++j)
+       {
+           for (label i = 0; i < A.m(); ++i)
+           {
+               flatMatrix[i + j*A.m()] = A(i, j);
+           }
+       }
+   }
+
+   return flatMatrix;
+}
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+template<class Type>
+typename std::enable_if
+<
+    std::is_same<floatScalar, Type>::value ||
+    std::is_same<doubleScalar, Type>::value ||
+    std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar absTol = 0,    //<! useful for cmps near zero
+    const scalar relTol = 1e-8  //<! are values the same within 8 decimals
+)
+{
+    Info<< msg << x << "?=" << y << endl;
+
+    unsigned nFail = 0;
+
+    if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+    {
+        ++nFail;
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Compare two containers elementwise, and print output.
+// Do ++nFail_ if two components are not equal within a given tolerance.
+// The function is converted from PEP-485
+template<class Type>
+typename std::enable_if
+<
+    !std::is_same<floatScalar, Type>::value &&
+    !std::is_same<doubleScalar, Type>::value &&
+    !std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar absTol = 0,
+    const scalar relTol = 1e-8
+)
+{
+    Info<< msg << x << "?=" << y << endl;
+
+    unsigned nFail = 0;
+
+    for (label i = 0; i < x.size(); ++i)
+    {
+        if (max(absTol, relTol*max(mag(x[i]), mag(y[i]))) < mag(x[i] - y[i]))
+        {
+            ++nFail;
+        }
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Compare two containers elementwise, and print output.
+// Do ++nFail_ if two components are not equal within a given tolerance.
+// The function is converted from PEP-485
+template<class Type1, class Type2>
+typename std::enable_if
+<
+    !std::is_same<floatScalar, Type1>::value &&
+    !std::is_same<doubleScalar, Type1>::value &&
+    !std::is_same<complex, Type1>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type1& x,
+    const Type2& y,
+    const scalar absTol = 0,
+    const scalar relTol = 1e-8
+)
+{
+    Info<< msg << x << "?=" << y << endl;
+
+    unsigned nFail = 0;
+
+    for (label i = 0; i < x.size(); ++i)
+    {
+        if (max(absTol, relTol*max(mag(x[i]), mag(y[i]))) < mag(x[i] - y[i]))
+        {
+            ++nFail;
+        }
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Compare two Booleans, and print output.
+// Do ++nFail_ if two Booleans are not equal.
+void cmp
+(
+    const word& msg,
+    const bool x,
+    const bool y
+)
+{
+    Info<< msg << x << "?=" << y << endl;
+
+    unsigned nFail = 0;
+
+    if (x != y)
+    {
+        ++nFail;
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// ************************************************************************* //

applications/test/matrices/DiagonalMatrix/Make/files

@@ -0,0 +1,3 @@
+Test-DiagonalMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-DiagonalMatrix

applications/test/matrices/DiagonalMatrix/Make/options

@@ -0,0 +1,2 @@
+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */

applications/test/matrices/DiagonalMatrix/Test-DiagonalMatrix.C

@@ -0,0 +1,231 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-DiagonalMatrix
+
+Description
+    Tests for \c DiagonalMatrix constructors, member functions and global
+    functions using \c floatScalar, \c doubleScalar, and \c complex base types.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' if no theoretical
+    cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+\*---------------------------------------------------------------------------*/
+
+#include "DiagonalMatrix.H"
+#include "RectangularMatrix.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+#include "TestTools.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Create each constructor of DiagonalMatrix<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct empty from size:" << nl;
+        const DiagonalMatrix<Type> A(5);
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct from size and initialise all elems to zero:" << nl;
+        const DiagonalMatrix<Type> A(5, Zero);
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct from size and initialise all elems to value" << nl;
+        const DiagonalMatrix<Type> A(5, Type(8));
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct from the diagonal of a Matrix" << nl;
+        const RectangularMatrix<Type> M(3, 5, Zero);
+        const DiagonalMatrix<Type> A(M);
+        Info<< A << endl;
+    }
+}
+
+
+// Execute each member function of DiagonalMatrix<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    DiagonalMatrix<Type> A(3, Zero);
+    assignMatrix(A, {Type(1), Type(2), Type(-3)});
+
+    Info<< "# Operand: " << nl
+        << "  DiagonalMatrix = " << A << endl;
+
+
+    {
+        Info<< "# Return the matrix inverse into itself:" << nl;
+        A.invert();
+        cmp
+        (
+            "  DiagonalMatrix<Type>.invert() = ",
+            A,
+            List<Type>({Type(1), Type(0.5), Type(-0.333333)}),
+            1e-6
+        );
+    }
+    {
+        Info<< "# Sort:" << nl;
+
+        DiagonalMatrix<Type> B(5, Zero);
+        assignMatrix(B, {Type(1), Type(2), Type(-3), Type(5), Type(1.01)});
+
+        auto descend = [&](Type a, Type b){ return mag(a) > mag(b); };
+        const List<label> sortPermutation(B.sortPermutation(descend));
+        cmp
+        (
+            "  Return a sort permutation labelList according to "
+            "a given comparison on the diagonal entries",
+            sortPermutation,
+            List<label>({3, 2, 1, 4, 0})
+        );
+
+        DiagonalMatrix<Type> sortedB0(5, Zero);
+        assignMatrix
+        (
+            sortedB0,
+            {
+                Type(5),
+                        Type(-3),
+                                 Type(2),
+                                         Type(1.01),
+                                                    Type(1)
+            }
+        );
+        const DiagonalMatrix<Type> sortedB1
+        (
+            applyPermutation(B, sortPermutation)
+        );
+        cmp
+        (
+            "  Return Matrix column-reordered according to "
+            "a given permutation labelList",
+            sortedB0,
+            sortedB1
+        );
+
+        DiagonalMatrix<Type> cpB(B);
+        cpB.applyPermutation(sortPermutation);
+        cmp
+        (
+            "  Column-reorder this Matrix according to "
+            "a given permutation labelList",
+            sortedB0,
+            cpB
+        );
+    }
+
+}
+
+
+// Execute each global function of DiagonalMatrix<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+    DiagonalMatrix<Type> A(3, Zero);
+    assignMatrix(A, {Type(1), Type(2), Type(-3)});
+
+    Info<< "# Operand: " << nl
+        << "  DiagonalMatrix = " << A << nl << endl;
+
+
+    cmp
+    (
+        "  Inverse = ",
+        inv(A),
+        List<Type>({Type(1), Type(0.5), Type(-0.333333)}),
+        1e-6
+    );
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global functions: "<< typeID[I] << " ##" << nl;
+    test_global_funcs(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main()
+{
+    const std::tuple<floatScalar, doubleScalar, complex> types
+    (
+        std::make_tuple(Zero, Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "DiagonalMatrix<floatScalar>",
+        "DiagonalMatrix<doubleScalar>",
+        "DiagonalMatrix<complex>"
+    });
+
+    run_tests(types, typeID);
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/matrices/EigenMatrix/Make/files

@@ -0,0 +1,3 @@
+Test-EigenMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-EigenMatrix

applications/test/matrices/EigenMatrix/Make/options

@@ -0,0 +1,2 @@
+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */

applications/test/matrices/EigenMatrix/Test-EigenMatrix.C

@@ -0,0 +1,556 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is derivative work of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-EigenMatrix
+
+Description
+    Tests for \c EigenMatrix constructors, and member functions
+    using \c floatScalar, and \c doubleScalar base types.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' if no theoretical
+    cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+\*---------------------------------------------------------------------------*/
+
+#include "scalarMatrices.H"
+#include "RectangularMatrix.H"
+#include "SquareMatrix.H"
+#include "complex.H"
+#include "IOmanip.H"
+#include "EigenMatrix.H"
+#include "TestTools.H"
+#include <algorithm>
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Return the absolute tolerance value for bitwise comparisons of floatScalars
+floatScalar getTol(floatScalar)
+{
+    return 1e-2;
+}
+
+
+// Return the absolute tolerance value for bitwise comparisons of doubleScalars
+doubleScalar getTol(doubleScalar)
+{
+    return 1e-10;
+}
+
+
+// Create each constructor of EigenMatrix<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct from a SquareMatrix<Type>" << nl;
+        const SquareMatrix<Type> A(5, Zero);
+        const EigenMatrix<Type> EM(A);
+    }
+    {
+        Info<< "# Construct from a SquareMatrix<Type> and symmetry flag" << nl;
+        const SquareMatrix<Type> A(5, Zero);
+        const EigenMatrix<Type> EM(A, true);
+    }
+}
+
+
+// Execute each member function of EigenMatrix<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    SquareMatrix<Type> A(3, Zero);
+    assignMatrix
+    (
+        A,
+        {
+            Type(1), Type(2), Type(3),
+            Type(4), Type(5), Type(6),
+            Type(7), Type(8), Type(9)
+        }
+    );
+
+    Info<< "# Operand: " << nl
+        << "  SquareMatrix = " << A << endl;
+
+
+    // Since eigenvalues are unique and eigenvectors are not unique,
+    // the bitwise comparisons are limited to eigenvalue computations.
+    // Here, only the execution of the functions is tested, rather than
+    // the verification of the eigendecomposition through theoretical relations.
+    {
+        const EigenMatrix<Type> EM(A);
+        const DiagonalMatrix<Type> EValsRe = EM.EValsRe();
+        const DiagonalMatrix<Type>& EValsIm = EM.EValsIm();
+        const SquareMatrix<Type>& EVecs = EM.EVecs();
+
+        cmp
+        (
+            "  Return real eigenvalues or real part of complex eigenvalues = ",
+            EValsRe,
+            List<Type>
+            ({
+                Type(16.116844),
+                Type(-1.116844),
+                Type(0)
+            }),
+            getTol(Type(0)),
+            1e-6
+        );
+        cmp
+        (
+            "  Return zero-matrix for real eigenvalues "
+            "or imaginary part of complex eigenvalues = ",
+            EValsIm,
+            List<Type>(3, Zero)
+        );
+        Info<< "  Return eigenvectors matrix = " << EVecs << endl;
+    }
+}
+
+
+// Test the relation: "sum(eigenvalues) = trace(A)"
+// w.wiki/4zs (Retrieved: 16-06-19) # Item-1
+template<class Type>
+void test_eigenvalues_sum
+(
+    const SquareMatrix<Type>& A,
+    const DiagonalMatrix<Type>& EValsRe
+)
+{
+    const Type trace = A.trace();
+    // Imaginary part of complex conjugates cancel each other
+    const Type EValsSum = sum(EValsRe);
+    Info<< "  # A.mRows = " << A.m() << nl;
+    cmp
+    (
+        "  # sum(eigenvalues) = trace(A) = ",
+        EValsSum,
+        trace,
+        getTol(Type(0))
+    );
+}
+
+
+// Test the relation: "prod(eigenvalues) = det(A)"
+// w.wiki/4zs (Retrieved: 16-06-19) # Item-2
+// Note that the determinant computation may fail
+// which is not a suggestion that eigendecomposition fails
+template<class Type>
+void test_eigenvalues_prod
+(
+    const SquareMatrix<Type>& A,
+    const DiagonalMatrix<Type>& EValsRe,
+    const DiagonalMatrix<Type>& EValsIm
+)
+{
+    const Type determinant = mag(det(A));
+    Type EValsProd = Type(1);
+
+    if (EValsIm.empty())
+    {
+        for (label i = 0; i < EValsRe.size(); ++i)
+        {
+            EValsProd *= Foam::sqrt(sqr(EValsRe[i]));
+        }
+    }
+    else
+    {
+        for (label i = 0; i < EValsRe.size(); ++i)
+        {
+            EValsProd *= Foam::sqrt(sqr(EValsRe[i]) + sqr(EValsIm[i]));
+        }
+    }
+    cmp
+    (
+        "  # prod(eigenvalues) = det(A) = ",
+        EValsProd,
+        determinant,
+        getTol(Type(0))
+    );
+}
+
+
+// Test eigenvalues in eigendecomposition relations
+// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307)
+template<class Type>
+void test_eigenvalues(Type)
+{
+    Random rndGen(1234);
+    const label numberOfTests = 20;
+
+    // Non-symmetric
+    for (label i = 0; i < numberOfTests; ++i)
+    {
+        const label mRows = rndGen.position(100, 200);
+        const labelPair m(mRows, mRows);
+
+        const SquareMatrix<Type> A
+        (
+            makeRandomMatrix<SquareMatrix<Type>>(m, rndGen)
+        );
+
+        const EigenMatrix<Type> EM(A);
+        const DiagonalMatrix<Type>& EValsRe = EM.EValsRe();
+
+        test_eigenvalues_sum(A, EValsRe);
+
+        // LUDecompose does not work with floatScalar at the time of writing,
+        // hence det function. Once LUDecompose is refactored, comment out below
+        // const DiagonalMatrix<Type>& EValsIm = EM.EValsIm();
+        // test_eigenvalues_prod(A, EValsRe, EValsIm);
+    }
+
+    // Symmetric
+    for (label i = 0; i < numberOfTests; ++i)
+    {
+        const label mRows = rndGen.position(100, 200);
+        const labelPair m(mRows, mRows);
+
+        SquareMatrix<Type> A
+        (
+            makeRandomMatrix<SquareMatrix<Type>>(m, rndGen)
+        );
+        // Symmetrise with noise
+        for (label n = 0; n < A.n() - 1; ++n)
+        {
+            for (label m = A.m() - 1; m > n; --m)
+            {
+                A(n, m) = A(m, n) + SMALL;
+            }
+        }
+
+        const bool symmetric = true;
+        const EigenMatrix<Type> EM(A, symmetric);
+        const DiagonalMatrix<Type>& EValsRe = EM.EValsRe();
+
+        test_eigenvalues_sum(A, EValsRe);
+    }
+}
+
+
+// Test the relation: "(A & EVec - EVal*EVec) = 0"
+template<class Type>
+void test_characteristic_eq
+(
+    const SquareMatrix<Type>& Aorig,
+    const DiagonalMatrix<Type>& EValsRe,
+    const DiagonalMatrix<Type>& EValsIm,
+    const SquareMatrix<complex>& EVecs
+)
+{
+    SquareMatrix<complex> A(Aorig.m());
+    auto convertToComplex = [&](const scalar& val) { return complex(val); };
+    std::transform
+    (
+        Aorig.cbegin(),
+        Aorig.cend(),
+        A.begin(),
+        convertToComplex
+    );
+
+    for (label i = 0; i < A.m(); ++i)
+    {
+        const RectangularMatrix<complex>& EVec(EVecs.subColumn(i));
+        const complex EVal(EValsRe[i], EValsIm[i]);
+        const RectangularMatrix<complex> leftSide(A*EVec);
+        const RectangularMatrix<complex> rightSide(EVal*EVec);
+
+        cmp
+        (
+            "  # (A & EVec - EVal*EVec) = 0:",
+            flt(leftSide),
+            flt(rightSide),
+            getTol(Type(0))
+        );
+    }
+}
+
+
+// Test eigenvectors in eigendecomposition relations
+template<class Type>
+void test_eigenvectors(Type)
+{
+    Random rndGen(1234);
+    const label numberOfTests = 20;
+
+    // Non-symmetric
+    for (label i = 0; i < numberOfTests; ++i)
+    {
+        const label mRows = rndGen.position(100, 200);
+        const labelPair m(mRows, mRows);
+
+        const SquareMatrix<Type> A
+        (
+            makeRandomMatrix<SquareMatrix<Type>>(m, rndGen)
+        );
+
+        const EigenMatrix<Type> EM(A);
+        const DiagonalMatrix<Type>& EValsRe = EM.EValsRe();
+        const DiagonalMatrix<Type>& EValsIm = EM.EValsIm();
+        const SquareMatrix<complex> EVecs(EM.complexEVecs());
+
+        test_characteristic_eq(A, EValsRe, EValsIm, EVecs);
+    }
+
+    // Symmetric
+    for (label i = 0; i < numberOfTests; ++i)
+    {
+        const label mRows = rndGen.position(100, 200);
+        const labelPair m(mRows, mRows);
+
+        SquareMatrix<Type> A
+        (
+            makeRandomMatrix<SquareMatrix<Type>>(m, rndGen)
+        );
+        // Symmetrise with noise
+        for (label n = 0; n < A.n() - 1; ++n)
+        {
+            for (label m = A.m() - 1; m > n; --m)
+            {
+                A(n, m) = A(m, n) + SMALL;
+            }
+        }
+
+        const bool symmetric = true;
+        const EigenMatrix<Type> EM(A, symmetric);
+        const DiagonalMatrix<Type>& EValsRe = EM.EValsRe();
+        const DiagonalMatrix<Type>& EValsIm = EM.EValsIm();
+        const SquareMatrix<complex> EVecs(EM.complexEVecs());
+
+        test_characteristic_eq(A, EValsRe, EValsIm, EVecs);
+    }
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test eigenvalues: "<< typeID[I] <<" ##" << nl;
+    test_eigenvalues(std::get<I>(types));
+
+    Info<< nl << "    ## Test eigenvectors: "<< typeID[I] <<" ##" << nl;
+    test_eigenvectors(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main()
+{
+    Info<< setprecision(15);
+
+    const std::tuple<floatScalar, doubleScalar> types
+    (
+        std::make_tuple(Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "SquareMatrix<floatScalar>",
+        "SquareMatrix<doubleScalar>"
+    });
+
+    run_tests(types, typeID);
+
+
+    Info<< nl << "    ## Test corner cases ##" << endl;
+    {
+        Info<< nl  << "  ## Rosser et al. (1951) matrix: ##" << nl;
+        // Rosser, J. B., Lanczos, C., Hestenes, M. R., & Karush, W. (1951).
+        // Separation of close eigenvalues of a real symmetric matrix.
+        // Jour. Research of the National Bureau of Standards, 47(4), 291-297.
+        // DOI:10.6028/jres.047.037
+        // 
+        // 8x8 symmetric square matrix consisting of close real eigenvalues
+        // ibid, p. 294
+        // {
+        //      1020.04901843, 1020.000, 1019.90195136, 1000.000,
+        //      1000.000, 0.09804864072, 0.000, -1020.0490
+        // }
+        // Note that prod(eigenvalues) != determinant(A) for this matrix
+        // via the LAPACK routine z/dgetrf
+
+        SquareMatrix<doubleScalar> A(8, Zero);
+        assignMatrix
+        (
+            A,
+            {
+                611, 196, -192, 407, -8, -52, -49, 29,
+                196, 899, 113, -192, -71, -43, -8, -44,
+                -192, 113, 899, 196, 61, 49, 8, 52,
+                407, -192, 196, 611, 8, 44, 59, -23,
+                -8, -71, 61, 8, 411, -599, 208, 208,
+                -52, -43, 49, 44, -599, 411, 208, 208,
+                -49, -8, 8, 59, 208, 208, 99, -911,
+                29, -44, 52, -23, 208, 208, -911, 99
+            }
+        );
+        const EigenMatrix<doubleScalar> EM(A);
+        const DiagonalMatrix<doubleScalar>& EValsRe = EM.EValsRe();
+        const DiagonalMatrix<doubleScalar>& EValsIm = EM.EValsIm();
+
+        test_eigenvalues_sum(A, EValsRe);
+
+        cmp
+        (
+            "  # Rosser et al. (1951) case, EValsRe = ",
+            EValsRe,
+            List<doubleScalar> // theoretical EValsRe
+            ({
+                -1020.0490, 0.000, 0.09804864072, 1000.000,
+                1000.000, 1019.90195136, 1020.000, 1020.04901843
+            }),
+            1e-3
+        );
+        cmp
+        (
+            "  # Rosser et al. (1951) case, EValsIm = ",
+            EValsIm,
+            List<doubleScalar>(8, Zero)
+        );
+    }
+    {
+        Info<< nl  << "  ## Test eigenvector unpacking: ##" << nl;
+
+        SquareMatrix<doubleScalar> A(3, Zero);
+        assignMatrix
+        (
+            A,
+            {
+                1, 2, 3,
+                -4, -5, 6,
+                7, -8, 9
+            }
+        );
+        const EigenMatrix<doubleScalar> EM(A);
+        const SquareMatrix<complex> complexEVecs(EM.complexEVecs());
+
+        SquareMatrix<complex> B(3, Zero);
+        assignMatrix
+        (
+            B,
+            {
+                complex(-0.373220280),
+                complex(0.417439996, 0.642691344),
+                complex(0.417439996, -0.642691344),
+                complex(-0.263919251),
+                complex(-1.165275867, 0.685068715),
+                complex(-1.165275867, -0.685068715),
+                complex(-0.889411744),
+                complex(-0.89990601, -0.3672785281),
+                complex(-0.89990601, 0.3672785281),
+            }
+        );
+
+        cmp
+        (
+            "  # ",
+            flt(complexEVecs),
+            flt(B)
+        );
+    }
+    {
+        Info<< nl << "  ## Test matrices with small values: ##" << nl;
+
+        const List<doubleScalar> epsilons
+        ({
+            0, SMALL, Foam::sqrt(SMALL), sqr(SMALL), Foam::cbrt(SMALL),
+            -SMALL, -Foam::sqrt(SMALL), -sqr(SMALL), -Foam::cbrt(SMALL)
+        });
+
+        Random rndGen(1234);
+        const label numberOfTests = 20;
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            const label mRows = rndGen.position(100, 200);
+ 
+            for (const auto& eps : epsilons)
+            {
+                const SquareMatrix<doubleScalar> A(mRows, eps);
+
+                const EigenMatrix<doubleScalar> EM(A);
+                const DiagonalMatrix<doubleScalar>& EValsRe = EM.EValsRe();
+                const DiagonalMatrix<doubleScalar>& EValsIm = EM.EValsIm();
+                const SquareMatrix<complex> EVecs(EM.complexEVecs());
+
+                test_eigenvalues_sum(A, EValsRe);
+                test_characteristic_eq(A, EValsRe, EValsIm, EVecs);
+            }
+        }
+    }
+    {
+        Info<< nl << "  ## Test matrices with repeating eigenvalues: ##" << nl;
+        SquareMatrix<doubleScalar> A(3, Zero);
+        assignMatrix
+        (
+            A,
+            {
+                0, 1, 1,
+                1, 0, 1,
+                1, 1, 0
+            }
+        );
+        const EigenMatrix<doubleScalar> EM(A);
+        const DiagonalMatrix<doubleScalar>& EValsRe = EM.EValsRe();
+        const DiagonalMatrix<doubleScalar>& EValsIm = EM.EValsIm();
+        const SquareMatrix<complex> EVecs(EM.complexEVecs());
+
+        test_eigenvalues_sum(A, EValsRe);
+        test_characteristic_eq(A, EValsRe, EValsIm, EVecs);
+    }
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+}

applications/test/matrices/Matrix/Test-Matrix.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

applications/test/matrices/QRMatrix/Test-QRMatrix.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,7 +37,6 @@ using namespace Foam::MatrixTools;
 #define RUNALL true
 const bool verbose = true;
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 void horizontalLine()

applications/test/matrices/RectangularMatrix/Make/files

@@ -0,0 +1,3 @@
+Test-RectangularMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-RectangularMatrix

applications/test/matrices/RectangularMatrix/Make/options

@@ -0,0 +1,2 @@
+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */

applications/test/matrices/RectangularMatrix/Test-RectangularMatrix.C

@@ -0,0 +1,878 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-RectangularMatrix
+
+Description
+    Tests for \c RectangularMatrix constructors, member functions, member
+    operators, global functions, global operators, and friend functions using
+    \c floatScalar, \c doubleScalar, and \c complex base types.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' if no theoretical
+    cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+    For \c complex base type, the cross-checks do only involve zero imag part.
+
+Note
+    Pending tests were tagged as "## Pending ##".
+
+\*---------------------------------------------------------------------------*/
+
+#include "RectangularMatrix.H"
+#include "SquareMatrix.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+#include "IOmanip.H"
+#include "TestTools.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Create each constructor of RectangularMatrix<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct a square matrix (rows == columns):" << nl;
+        const RectangularMatrix<Type> A(5);
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct given number of rows/columns:" << nl;
+        const RectangularMatrix<Type> A(2, 4);
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct given number of rows/columns "
+            << "initializing all elements to zero:" << nl;
+        const RectangularMatrix<Type> A(2, 4, Zero);
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct given number of rows/columns "
+            << "initializing all elements to a value:" << nl;
+        const RectangularMatrix<Type> A(2, 4, Type(3));
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct given number of rows/columns "
+            << "initializing all elements to zero, and diagonal to one:" << nl;
+        const RectangularMatrix<Type> A(labelPair(2, 4), I);
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct given number of rows/columns "
+            << "by using a label pair:" << nl;
+        const RectangularMatrix<Type> A(labelPair(2, 4));
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct given number of rows/columns by using a label pair "
+            << "and initializing all elements to zero:" << nl;
+        const RectangularMatrix<Type> A(labelPair(2, 4), Zero);
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct given number of rows/columns by using a label pair "
+            << "and initializing all elements to the given value:" << nl;
+        const RectangularMatrix<Type> A(labelPair(2, 4), Type(3));
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct from a block of another matrix:" << nl;
+        const RectangularMatrix<Type> B(labelPair(2, 4), Type(3));
+        const RectangularMatrix<Type> A(B.subMatrix(1, 1));
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct from a block of another matrix:" << nl;
+        RectangularMatrix<Type> B(labelPair(2, 4), Type(3));
+        const RectangularMatrix<Type> A(B.subMatrix(1, 1));
+        Info<< A << endl;
+    }
+    {
+        Info<< "#: Construct as copy of a square matrix:" << nl;
+        const SquareMatrix<Type> B(5);
+        const RectangularMatrix<Type> A(B);
+        Info<< A << endl;
+    }
+}
+
+
+// Execute each member function of RectangularMatrix<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    RectangularMatrix<Type> A(labelPair(2, 3), Zero);
+    assignMatrix
+    (
+        A,
+        {
+            Type(1), Type(-2.2), Type(-3.4),
+            Type(-0.35), Type(1), Type(5)
+        }
+    );
+
+    Info<< "# Operand: " << nl
+        << "  RectangularMatrix = " << A << endl;
+
+
+    // Matrix.H
+    {
+        Info<< "# Access:" << nl;
+
+        cmp("  The number of rows = ", Type(A.m()), Type(2));
+        cmp("  The number of columns = ", Type(A.n()), Type(3));
+        cmp("  The number of elements in Matrix = ", Type(A.size()), Type(6));
+        cmp("  Return row/column sizes = ", A.sizes(), labelPair(2, 3));
+        cmp("  Return true if Matrix is empty = ", A.empty(), false);
+        cmp
+        (
+            "  Return const pointer to the first data elem = ",
+            *(A.cdata()),
+            Type(1)
+        );
+        cmp
+        (
+            "  Return pointer to the first data elem = ",
+            *(A.data()),
+            Type(1)
+        );
+        cmp
+        (
+            "  Return const pointer to data in the specified row = ",
+            *(A.rowData(1)),
+            Type(-0.35)
+        );
+        cmp
+        (
+            "  Return pointer to data in the specified row = ",
+            *(A.rowData(1)),
+            Type(-0.35)
+        );
+        const Type& a = A.at(4);
+        cmp("  Linear addressing const element access = ", a, Type(1));
+        Type& b = A.at(4);
+        cmp("  Linear addressing element access = ", b, Type(1));
+    }
+    {
+        Info<< "# Block access (const):" << nl;
+
+        const RectangularMatrix<Type> Acol(A.subColumn(1));
+        cmp
+        (
+            "  Return const column or column's subset of Matrix = ",
+            flt(Acol),
+            List<Type>({Type(-2.2), Type(1)})
+        );
+
+        const RectangularMatrix<Type> Arow(A.subRow(1));
+        cmp
+        (
+            "  Return const row or row's subset of Matrix = ",
+            flt(Arow),
+            List<Type>({Type(-0.35), Type(1), Type(5)})
+        );
+
+        const RectangularMatrix<Type> Amat(A.subMatrix(1, 1));
+        cmp
+        (
+            "  Return const sub-block of Matrix = ",
+            flt(Amat),
+            List<Type>({Type(1), Type(5)})
+        );
+    }
+    {
+        Info<< "# Block access (non-const):" << nl;
+
+        RectangularMatrix<Type> Acol(A.subColumn(1));
+        cmp
+        (
+            "  Return column or column's subset of Matrix = ",
+            flt(Acol),
+            List<Type>({Type(-2.2), Type(1)})
+        );
+
+        RectangularMatrix<Type> Arow(A.subRow(1));
+        cmp
+        (
+            "  Return row or row's subset of Matrix = ",
+            flt(Arow),
+            List<Type>({Type(-0.35), Type(1), Type(5)})
+        );
+
+        RectangularMatrix<Type> Amat(A.subMatrix(1, 1));
+        cmp
+        (
+            "  Return sub-block of Matrix = ",
+            flt(Amat),
+            List<Type>({Type(1), Type(5)})
+        );
+    }
+    {
+        Info<< "# Check:" << nl;
+
+        A.checki(0);
+        A.checkj(1);
+        A.checkSize();
+        cmp("  Check all entries have identical values = ", A.uniform(), false);
+    }
+    {
+        Info<< "# Edit:" << nl;
+
+        RectangularMatrix<Type> cpA(A);
+        cpA.clear();
+        cmp("  Clear Matrix, i.e. set sizes to zero = ", cpA.empty(), true);
+
+        RectangularMatrix<Type> cpA1(A);
+        cmp
+        (
+            "  Release storage management of Matrix contents by transferring "
+            "management to a List = ",
+            (cpA1.release()),
+            List<Type>
+            ({
+                Type(1), Type(-2.2), Type(-3.4), Type(-0.35), Type(1), Type(5)
+            })
+        );
+
+        RectangularMatrix<Type> cpA2(A);
+        cpA2.swap(cpA1);
+        cmp("  Swap contents = ", flt(cpA1), flt(A));
+
+        cpA2.transfer(cpA1);
+        cmp
+        (
+            "  Transfer the contents of the argument Matrix into this Matrix "
+            "and annul the argument MatrixSwap contents = ",
+            flt(cpA2),
+            flt(A)
+        );
+
+        cpA2.resize(1, 2);
+        cmp
+        (
+            "  Change Matrix dimensions, preserving the elements = ",
+            flt(cpA2),
+            List<Type>({Type(1), Type(-2.2)})
+        );
+
+        RectangularMatrix<Type> cpA3(A);
+        cpA3.setSize(1, 2);
+        cmp
+        (
+            "  Change Matrix dimensions, preserving the elements = ",
+            flt(cpA3),
+            List<Type>({Type(1), Type(-2.2)})
+        );
+
+        RectangularMatrix<Type> cpA4(A);
+        cpA4.shallowResize(1, 2);
+        cmp
+        (
+            "  Resize Matrix without reallocating storage (unsafe) = ",
+            flt(cpA4),
+            List<Type>({Type(1), Type(-2.2)})
+        );
+
+        RectangularMatrix<Type> smallA(labelPair(2, 3), Type(VSMALL));
+        smallA.round();
+        cmp
+        (
+            "  Round elements with mag smaller than tol (SMALL) to zero = ",
+            flt(smallA),
+            List<Type>(6, Zero)
+        );
+    }
+    {
+        Info<< "# Operations:" << nl;
+
+        cmp("  Transpose = ", flt((A.T()).T()), flt(A));
+
+        RectangularMatrix<complex> cA(labelPair(2, 2), Zero);
+        assignMatrix
+        (
+            cA,
+            {
+                complex(2, -1.2), complex(4.1, 0.4),
+                complex(8.1, -1.25), complex(7.3, -1.4)
+            }
+        );
+        RectangularMatrix<complex> cAT(labelPair(2, 2), Zero);
+        assignMatrix
+        (
+            cAT,
+            {
+                complex(2, 1.2), complex(8.1, 1.25),
+                complex(4.1, -0.4), complex(7.3, 1.4)
+            }
+        );
+        cmp("  Hermitian transpose = ", flt(cA.T()), flt(cAT));
+
+        RectangularMatrix<Type> B(3, 3, Zero);
+        assignMatrix
+        (
+            B,
+            {
+                Type(4.1), Type(12.5), Type(-16.3),
+                Type(-192.3), Type(-9.1), Type(-3.0),
+                Type(1.0), Type(5.02), Type(-4.4)
+            }
+        );
+        const RectangularMatrix<Type> cpB(B);
+        const List<Type> lst({Type(1), Type(2), Type(3)});
+        const Field<Type> out1(B*lst);
+        const Field<Type> out2(B.Amul(lst));
+        const Field<Type> out3(lst*B);
+        const Field<Type> out4(B.Tmul(lst));
+        const Field<Type> rsult1({Type(-19.8), Type(-219.5), Type(-2.16)});
+        const Field<Type> rsult2({Type(-377.5), Type(9.36), Type(-35.5)});
+        cmp
+        (
+            "  Right-multiply Matrix by a List (A * x) = ",
+            out1,
+            rsult1,
+            1e-4
+        );
+        cmp
+        (
+            "  Right-multiply Matrix by a column vector A.Amul(x) = ",
+            out1,
+            out2,
+            1e-4
+        );
+        cmp
+        (
+            "  Left-multiply Matrix by a List (x * A) = ",
+            out3,
+            rsult2,
+            1e-4
+        );
+        cmp
+        (
+            "  Left-multiply Matrix by a row vector A.Tmul(x) = ",
+            out4,
+            rsult2,
+            1e-4
+        );
+
+        List<Type> diagB1({Type(4.1), Type(-9.1), Type(-4.4)});
+        cmp
+        (
+            "  Extract the diagonal elements = ",
+            B.diag(),
+            diagB1
+        );
+
+        List<Type> diagB2({Type(-100), Type(-100), Type(-100)});
+        B.diag(diagB2);
+        cmp
+        (
+            "  Assign diagonal of Matrix = ",
+            B.diag(),
+            diagB2
+        );
+
+        B = cpB;
+        cmp("  Trace = ", B.trace(), Type(-9.4), 1e-4);
+        cmp
+        (
+            "  Return L2-Norm of chosen column = ",
+            B.columnNorm(0),
+            192.34630227794867,
+            1e-4
+        );
+        cmp
+        (
+            "  Return L2-Norm of chosen column (noSqrt=true) = ",
+            B.columnNorm(0, true),
+            36997.1,
+            1e-2
+        );
+        cmp
+        (
+            "  Return Frobenius norm of Matrix = ",
+            B.norm(),
+            193.7921835369012,
+            1e-4
+        );
+    }
+}
+
+
+// Execute each member operators of RectangularMatrix<Type>, and print output
+template<class Type>
+void test_member_opers(Type)
+{
+    RectangularMatrix<Type> A(labelPair(2, 4), I);
+    const RectangularMatrix<Type> cpA(A);
+
+    Info<< "# Operand: " << nl
+        << "  RectangularMatrix = " << A << endl;
+
+
+    A = Zero;
+    cmp("  Assign all elements to zero =", flt(A), List<Type>(8, Zero));
+    A = cpA;
+
+    A = Type(5);
+    cmp("  Assign all elements to value =", flt(A), List<Type>(8, Type(5)));
+    A = cpA;
+
+    const Type* a = A[1];
+    cmp
+    (
+        "  Return const pointer to data in the specified row = ",
+        *a,
+        Type(0)
+    );
+
+    Type* b = A[1];
+    cmp
+    (
+        "  Return pointer to data in the specified row = ",
+        *b,
+        Type(0)
+    );
+
+    const Type& c = A(1, 1);
+    cmp
+    (
+        "  (i, j) const element access operator = ",
+        c,
+        Type(1)
+    );
+
+    Type& d = A(1, 1);
+    cmp
+    (
+        "  (i, j) element access operator = ",
+        d,
+        Type(1)
+    );
+
+    //  ## Pending ##
+    //  Copy assignment
+    //  Move assignment
+    //  Assignment to a block of another Matrix
+    //  Assignment to a block of another Matrix
+    //  #############
+
+    A = Zero;
+    cmp
+    (
+        "  Assignment of all elements to zero = ",
+        flt(A),
+        flt(RectangularMatrix<Type>(2, 4, Zero))
+    );
+    A = cpA;
+
+    A = Type(-1.2);
+    cmp
+    (
+        "  Assignment of all elements to the given value = ",
+        flt(A),
+        flt(RectangularMatrix<Type>(2, 4, Type(-1.2)))
+    );
+    A = cpA;
+
+    A += cpA;
+    cmp
+    (
+        "  Matrix addition =",
+        flt(A),
+        flt(Type(2)*cpA)
+    );
+    A = cpA;
+
+    A -= cpA;
+    cmp
+    (
+        "  Matrix subtraction = ",
+        flt(A),
+        flt(RectangularMatrix<Type>(2, 4, Zero))
+    );
+    A = cpA;
+
+    A = Zero;
+    A += Type(5);
+    cmp
+    (
+        "  Matrix scalar addition = ",
+        flt(A),
+        flt(RectangularMatrix<Type>(2, 4, Type(5)))
+    );
+    A = cpA;
+
+    A = Zero;
+    A -= Type(5);
+    cmp
+    (
+        "  Matrix scalar subtraction = ",
+        flt(A),
+        flt(RectangularMatrix<Type>(2, 4, Type(-5)))
+    );
+    A = cpA;
+
+    A = Type(1);
+    A *= Type(5);
+    cmp
+    (
+        "  Matrix scalar multiplication = ",
+        flt(A),
+        flt(RectangularMatrix<Type>(2, 4, Type(5)))
+    );
+    A = cpA;
+
+    A = Type(1);
+    A /= Type(5);
+    cmp
+    (
+        "  Matrix scalar division = ",
+        flt(A),
+        flt(RectangularMatrix<Type>(2, 4, Type(0.2)))
+    );
+    A = cpA;
+
+    A(0,0) = Type(-10.5);
+    Type* i0 = A.begin();
+    cmp
+    (
+        "  Return an iterator to begin traversing a Matrix = ",
+        *i0,
+        Type(-10.5)
+    );
+
+    A(1,3) = Type(20);
+    Type* i1 = A.end();
+    cmp
+    (
+        "  Return an iterator to end traversing a Matrix = ",
+        *(--i1),
+        Type(20)
+    );
+
+    const Type* i2 = A.cbegin();
+    cmp
+    (
+        "  Return const iterator to begin traversing a Matrix = ",
+        *i2,
+        Type(-10.5)
+    );
+
+    const Type* i3 = A.cend();
+    cmp
+    (
+        "  Return const iterator to end traversing a Matrix = ",
+        *(--i3),
+        Type(20)
+    );
+
+    const Type* i4 = A.begin();
+    cmp
+    (
+        "  Return const iterator to begin traversing a Matrix = ",
+        *i4,
+        Type(-10.5)
+    );
+
+    const Type* i5 = A.end();
+    cmp
+    (
+        "  Return const iterator to end traversing a Matrix = ",
+        *(--i5),
+        Type(20)
+    );
+
+    //  ## Pending ##
+    //  Read Matrix from Istream, discarding existing contents
+    //  Write Matrix, with line-breaks in ASCII when length exceeds shortLen
+    //  #############
+}
+
+
+// Execute each friend function of RectangularMatrix<Type>, and print output
+template<class Type>
+void test_friend_funcs(Type)
+{
+    const Field<Type> F
+    ({
+        Type(1), Type(-2.2),
+        Type(10), Type(-0.1)
+    });
+
+    Info<< "# Operand: " << nl
+        << "  Field = " << F << endl;
+
+    {
+        RectangularMatrix<Type> A(labelPair(4, 4), Zero);
+        assignMatrix
+        (
+            A,
+            {
+                Type(1), Type(-2.2), Type(10), Type(-0.1),
+                Type(-2.2), Type(4.84), Type(-22), Type(0.22),
+                Type(10), Type(-22), Type(100), Type(-1),
+                Type(-0.1), Type(0.22), Type(-1), Type(0.01)
+            }
+        );
+        cmp
+        (
+            "  Return the outer product of Field-Field as RectangularMatrix = ",
+            flt(outer(F, F)),
+            flt(A),
+            1e-6
+        );
+    }
+}
+
+
+// Execute each global function of RectangularMatrix<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+    RectangularMatrix<Type> A(labelPair(2, 3), Zero);
+    assignMatrix
+    (
+        A,
+        {
+            Type(1), Type(-2.2), Type(-3.4),
+            Type(-0.35), Type(10.32), Type(5)
+        }
+    );
+    const RectangularMatrix<Type> cpA(A);
+
+    Info<< "# Operand: " << nl
+        << "  RectangularMatrix = " << A << endl;
+
+
+    //  ## Pending ##
+    //  cmp("  Find max value in Matrix = ", max(A), Type(10.32));
+    //  cmp("  Find min value in Matrix = ", min(A), Type(-3.4));
+    //  cmp
+    //  (
+    //      "  Find min-max value in Matrix = ",
+    //      minMax(A),
+    //      MinMax<Type>(10.32, -3.4)
+    //  );
+    //  #############
+}
+
+
+// Execute each global operators of RectangularMatrix<Type>, and print output
+template<class Type>
+void test_global_opers(Type)
+{
+    RectangularMatrix<Type> A(labelPair(2, 3), Zero);
+    assignMatrix
+    (
+        A,
+        {
+            Type(1), Type(-2.2), Type(-3.4),
+            Type(-0.35), Type(10.32), Type(5)
+        }
+    );
+    const RectangularMatrix<Type> cpA(A);
+
+    Info<< "# Operand: " << nl
+        << "  RectangularMatrix = " << A << endl;
+
+
+    // Matrix.C
+    cmp("  Matrix negation = ", flt(-A), flt(Type(-1)*A));
+    cmp("  Matrix addition = ", flt(A + A), flt(Type(2)*A));
+    cmp("  Matrix subtraction = ", flt(A - A), flt(Type(0)*A));
+    cmp("  Scalar multiplication = ", flt(Type(2)*A), flt(A + A));
+    cmp("  Scalar multiplication = ", flt(A*Type(2)), flt(A + A));
+    cmp
+    (
+        "  Scalar addition = ",
+        flt(Type(0)*A + Type(1)),
+        flt(RectangularMatrix<Type>(A.sizes(), Type(1)))
+    );
+    cmp
+    (
+        "  Scalar addition = ",
+        flt(Type(1) + Type(0)*A),
+        flt(RectangularMatrix<Type>(A.sizes(), Type(1)))
+    );
+    cmp
+    (
+        "  Scalar subtraction = ",
+        flt(Type(0)*A - Type(1)),
+        flt(RectangularMatrix<Type>(A.sizes(), Type(-1)))
+    );
+    cmp
+    (
+        "  Scalar subtraction = ",
+        flt(Type(1) - Type(0)*A),
+        flt(RectangularMatrix<Type>(A.sizes(), Type(1)))
+    );
+    cmp
+    (
+        "  Scalar division = ",
+        flt((Type(1) + Type(0)*A)/Type(2.5)),
+        flt(RectangularMatrix<Type>(A.sizes(), Type(0.4)))
+    );
+
+    RectangularMatrix<Type> innerA(2, 2, Zero);
+    assignMatrix
+    (
+        innerA,
+        {
+            Type(17.4), Type(-40.054),
+            Type(-40.054), Type(131.6249)
+        }
+    );
+    const RectangularMatrix<Type> A1(A);
+    const RectangularMatrix<Type> A2(A.T());
+    cmp
+    (
+        "  Matrix-Matrix multiplication using ikj-algorithm = ",
+        flt(A1*A2),
+        flt(innerA),
+        1e-1
+    );
+    cmp
+    (
+        "  Implicit A.T()*B = ",
+        flt(A2&A2),
+        flt(innerA),
+        1e-1
+    );
+
+    RectangularMatrix<Type> outerA(3, 3, Zero);
+    assignMatrix
+    (
+        outerA,
+        {
+            Type(1.1225), Type(-5.812), Type(-5.15),
+            Type(-5.812), Type(111.3424), Type(59.08),
+            Type(-5.15), Type(59.08), Type(36.56)
+        }
+    );
+    cmp
+    (
+        "  Implicit A*B.T() = ",
+        flt(A2^A2),
+        flt(outerA),
+        1e-1
+    );
+
+    // MatrixI.H
+    const List<Type> lst({Type(1), Type(2), Type(-3)});
+    const List<Type> out1(A*lst);
+    const List<Type> out2(lst*A.T());
+    const List<Type> rslt1({Type(6.8), Type(5.29)});
+    const List<Type> rslt2({Type(6.8), Type(5.29)});
+    cmp
+    (
+        "  Matrix-vector multiplication (A * x), where x is a column vector = ",
+        out1,
+        rslt1,
+        1e-1
+    );
+    cmp
+    (
+        "  Vector-matrix multiplication (x * A), where x is a row vector = ",
+        out2,
+        rslt2,
+        1e-1
+    );
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test member opers: "<< typeID[I] <<" ##" << nl;
+    test_member_opers(std::get<I>(types));
+
+    Info<< nl << "    ## Test global functions: "<< typeID[I] << " ##" << nl;
+    test_global_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global operators: "<< typeID[I] << " ##" << nl;
+    test_global_opers(std::get<I>(types));
+
+    Info<< nl << "    ## Test friend funcs: "<< typeID[I] <<" ##" << nl;
+    test_friend_funcs(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main()
+{
+    Info<< setprecision(15);
+
+    const std::tuple<floatScalar, doubleScalar, complex> types
+    (
+        std::make_tuple(Zero, Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "RectangularMatrix<floatScalar>",
+        "RectangularMatrix<doubleScalar>",
+        "RectangularMatrix<complex>"
+    });
+
+    run_tests(types, typeID);
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/matrices/SquareMatrix/Make/files

@@ -0,0 +1,3 @@
+Test-SquareMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-SquareMatrix

applications/test/matrices/SquareMatrix/Make/options

@@ -0,0 +1,2 @@
+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */

applications/test/matrices/SymmetricSquareMatrix/Make/files

@@ -0,0 +1,3 @@
+Test-SymmetricSquareMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-SymmetricSquareMatrix

applications/test/matrices/SymmetricSquareMatrix/Make/options

@@ -0,0 +1,2 @@
+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */

applications/test/matrices/SymmetricSquareMatrix/Test-SymmetricSquareMatrix.C

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-SymmetricSquareMatrix
+
+Description
+    Tests for \c SymmetricSquareMatrix constructors, member functions, member
+    operators, global functions, global operators, and friend functions using
+    \c floatScalar, \c doubleScalar, and \c complex base types.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' if no theoretical
+    cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+    For \c complex base type, the cross-checks do only involve zero imag part.
+
+Note
+    Pending tests were tagged as "## Pending ##".
+
+\*---------------------------------------------------------------------------*/
+
+#include "scalarMatrices.H"
+#include "RectangularMatrix.H"
+#include "SquareMatrix.H"
+#include "SymmetricSquareMatrix.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+#include "IOmanip.H"
+#include "Random.H"
+#include "TestTools.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Create each constructor of SymmetricSquareMatrix<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct for given size (rows == cols):" << nl;
+        const SymmetricSquareMatrix<Type> A(2);
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct for given size (rows == cols) "
+            << "initializing all elements to zero:" << nl;
+        const SymmetricSquareMatrix<Type> A(2, Zero);
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct for given size (rows == cols) "
+            << "initializing all elements to a value:" << nl;
+        const SymmetricSquareMatrix<Type> A(2, Type(3));
+        Info<< A << endl;
+    }
+    {
+        Info<< "# Construct for given size (rows == cols) "
+            << "initializing to the identity matrix:" << nl;
+        const SymmetricSquareMatrix<Type> A(2, I);
+        Info<< A << endl;
+    }
+
+    //  ## Pending ##
+    //  Construct from Istream
+    //  Clone
+    //  #############
+}
+
+
+// Execute each member function of SymmetricSquareMatrix<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    //  ## Pending ##
+    //  #############
+}
+
+
+// Execute each member operators of SymmetricSquareMatrix<Type>, and print output
+template<class Type>
+void test_member_opers(Type)
+{
+    //  ## Pending ##
+    //  #############
+}
+
+
+// Execute each friend function of SymmetricSquareMatrix<Type>, and print output
+template<class Type>
+void test_friend_funcs(Type)
+{
+    //  ## Pending ##
+    //  #############
+}
+
+
+// Execute each global function of SymmetricSquareMatrix<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+    //  ## Pending ##
+    //  #############
+}
+
+
+// Execute each global operators of SymmetricSquareMatrix<Type>, and print output
+template<class Type>
+void test_global_opers(Type)
+{
+    //  ## Pending ##
+    //  #############
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test member opers: "<< typeID[I] <<" ##" << nl;
+    test_member_opers(std::get<I>(types));
+
+    Info<< nl << "    ## Test global functions: "<< typeID[I] << " ##" << nl;
+    test_global_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global operators: "<< typeID[I] << " ##" << nl;
+    test_global_opers(std::get<I>(types));
+
+    Info<< nl << "    ## Test friend funcs: "<< typeID[I] <<" ##" << nl;
+    test_friend_funcs(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main()
+{
+    Info<< setprecision(15);
+
+    const std::tuple<floatScalar, doubleScalar, complex> types
+    (
+        std::make_tuple(Zero, Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "SymmetricSquareMatrix<floatScalar>",
+        "SymmetricSquareMatrix<doubleScalar>",
+        "SymmetricSquareMatrix<complex>"
+    });
+
+    run_tests(types, typeID);
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,14 +28,7 @@ License
 
 #include "DiagonalMatrix.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type>
-inline Foam::DiagonalMatrix<Type>::DiagonalMatrix()
-:
-    List<Type>()
-{}
-
 
 template<class Type>
 Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label n)
@@ -74,8 +67,10 @@ Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& mat)
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 template<class Type>
-Foam::DiagonalMatrix<Type>& Foam::DiagonalMatrix<Type>::invert()
+void Foam::DiagonalMatrix<Type>::invert()
 {
     for (Type& val : *this)
     {
@@ -88,13 +83,75 @@ Foam::DiagonalMatrix<Type>& Foam::DiagonalMatrix<Type>::invert()
             val = Type(1)/val;
         }
     }
-
-    return this;
 }
 
 
 template<class Type>
-Foam::DiagonalMatrix<Type> Foam::inv(const DiagonalMatrix<Type>& mat)
+template<class CompOp>
+Foam::List<Foam::label> Foam::DiagonalMatrix<Type>::sortPermutation
+(
+    CompOp& compare
+) const
+{
+    List<label> p(this->size());
+    std::iota(p.begin(), p.end(), 0);
+    std::sort
+    (
+        p.begin(),
+        p.end(),
+        [&](label i, label j){ return compare((*this)[i], (*this)[j]); }
+    );
+
+    return p;
+}
+
+
+template<class Type>
+void Foam::DiagonalMatrix<Type>::applyPermutation(const List<label>& p)
+{
+    #ifdef FULLDEBUG
+    if (this->size() != p.size())
+    {
+        FatalErrorInFunction
+            << "Attempt to column-reorder according to an uneven list: " << nl
+            << "DiagonalMatrix diagonal size = " << this->size() << nl
+            << "Permutation list size = " << p.size() << nl
+            << abort(FatalError);
+    }
+    #endif
+
+    List<bool> pass(p.size(), false);
+
+    for (label i = 0; i < p.size(); ++i)
+    {
+        if (pass[i])
+        {
+            continue;
+        }
+        pass[i] = true;
+        label prev = i;
+        label j = p[i];
+        while (i != j)
+        {
+            Swap((*this)[prev], (*this)[j]);
+            pass[j] = true;
+            prev = j;
+            j = p[j];
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+//- Return the matrix inverse as a DiagonalMatrix if no elem is equal to zero
+template<class Type>
+DiagonalMatrix<Type> inv(const DiagonalMatrix<Type>& mat)
 {
     DiagonalMatrix<Type> Ainv(mat.size());
 
@@ -118,4 +175,41 @@ Foam::DiagonalMatrix<Type> Foam::inv(const DiagonalMatrix<Type>& mat)
 }
 
 
+//- Return Matrix column-reordered according to
+//- a given permutation labelList
+template<class Type>
+DiagonalMatrix<Type> applyPermutation
+(
+    const DiagonalMatrix<Type>& mat,
+    const List<label>& p
+)
+{
+    #ifdef FULLDEBUG
+    if (mat.size() != p.size())
+    {
+        FatalErrorInFunction
+            << "Attempt to column-reorder according to an uneven list: " << nl
+            << "DiagonalMatrix diagonal size = " << mat.size() << nl
+            << "Permutation list size = " << p.size() << nl
+            << abort(FatalError);
+    }
+    #endif
+
+    DiagonalMatrix<Type> reordered(mat.size());
+
+    label j = 0;
+    for (const label i : p)
+    {
+        reordered[j] = mat[i];
+        ++j;
+    }
+
+    return reordered;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 // ************************************************************************* //

src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,11 @@ Class
     Foam::DiagonalMatrix
 
 Description
-    A 2D diagonal matrix of objects of type \<Type\>, size (N x N)
+    A templated (N x N) diagonal matrix of objects of \<Type\>, effectively
+    containing N elements, derived from List.
+
+See also
+    Test-DiagonalMatrix.C
 
 SourceFiles
     DiagonalMatrix.C
@@ -39,6 +43,7 @@ SourceFiles
 #define DiagonalMatrix_H
 
 #include "List.H"
+#include <numeric>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,19 +64,27 @@ class DiagonalMatrix
 {
 public:
 
-    // Constructors
+    // Generated Methods
 
-        //- Construct null.
+        //- Default construct
-        DiagonalMatrix<Type>();
+        DiagonalMatrix() = default;
+
+        //- Copy construct
+        DiagonalMatrix(const DiagonalMatrix&) = default;
+
+        //- Copy assignment
+        DiagonalMatrix& operator=(const DiagonalMatrix&) = default;
+
+
+    // Constructors
 
         //- Construct empty from size
         explicit DiagonalMatrix<Type>(const label n);
 
-        //- Construct from size
+        //- Construct from size and initialise all elems to zero
-        //- initializing all elements to the given value
         DiagonalMatrix<Type>(const label n, const zero);
 
-        //- Construct from size and a value
+        //- Construct from size and initialise all elems to value
         DiagonalMatrix<Type>(const label n, const Type& val);
 
         //- Construct from the diagonal of a Matrix
@@ -81,19 +94,20 @@ public:
 
     // Member Functions
 
-        //- Invert the diagonal matrix and return itself
+        //- Return the matrix inverse into itself if no elem is equal to zero
-        DiagonalMatrix<Type>& invert();
+        void invert();
 
+        //- Return a sort permutation labelList according to
+        //- a given comparison on the diagonal entries
+        template<class CompOp>
+        List<label> sortPermutation(CompOp& compare) const;
+
+        //- Column-reorder this Matrix according to
+        //- a given permutation labelList
+        void applyPermutation(const List<label>& p);
 };
 
 
-// Global functions
-
-//- Return the diagonal Matrix inverse
-template<class Type>
-DiagonalMatrix<Type> inv(const DiagonalMatrix<Type>& mat);
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/matrices/EigenMatrix/EigenMatrix.H

@@ -0,0 +1,267 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Code created 2014-2018 by Alberto Passalacqua
+    Contributed 2018-07-31 to the OpenFOAM Foundation
+    Copyright (C) 2018 OpenFOAM Foundation
+    Copyright (C) 2019-2020 Alberto Passalacqua
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::EigenMatrix
+
+Description
+    EigenMatrix (i.e. eigendecomposition or spectral decomposition) decomposes
+    a diagonalisable nonsymmetric real square matrix into its canonical form,
+    whereby the matrix is represented in terms of its eigenvalues and
+    eigenvectors.
+    
+    The eigenvalue equation (i.e. eigenvalue problem) is:
+
+    \f[
+        A v = \lambda v
+    \f]
+
+    where
+    \vartable
+      A       | a diagonalisable square matrix of dimension m-by-m
+      v       | a (non-zero) vector of dimension m (right eigenvector)
+      \lambda | a scalar corresponding to v (eigenvalue)
+    \endvartable
+
+
+    If \c A is symmetric, the following relation is satisfied:
+
+    \f[
+        A = v*D*v^T
+    \f]
+
+    where
+    \vartable
+      D | diagonal real eigenvalue matrix
+      v | orthogonal eigenvector matrix
+    \endvartable
+
+    If \c A is not symmetric, \c D becomes a block diagonal matrix wherein
+    the real eigenvalues are present on the diagonal within 1-by-1 blocks, and
+    complex eigenvalues within 2-by-2 blocks, i.e. \f$\lambda + i \mu\f$ with
+    \f$[\lambda, \mu; -\mu, \lambda]\f$.
+
+    The columns of \c v represent eigenvectors corresponding to eigenvalues,
+    satisfying the eigenvalue equation. Even though eigenvalues of a matrix
+    are unique, eigenvectors of the matrix are not. For the same eigenvalue,
+    the corresponding eigenvector can be real or complex with non-unique
+    entries. In addition, the validity of the equation \f$A = v*D*v^T\f$
+    depends on the condition number of \c v, which can be ill-conditioned,
+    or singular for invalidated equations.
+
+    References:
+    \verbatim
+        OpenFOAM-compatible implementation:
+            Passalacqua, A., Heylmun, J., Icardi, M.,
+            Madadi, E., Bachant, P., & Hu, X. (2019).
+            OpenQBMM 5.0.1 for OpenFOAM 7, Zenodo.
+            DOI:10.5281/zenodo.3471804
+
+        Implementations for the functions:
+        'tridiagonaliseSymmMatrix', 'symmTridiagonalQL',
+        'Hessenberg' and 'realSchur' (based on ALGOL-procedure:tred2):
+            Wilkinson, J. H., & Reinsch, C. (1971).
+            In Bauer, F. L. & Householder A. S. (Eds.),
+            Handbook for Automatic Computation: Volume II: Linear Algebra.
+            (Vol. 186), Springer-Verlag Berlin Heidelberg.
+            DOI: 10.1007/978-3-642-86940-2
+
+        Explanations on how real eigenvectors
+        can be unpacked into complex domain:
+            Moler, C. (1998).
+            Re: Eigenvectors.
+            Retrieved from https://bit.ly/3ao4Wat
+
+        TNT/JAMA implementation:
+            Pozo, R. (1997).
+            Template Numerical Toolkit for linear algebra:
+            High performance programming with C++
+            and the Standard Template Library.
+            The International Journal of Supercomputer Applications
+            and High Performance Computing, 11(3), 251-263.
+            DOI:10.1177/109434209701100307
+
+            (No particular order) Hicklin, J., Moler, C., Webb, P.,
+            Boisvert, R. F., Miller, B., Pozo, R., & Remington, K. (2012).
+            JAMA: A Java Matrix Package 1.0.3.
+            Retrived from https://math.nist.gov/javanumerics/jama/
+    \endverbatim
+
+Note
+    - This implementation is an integration of the \c OpenQBMM \c eigenSolver
+    class (2019) without any changes to its internal mechanisms. Therefore, no
+    differences between EigenMatrix and \c eigenSolver (2019) classes should be
+    expected in terms of input-process-output operations.
+    - The \c OpenQBMM \c eigenSolver class derives almost completely from the
+    \c TNT/JAMA implementation, a public-domain library developed by \c NIST
+    and \c MathWorks from 1998 to 2012, available at
+    http://math.nist.gov/tnt/index.html (Retrieved June 6, 2020). Their
+    implementation was based upon \c EISPACK.
+    - The \c tridiagonaliseSymmMatrix, \c symmTridiagonalQL, \c Hessenberg and
+    \c realSchur methods are based on the \c Algol procedures \c tred2 by
+    Bowdler, Martin, Reinsch, and Wilkinson, Handbook for Auto. Comp.,
+    Vol. II-Linear Algebra, and the corresponding \c FORTRAN subroutine
+    in \c EISPACK.
+
+See also
+    Test-EigenMatrix.C
+
+SourceFiles
+    EigenMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef EigenMatrix_H
+#define EigenMatrix_H
+
+#include "scalarMatrices.H"
+#include "DiagonalMatrix.H"
+#include "SquareMatrix.H"
+#include "complex.H"
+#include <algorithm>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class EigenMatrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class cmptType>
+class EigenMatrix
+{
+    static_assert
+    (
+        std::is_floating_point<cmptType>::value,
+        "EigenMatrix operates only with scalar base type."
+    );
+
+    // Private Data
+
+        //- Number of rows and columns in input matrix
+        const label n_;
+
+        //- Real eigenvalues or real part of complex eigenvalues
+        DiagonalMatrix<cmptType> EValsRe_;
+
+        //- Zero-matrix for real eigenvalues
+        //- or imaginary part of complex eigenvalues
+        DiagonalMatrix<cmptType> EValsIm_;
+
+        //- Right eigenvectors matrix where each column is
+        //- a right eigenvector that corresponds to an eigenvalue
+        SquareMatrix<cmptType> EVecs_;
+
+        //- Copy of nonsymmetric input matrix evolving to eigenvectors matrix
+        SquareMatrix<cmptType> H_;
+
+
+    // Private Member Functions
+
+        //- Householder transform of a symmetric matrix to tri-diagonal form
+        void tridiagonaliseSymmMatrix();
+
+        //- Symmetric tri-diagonal QL algorithm
+        void symmTridiagonalQL();
+
+        //- Reduce non-symmetric matrix to upper-Hessenberg form
+        void Hessenberg();
+
+        //- Reduce matrix from Hessenberg to real Schur form
+        void realSchur();
+
+
+public:
+
+    // Generated Methods
+
+        //- No default construct
+        EigenMatrix() = delete;
+
+        //- No copy construct
+        EigenMatrix(const EigenMatrix&) = delete;
+
+        //- No copy assignment
+        EigenMatrix& operator=(const EigenMatrix&) = delete;
+
+
+    // Constructors
+
+        //- Construct from a SquareMatrix<cmptType>
+        explicit EigenMatrix(const SquareMatrix<cmptType>& A);
+
+        //- Construct from a SquareMatrix<cmptType> and symmetry flag
+        //  Does not perform symmetric check
+        EigenMatrix(const SquareMatrix<cmptType>& A, bool symmetric);
+
+
+    // Access
+
+        //- Return real eigenvalues or real part of complex eigenvalues
+        const DiagonalMatrix<cmptType>& EValsRe() const
+        {
+            return EValsRe_;
+        }
+
+        //- Return zero-matrix for real eigenvalues
+        //- or imaginary part of complex eigenvalues
+        const DiagonalMatrix<cmptType>& EValsIm() const
+        {
+            return EValsIm_;
+        }
+
+        //- Return right eigenvectors matrix where each column is
+        //- a right eigenvector that corresponds to an eigenvalue
+        const SquareMatrix<cmptType>& EVecs() const
+        {
+            return EVecs_;
+        }
+
+        //- Return right eigenvectors in unpacked form
+        const SquareMatrix<complex> complexEVecs() const;
+
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "EigenMatrix.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/Matrix/Matrix.C

@@ -644,10 +644,16 @@ void Foam::Matrix<Form, Type>::operator/=(const Type& s)
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
 
+//- Find max value in Matrix
 template<class Form, class Type>
-const Type& Foam::max(const Matrix<Form, Type>& mat)
+const Type& max(const Matrix<Form, Type>& mat)
 {
     if (mat.empty())
     {
@@ -659,8 +665,9 @@ const Type& Foam::max(const Matrix<Form, Type>& mat)
 }
 
 
+//- Find min value in Matrix
 template<class Form, class Type>
-const Type& Foam::min(const Matrix<Form, Type>& mat)
+const Type& min(const Matrix<Form, Type>& mat)
 {
     if (mat.empty())
     {
@@ -672,8 +679,9 @@ const Type& Foam::min(const Matrix<Form, Type>& mat)
 }
 
 
+//- Find the min/max values of Matrix
 template<class Form, class Type>
-Foam::MinMax<Type> Foam::minMax(const Matrix<Form, Type>& mat)
+MinMax<Type> minMax(const Matrix<Form, Type>& mat)
 {
     MinMax<Type> result;
 
@@ -686,10 +694,12 @@ Foam::MinMax<Type> Foam::minMax(const Matrix<Form, Type>& mat)
 }
 
 
+
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
+//- Matrix negation
 template<class Form, class Type>
-Form Foam::operator-(const Matrix<Form, Type>& mat)
+Form operator-(const Matrix<Form, Type>& mat)
 {
     Form result(mat.sizes());
 
@@ -705,8 +715,9 @@ Form Foam::operator-(const Matrix<Form, Type>& mat)
 }
 
 
+//- Matrix addition. Returns Matrix of the same form as the first parameter.
 template<class Form1, class Form2, class Type>
-Form1 Foam::operator+
+Form1 operator+
 (
     const Matrix<Form1, Type>& A,
     const Matrix<Form2, Type>& B
@@ -738,8 +749,9 @@ Form1 Foam::operator+
 }
 
 
+//- Matrix subtraction. Returns Matrix of the same form as the first parameter.
 template<class Form1, class Form2, class Type>
-Form1 Foam::operator-
+Form1 operator-
 (
     const Matrix<Form1, Type>& A,
     const Matrix<Form2, Type>& B
@@ -771,8 +783,9 @@ Form1 Foam::operator-
 }
 
 
+//- Scalar multiplication of Matrix
 template<class Form, class Type>
-Form Foam::operator*(const Type& s, const Matrix<Form, Type>& mat)
+Form operator*(const Type& s, const Matrix<Form, Type>& mat)
 {
     Form result(mat.sizes());
 
@@ -788,8 +801,17 @@ Form Foam::operator*(const Type& s, const Matrix<Form, Type>& mat)
 }
 
 
+//- Scalar multiplication of Matrix
 template<class Form, class Type>
-Form Foam::operator+(const Type& s, const Matrix<Form, Type>& mat)
+Form operator*(const Matrix<Form, Type>& mat, const Type& s)
+{
+    return s*mat;
+}
+
+
+//- Scalar addition of Matrix
+template<class Form, class Type>
+Form operator+(const Type& s, const Matrix<Form, Type>& mat)
 {
     Form result(mat.sizes());
 
@@ -805,8 +827,17 @@ Form Foam::operator+(const Type& s, const Matrix<Form, Type>& mat)
 }
 
 
+//- Scalar addition of Matrix
 template<class Form, class Type>
-Form Foam::operator-(const Type& s, const Matrix<Form, Type>& mat)
+Form operator+(const Matrix<Form, Type>& mat, const Type& s)
+{
+    return s + mat;
+}
+
+
+//- Scalar subtraction of Matrix
+template<class Form, class Type>
+Form operator-(const Type& s, const Matrix<Form, Type>& mat)
 {
     Form result(mat.sizes());
 
@@ -822,22 +853,9 @@ Form Foam::operator-(const Type& s, const Matrix<Form, Type>& mat)
 }
 
 
+//- Scalar subtraction of Matrix
 template<class Form, class Type>
-Form Foam::operator*(const Matrix<Form, Type>& mat, const Type& s)
+Form operator-(const Matrix<Form, Type>& mat, const Type& s)
-{
-    return s*mat;
-}
-
-
-template<class Form, class Type>
-Form Foam::operator+(const Matrix<Form, Type>& mat, const Type& s)
-{
-    return s + mat;
-}
-
-
-template<class Form, class Type>
-Form Foam::operator-(const Matrix<Form, Type>& mat, const Type& s)
 {
     Form result(mat.sizes());
 
@@ -853,8 +871,9 @@ Form Foam::operator-(const Matrix<Form, Type>& mat, const Type& s)
 }
 
 
+//- Scalar division of Matrix
 template<class Form, class Type>
-Form Foam::operator/(const Matrix<Form, Type>& mat, const Type& s)
+Form operator/(const Matrix<Form, Type>& mat, const Type& s)
 {
     Form result(mat.sizes());
 
@@ -870,9 +889,10 @@ Form Foam::operator/(const Matrix<Form, Type>& mat, const Type& s)
 }
 
 
+//- Matrix-Matrix multiplication using ikj-algorithm
 template<class Form1, class Form2, class Type>
-typename Foam::typeOfInnerProduct<Type, Form1, Form2>::type
+typename typeOfInnerProduct<Type, Form1, Form2>::type
-Foam::operator*
+operator*
 (
     const Matrix<Form1, Type>& A,
     const Matrix<Form2, Type>& B
@@ -912,9 +932,10 @@ Foam::operator*
 }
 
 
+//- Implicit inner product of Matrix-Matrix, equivalent to A.T()*B
 template<class Form1, class Form2, class Type>
-typename Foam::typeOfInnerProduct<Type, Form1, Form2>::type
+typename typeOfInnerProduct<Type, Form1, Form2>::type
-Foam::operator&
+operator&
 (
     const Matrix<Form1, Type>& AT,
     const Matrix<Form2, Type>& B
@@ -954,9 +975,10 @@ Foam::operator&
 }
 
 
+//- Implicit outer product of Matrix-Matrix, equivalent to A*B.T()
 template<class Form1, class Form2, class Type>
-typename Foam::typeOfInnerProduct<Type, Form1, Form2>::type
+typename typeOfInnerProduct<Type, Form1, Form2>::type
-Foam::operator^
+operator^
 (
     const Matrix<Form1, Type>& A,
     const Matrix<Form2, Type>& BT
@@ -996,7 +1018,11 @@ Foam::operator^
 }
 
 
-// * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
 
 #include "MatrixIO.C"
 

src/OpenFOAM/matrices/Matrix/Matrix.H

@@ -220,7 +220,6 @@ public:
         //  Subscript checking only with FULLDEBUG
         inline Type& at(const label idx);
 
-
     // Block Access (const)
 
         //- Return const column or column's subset of Matrix
@@ -265,7 +264,6 @@ public:
             const label colIndex
         ) const;
 
-
     // Block Access (non-const)
 
         //- Return column or column's subset of Matrix
@@ -301,7 +299,6 @@ public:
             const label colIndex
         );
 
-
     // Check
 
         //- Check index i is within valid range [0, m)
@@ -316,7 +313,6 @@ public:
         //- True if all entries have identical values, and Matrix is non-empty
         inline bool uniform() const;
 
-
     // Edit
 
         //- Clear Matrix, i.e. set sizes to zero
@@ -342,10 +338,9 @@ public:
         //- Resize Matrix without reallocating storage (unsafe)
         inline void shallowResize(const label m, const label n);
 
-        //- Round to zero elements with magnitude smaller than tol (SMALL)
+        //- Round elements with magnitude smaller than tol (SMALL) to zero
         void round(const scalar tol = SMALL);
 
-
     // Operations
 
         //- Return (conjugate) transpose of Matrix
@@ -386,7 +381,6 @@ public:
 
         //- Return Frobenius norm of Matrix
         //  Optional without sqrt for parallel usage.
-        //  https://en.wikipedia.org/wiki/Matrix_norm#Frobenius_norm
         scalar norm(const bool noSqrt=false) const;
 
 
@@ -596,159 +590,6 @@ template<class Form, class Type>
 Ostream& operator<<(Ostream& os, const Matrix<Form, Type>& mat);
 
 
-// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
-
-//- Find max value in Matrix
-template<class Form, class Type>
-const Type& max(const Matrix<Form, Type>& mat);
-
-//- Find min value in Matrix
-template<class Form, class Type>
-const Type& min(const Matrix<Form, Type>& mat);
-
-//- Find the min/max values of Matrix
-template<class Form, class Type>
-MinMax<Type> minMax(const Matrix<Form, Type>& mat);
-
-
-// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * * //
-
-//- Matrix negation
-template<class Form, class Type>
-Form operator-(const Matrix<Form, Type>& mat);
-
-//- Matrix addition. Returns Matrix of the same form as the first parameter.
-template<class Form1, class Form2, class Type>
-Form1 operator+
-(
-    const Matrix<Form1, Type>& A,
-    const Matrix<Form2, Type>& B
-);
-
-//- Matrix subtraction. Returns Matrix of the same form as the first parameter.
-template<class Form1, class Form2, class Type>
-Form1 operator-
-(
-    const Matrix<Form1, Type>& A,
-    const Matrix<Form2, Type>& B
-);
-
-//- Scalar multiplication of Matrix
-template<class Form, class Type>
-Form operator*
-(
-    const Type& s,
-    const Matrix<Form, Type>& mat
-);
-
-//- Scalar addition of Matrix
-template<class Form, class Type>
-Form operator+
-(
-    const Type& s,
-    const Matrix<Form, Type>& mat
-);
-
-//- Scalar subtraction of Matrix
-template<class Form, class Type>
-Form operator-
-(
-    const Type& s,
-    const Matrix<Form, Type>& mat
-);
-
-//- Scalar multiplication of Matrix
-template<class Form, class Type>
-Form operator*
-(
-    const Matrix<Form, Type>& mat,
-    const Type& s
-);
-
-//- Scalar addition of Matrix
-template<class Form, class Type>
-Form operator+
-(
-    const Matrix<Form, Type>& mat,
-    const Type& s
-);
-
-//- Scalar subtraction of Matrix
-template<class Form, class Type>
-Form operator-
-(
-    const Matrix<Form, Type>& mat,
-    const Type& s
-);
-
-//- Scalar division of Matrix
-template<class Form, class Type>
-Form operator/
-(
-    const Matrix<Form, Type>& mat,
-    const Type& s
-);
-
-//- Matrix-Matrix multiplication using ikj-algorithm
-template<class Form1, class Form2, class Type>
-typename typeOfInnerProduct<Type, Form1, Form2>::type
-operator*
-(
-    const Matrix<Form1, Type>& A,
-    const Matrix<Form2, Type>& B
-);
-
-//- Matrix-vector multiplication (A * x), where x is a column vector
-template<class Form, class Type>
-inline tmp<Field<Type>> operator*
-(
-    const Matrix<Form, Type>& mat,
-    const UList<Type>& x
-);
-
-//- Matrix-vector multiplication (A * x), where x is a column vector
-template<class Form, class Type, class Addr>
-inline tmp<Field<Type>> operator*
-(
-    const Matrix<Form, Type>& mat,
-    const IndirectListBase<Type, Addr>& x
-);
-
-//- Vector-Matrix multiplication (x * A), where x is a row vector
-template<class Form, class Type>
-inline tmp<Field<Type>> operator*
-(
-    const UList<Type>& x,
-    const Matrix<Form, Type>& mat
-);
-
-//- Vector-Matrix multiplication (x * A), where x is a row vector
-template<class Form, class Type, class Addr>
-inline tmp<Field<Type>> operator*
-(
-    const IndirectListBase<Type, Addr>& x,
-    const Matrix<Form, Type>& mat
-);
-
-//- Implicit inner product of Matrix-Matrix, equivalent to A.T()*B
-template<class Form1, class Form2, class Type>
-typename typeOfInnerProduct<Type, Form1, Form2>::type
-operator&
-(
-    const Matrix<Form1, Type>& ATranspose,
-    const Matrix<Form2, Type>& B
-);
-
-//- Implicit outer product of Matrix-Matrix, equivalent to A*B.T()
-template<class Form1, class Form2, class Type>
-typename typeOfInnerProduct<Type, Form1, Form2>::type
-operator^
-(
-    const Matrix<Form1, Type>& A,
-    const Matrix<Form2, Type>& BTranspose
-);
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/matrices/Matrix/MatrixI.H

@@ -242,7 +242,7 @@ inline const Type& Foam::Matrix<Form, Type>::at(const label idx) const
             << abort(FatalError);
     }
     #endif
-    return (v_ + idx);
+    return *(v_ + idx);
 }
 
 
@@ -257,7 +257,7 @@ inline Type& Foam::Matrix<Form, Type>::at(const label idx)
             << abort(FatalError);
     }
     #endif
-    return (v_ + idx);
+    return *(v_ + idx);
 }
 
 
@@ -608,8 +608,16 @@ inline Type* Foam::Matrix<Form, Type>::operator[](const label irow)
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
+
+//- Matrix-vector multiplication (A * x), where x is a column vector
 template<class Form, class Type>
-inline Foam::tmp<Foam::Field<Type>> Foam::operator*
+inline tmp<Field<Type>> operator*
 (
     const Matrix<Form, Type>& mat,
     const UList<Type>& x
@@ -619,8 +627,9 @@ inline Foam::tmp<Foam::Field<Type>> Foam::operator*
 }
 
 
+//- Matrix-vector multiplication (A * x), where x is a column vector
 template<class Form, class Type, class Addr>
-inline Foam::tmp<Foam::Field<Type>> Foam::operator*
+inline tmp<Field<Type>> operator*
 (
     const Matrix<Form, Type>& mat,
     const IndirectListBase<Type, Addr>& x
@@ -630,8 +639,9 @@ inline Foam::tmp<Foam::Field<Type>> Foam::operator*
 }
 
 
+//- Vector-Matrix multiplication (x * A), where x is a row vector
 template<class Form, class Type>
-inline Foam::tmp<Foam::Field<Type>> Foam::operator*
+inline tmp<Field<Type>> operator*
 (
     const UList<Type>& x,
     const Matrix<Form, Type>& mat
@@ -641,8 +651,9 @@ inline Foam::tmp<Foam::Field<Type>> Foam::operator*
 }
 
 
+//- Vector-Matrix multiplication (x * A), where x is a row vector
 template<class Form, class Type, class Addr>
-inline Foam::tmp<Foam::Field<Type>> Foam::operator*
+inline tmp<Field<Type>> operator*
 (
     const IndirectListBase<Type, Addr>& x,
     const Matrix<Form, Type>& mat
@@ -651,5 +662,8 @@ inline Foam::tmp<Foam::Field<Type>> Foam::operator*
     return mat.Tmul(x);
 }
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
 
 // ************************************************************************* //

src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C

@@ -176,7 +176,7 @@ void Foam::MatrixBlock<MatrixType>::operator=
     const ConstMatrixBlock<MatrixType>& Mb
 )
 {
-    if (this != &Mb)
+    if (reinterpret_cast<const ConstMatrixBlock<MatrixType>*>(this) != &Mb)
     {
         if (mRows_ != Mb.m() || nCols_ != Mb.n())
         {
@@ -233,7 +233,7 @@ void Foam::MatrixBlock<MatrixType>::operator=
     const ConstMatrixBlock<MatrixType2>& Mb
 )
 {
-    if (this != &Mb)
+    if (reinterpret_cast<const ConstMatrixBlock<MatrixType2>*>(this) != &Mb)
     {
         if (mRows_ != Mb.m() || nCols_ != Mb.n())
         {

src/OpenFOAM/matrices/QRMatrix/QRMatrix.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -497,15 +497,17 @@ MatrixType Foam::pinv
 {
     typedef typename MatrixType::cmptType cmptType;
 
-    #if FULLDEBUG
     if (A.empty())
     {
         FatalErrorInFunction
-            << "Empty matrix." << abort(FatalError);
+            << "Empty matrix found."
+            << abort(FatalError);
     }
 
-    // Check if A is full-rank
+    if (A.size() == 1)
-    #endif
+    {
+        return MatrixType({1, 1}, cmptType(1.0)/(A(0,0) + cmptType(VSMALL)));
+    }
 
     QRMatrix<MatrixType> QRM
     (
@@ -534,9 +536,12 @@ MatrixType Foam::pinv
 
     if (firstZeroElemi == 0)
     {
-        FatalErrorInFunction
+        WarningInFunction
             << "The largest (magnitude) diagonal element is (almost) zero."
-            << abort(FatalError);
+            << nl << "Returning a zero matrix."
+            << endl;
+
+        return MatrixType(A.sizes(), Zero);
     }
 
     // Exclude the first (almost) zero diagonal row and the rows below

src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,13 +28,14 @@ Class
     Foam::RectangularMatrix
 
 Description
-    A templated 2D rectangular m x n matrix of objects of \<Type\>.
+    A templated (M x N) rectangular matrix of objects of \<Type\>,
+    containing M*N elements, derived from Matrix.
 
-    The matrix dimensions are used for subscript bounds checking etc.
+See also
+    Test-RectangularMatrix.C
 
 SourceFiles
     RectangularMatrixI.H
-    RectangularMatrix.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -61,10 +62,19 @@ class RectangularMatrix
 
 public:
 
-    // Constructors
+    // Generated Methods
 
-        //- Construct null
+        //- Default construct
-        inline RectangularMatrix();
+        RectangularMatrix() = default;
+
+        //- Copy construct
+        RectangularMatrix(const RectangularMatrix&) = default;
+
+        //- Copy assignment
+        RectangularMatrix& operator=(const RectangularMatrix&) = default;
+
+
+    // Constructors
 
         //- Construct a square matrix (rows == columns)
         inline explicit RectangularMatrix(const label n);
@@ -85,15 +95,15 @@ public:
         template<class AnyType>
         inline RectangularMatrix(const labelPair& dims, const Identity<AnyType>);
 
-        //- Construct given number of rows/columns
+        //- Construct given number of rows/columns by using a label pair
         inline explicit RectangularMatrix(const labelPair& dims);
 
-        //- Construct given number of rows/columns
+        //- Construct given number of rows/columns by using a label pair
-        //- initializing all elements to zero
+        //- and initializing all elements to zero
         inline RectangularMatrix(const labelPair& dims, const zero);
 
-        //- Construct given number of rows/columns
+        //- Construct given number of rows/columns by using a label pair
-        //- initializing all elements to the given value
+        //- and initializing all elements to the given value
         inline RectangularMatrix(const labelPair& dims, const Type& val);
 
         //- Construct from a block of another matrix
@@ -107,7 +117,7 @@ public:
         //- Construct as copy of a square matrix
         inline RectangularMatrix(const SquareMatrix<Type>& mat);
 
-        //- Construct from Istream.
+        //- Construct from Istream
         inline explicit RectangularMatrix(Istream& is);
 
         //- Clone
@@ -124,41 +134,6 @@ public:
 };
 
 
-// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-class typeOfInnerProduct<Type, RectangularMatrix<Type>, RectangularMatrix<Type>>
-{
-public:
-
-    typedef RectangularMatrix<Type> type;
-};
-
-
-template<class Type>
-class typeOfInnerProduct<Type, RectangularMatrix<Type>, SquareMatrix<Type>>
-{
-public:
-
-    typedef RectangularMatrix<Type> type;
-};
-
-
-template<class Type>
-class typeOfInnerProduct<Type, SquareMatrix<Type>, RectangularMatrix<Type>>
-{
-public:
-
-    typedef RectangularMatrix<Type> type;
-};
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-RectangularMatrix<Type> outer(const Field<Type>& f1, const Field<Type>& f2);
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,13 +28,6 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class Type>
-inline Foam::RectangularMatrix<Type>::RectangularMatrix()
-:
-    Matrix<RectangularMatrix<Type>, Type>()
-{}
-
-
 template<class Type>
 inline Foam::RectangularMatrix<Type>::RectangularMatrix
 (
@@ -197,8 +190,38 @@ inline void Foam::RectangularMatrix<Type>::operator=(const Type& val)
 namespace Foam
 {
 
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+class typeOfInnerProduct<Type, RectangularMatrix<Type>, RectangularMatrix<Type>>
+{
+public:
+
+    typedef RectangularMatrix<Type> type;
+};
+
+
+template<class Type>
+class typeOfInnerProduct<Type, RectangularMatrix<Type>, SquareMatrix<Type>>
+{
+public:
+
+    typedef RectangularMatrix<Type> type;
+};
+
+
+template<class Type>
+class typeOfInnerProduct<Type, SquareMatrix<Type>, RectangularMatrix<Type>>
+{
+public:
+
+    typedef RectangularMatrix<Type> type;
+};
+
+
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 
+// Return the outer product of Field-Field as RectangularMatrix
 template<class Type>
 inline Foam::RectangularMatrix<Type> outer
 (
@@ -224,5 +247,4 @@ inline Foam::RectangularMatrix<Type> outer
 
 } // End namespace Foam
 
-
 // ************************************************************************* //

src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,47 +30,52 @@ License
 #include "RectangularMatrix.H"
 #include "labelList.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-Foam::scalar Foam::detDecomposed
+template<class CompOp>
+Foam::List<Foam::label> Foam::SquareMatrix<Type>::sortPermutation
 (
-    const SquareMatrix<Type>& matrix,
+    CompOp& compare
-    const label sign
+) const
-)
 {
-    Type diagProduct = pTraits<Type>::one;
+    List<label> p(this->m());
+    std::iota(p.begin(), p.end(), 0);
+    std::sort
+    (
+        p.begin(),
+        p.end(),
+        [&](label i, label j){ return compare((*this)(i,i), (*this)(j,j)); }
+    );
 
-    for (label i = 0; i < matrix.m(); ++i)
+    return p;
-    {
-        diagProduct *= matrix(i, i);
-    }
-
-    return sign*diagProduct;
 }
 
 
 template<class Type>
-Foam::scalar Foam::det(const SquareMatrix<Type>& matrix)
+void Foam::SquareMatrix<Type>::applyPermutation(const List<label>& p)
 {
-    SquareMatrix<Type> matrixTmp = matrix;
+    #ifdef FULLDEBUG
+    if (this->m() != p.size())
+    {
+        FatalErrorInFunction
+            << "Attempt to column-reorder according to an uneven list: " << nl
+            << "SquareMatrix diagonal size = " << this->m() << nl
+            << "Permutation list size = " << p.size() << nl
+            << abort(FatalError);
+    }
+    #endif
 
-    labelList pivotIndices(matrix.m());
+    SquareMatrix<Type> reordered(this->sizes());
-    label sign;
-    LUDecompose(matrixTmp, pivotIndices, sign);
 
-    return detDecomposed(matrixTmp, sign);
+    label j = 0;
+    for (const label i : p)
+    {
+        reordered.subColumn(j) = this->subColumn(i);
+        ++j;
     }
 
-
+    this->transfer(reordered);
-template<class Type>
-Foam::scalar Foam::det(SquareMatrix<Type>& matrix)
-{
-    labelList pivotIndices(matrix.m());
-    label sign;
-    LUDecompose(matrix, pivotIndices, sign);
-
-    return detDecomposed(matrix, sign);
 }
 
 
@@ -90,3 +95,103 @@ void Foam::SquareMatrix<Type>::operator=(const Identity<AnyType>)
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+//- Return the determinant of the LU decomposed SquareMatrix
+template<class Type>
+scalar detDecomposed
+(
+    const SquareMatrix<Type>& matrix,
+    const label sign
+)
+{
+    Type diagProduct = pTraits<Type>::one;
+
+    for (label i = 0; i < matrix.m(); ++i)
+    {
+        diagProduct *= matrix(i, i);
+    }
+
+    return sign*diagProduct;
+}
+
+
+//- Return the determinant of SquareMatrix
+template<class Type>
+scalar det(const SquareMatrix<Type>& matrix)
+{
+    SquareMatrix<Type> matrixTmp = matrix;
+
+    labelList pivotIndices(matrix.m());
+    label sign;
+    LUDecompose(matrixTmp, pivotIndices, sign);
+
+    return detDecomposed(matrixTmp, sign);
+}
+
+
+//- Return the SquareMatrix det and the LU decomposition in the original matrix
+template<class Type>
+scalar det(SquareMatrix<Type>& matrix)
+{
+    labelList pivotIndices(matrix.m());
+    label sign;
+    LUDecompose(matrix, pivotIndices, sign);
+
+    return detDecomposed(matrix, sign);
+}
+
+
+//- Return Matrix column-reordered according to
+//- a given permutation labelList
+template<class Type>
+SquareMatrix<Type> applyPermutation
+(
+    const SquareMatrix<Type>& mat,
+    const List<label>& p
+)
+{
+    #ifdef FULLDEBUG
+    if (mat.m() != p.size())
+    {
+        FatalErrorInFunction
+            << "Attempt to column-reorder according to an uneven list: " << nl
+            << "SquareMatrix diagonal size = " << mat.m() << nl
+            << "Permutation list size = " << p.size() << nl
+            << abort(FatalError);
+    }
+    #endif
+
+    SquareMatrix<Type> reordered(mat.sizes());
+
+    label j = 0;
+    for (const label i : p)
+    {
+        reordered.subColumn(j) = mat.subColumn(i);
+        ++j;
+    }
+
+    return reordered;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+class typeOfInnerProduct<Type, SquareMatrix<Type>, SquareMatrix<Type>>
+{
+public:
+
+    typedef SquareMatrix<Type> type;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,8 +28,11 @@ Class
     Foam::SquareMatrix
 
 Description
-    A templated 2D square matrix of objects of \<T\>, where the n x n matrix
+    A templated (N x N) square matrix of objects of \<Type\>,
-    dimension is known and used for subscript bounds checking, etc.
+    containing N*N elements, derived from Matrix.
+
+See also
+    Test-SquareMatrix.C
 
 SourceFiles
     SquareMatrixI.H
@@ -42,6 +45,7 @@ SourceFiles
 
 #include "Matrix.H"
 #include "Identity.H"
+#include <numeric>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,10 +68,19 @@ class SquareMatrix
 
 public:
 
-    // Constructors
+    // Generated Methods
 
-        //- Construct null
+        //- Default construct
-        inline SquareMatrix();
+        SquareMatrix() = default;
+
+        //- Copy construct
+        SquareMatrix(const SquareMatrix&) = default;
+
+        //- Copy assignment
+        SquareMatrix& operator=(const SquareMatrix&) = default;
+
+
+    // Constructors
 
         //- Construct for given size (rows == cols)
         inline explicit SquareMatrix(const label n);
@@ -85,20 +98,25 @@ public:
         template<class AnyType>
         inline SquareMatrix(const label n, const Identity<AnyType>);
 
+        //- Construct for given size (rows == cols) by using a labelPair
+        //- initializing to the identity matrix
         template<class AnyType>
         inline SquareMatrix(const labelPair& dims, const Identity<AnyType>);
 
-        //- Construct given number of rows/columns (checked to be equal)
+        //- Construct given number of rows/columns
+        //- by using a labelPair (checked to be equal)
         //  For constructor consistency in rectangular matrices
         inline explicit SquareMatrix(const labelPair& dims);
 
-        //- Construct given number of rows/columns (checked to be equal)
+        //- Construct given number of rows/columns
-        //- initializing all elements to zero
+        //- by using a labelPair (checked to be equal)
+        //- and initializing all elements to zero
         //  For constructor consistency with rectangular matrices
         inline SquareMatrix(const labelPair& dims, const zero);
 
-        //- Construct given number of rows/columns (checked to be equal)
+        //- Construct given number of rows/columns
-        //- initializing all elements to the given value
+        //- by using a labelPair (checked to be equal)
+        //- and initializing all elements to the given value
         //  For constructor consistency with rectangular matrices
         inline SquareMatrix(const labelPair& dims, const Type& val);
 
@@ -106,7 +124,7 @@ public:
         //- initializing all elements to zero
         inline SquareMatrix(const label m, const label n, const zero);
 
-        //- Construct from sub-matrix block
+        //- Construct from const sub-matrix block
         template<class MatrixType>
         inline SquareMatrix(const ConstMatrixBlock<MatrixType>& mat);
 
@@ -149,6 +167,17 @@ public:
         //- Return true if the square matrix is reduced tridiagonal
         inline bool tridiagonal() const;
 
+    // Sort
+
+        //- Return a sort permutation labelList according to
+        //- a given comparison on the diagonal entries
+        template<class CompOp>
+        List<label> sortPermutation(CompOp& compare) const;
+
+        //- Column-reorder this Matrix according to
+        //- a given permutation labelList
+        void applyPermutation(const List<label>& p);
+
 
     // Member Operators
 
@@ -164,28 +193,6 @@ public:
 };
 
 
-// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
-
-//- Return the LU decomposed SquareMatrix det
-template<class Type>
-scalar detDecomposed(const SquareMatrix<Type>&, const label sign);
-
-//- Return the SquareMatrix det
-template<class Type>
-scalar det(const SquareMatrix<Type>&);
-
-//- Return the SquareMatrix det and the LU decomposition in the original matrix
-template<class Type>
-scalar det(SquareMatrix<Type>&);
-
-template<class Type>
-class typeOfInnerProduct<Type, SquareMatrix<Type>, SquareMatrix<Type>>
-{
-public:
-
-    typedef SquareMatrix<Type> type;
-};
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,13 +37,6 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class Type>
-inline Foam::SquareMatrix<Type>::SquareMatrix()
-:
-    Matrix<SquareMatrix<Type>, Type>()
-{}
-
-
 template<class Type>
 inline Foam::SquareMatrix<Type>::SquareMatrix(const label n)
 :
@@ -318,6 +311,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 
+// Return the outer product of Field-Field as SquareMatrix
 template<class Type>
 inline Foam::SquareMatrix<Type> symmOuter
 (

src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -126,4 +126,4 @@ void Foam::SymmetricSquareMatrix<Type>::operator=(const Identity<AnyType>)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// ************************************************************************* //

src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,8 +28,11 @@ Class
     Foam::SymmetricSquareMatrix
 
 Description
-    A templated 2D square symmetric matrix of objects of \<T\>, where the
+    A templated (N x N) square matrix of objects of \<Type\>,
-    n x n matrix dimension is known and used for subscript bounds checking, etc.
+    containing N*N elements, derived from Matrix.
+
+See also
+    Test-SymmetricSquareMatrix.C
 
 SourceFiles
     SymmetricSquareMatrixI.H
@@ -59,10 +62,19 @@ class SymmetricSquareMatrix
 
 public:
 
-    // Constructors
+    // Generated Methods
 
-        //- Construct null
+        //- Default construct
-        inline SymmetricSquareMatrix();
+        SymmetricSquareMatrix() = default;
+
+        //- Copy construct
+        SymmetricSquareMatrix(const SymmetricSquareMatrix&) = default;
+
+        //- Copy assignment
+        SymmetricSquareMatrix& operator=(const SymmetricSquareMatrix&) = default;
+
+
+    // Constructors
 
         //- Construct for given size (rows == cols)
         inline explicit SymmetricSquareMatrix(const label n);
@@ -76,7 +88,7 @@ public:
         inline SymmetricSquareMatrix(const label n, const Type& val);
 
         //- Construct for given size (rows == cols)
-        //- setting the diagonal to one and off-diagonals to zero
+        //- initializing to the identity matrix
         template<class AnyType>
         inline SymmetricSquareMatrix(const label n, const Identity<AnyType>);
 

src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,13 +37,6 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class Type>
-inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix()
-:
-    Matrix<SymmetricSquareMatrix<Type>, Type>()
-{}
-
-
 template<class Type>
 inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix(const label n)
 :

src/functionObjects/field/Make/files

@@ -118,4 +118,6 @@ stabilityBlendingFactor/stabilityBlendingFactor.C
 
 interfaceHeight/interfaceHeight.C
 
+STDMD/STDMD.C
+
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects

src/functionObjects/field/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/OpenFOAM/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
@@ -22,6 +23,7 @@ EXE_INC = \
     -I$(LIB_SRC)/phaseSystemModels/reactingEulerFoam/phaseSystems/lnInclude
 
 LIB_LIBS = \
+    -lOpenFOAM \
     -lfiniteVolume \
     -lfileFormats \
     -lsurfMesh \

src/functionObjects/field/STDMD/STDMD.H

@@ -0,0 +1,644 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::STDMD
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    (Beta Release) STDMD (i.e. Streaming Total Dynamic Mode Decomposition) is
+    a variant of a data-driven dimensionality reduction method.
+
+    STDMD is being used as a mathematical post-processing tool to compute
+    a set of dominant modes out of a given flow (or dataset) each of which is
+    associated with a constant frequency and decay rate, so that dynamic
+    features of a given flow may become interpretable, and tractable.
+    Among other Dynamic Mode Decomposition (DMD) variants, STDMD is presumed
+    to provide the general DMD method capabilities alongside economised and
+    feasible memory and CPU usage.
+
+    The code implementation corresponds to Figs. 15-16 of the first citation
+    below, more broadly to Section 2.4.
+
+    References:
+    \verbatim
+    STDMD and mode-sorting algorithms (tags:K, HRDC, KZ, HWR):
+        Kiewat, M. (2019).
+        Streaming modal decomposition approaches for vehicle aerodynamics.
+        PhD thesis. Munich: Technical University of Munich.
+        URL:mediatum.ub.tum.de/doc/1482652/1482652.pdf
+
+        Hemati, M. S., Rowley, C. W., Deem, E. A., & Cattafesta, L. N. (2017).
+        De-biasing the dynamic mode decomposition for applied Koopman
+        spectral analysis of noisy datasets.
+        Theoretical and Computational Fluid Dynamics, 31(4), 349-368.
+        DOI:10.1007/s00162-017-0432-2
+
+        Kou, J., & Zhang, W. (2017).
+        An improved criterion to select dominant modes from dynamic mode
+        decomposition.
+        European Journal of Mechanics-B/Fluids, 62, 109-129.
+        DOI:10.1016/j.euromechflu.2016.11.015
+
+        Hemati, M. S., Williams, M. O., & Rowley, C. W. (2014).
+        Dynamic mode decomposition for large and streaming datasets.
+        Physics of Fluids, 26(11), 111701.
+        DOI:10.1063/1.4901016
+
+    Parallel classical Gram-Schmidt process (tag:Ka):
+        Katagiri, T. (2003).
+        Performance evaluation of parallel Gram-Schmidt re-orthogonalization
+        methods.
+        In: Palma J. M. L. M., Sousa A. A., Dongarra J., Hernández V. (eds)
+        High Performance Computing for Computational Science — VECPAR 2002.
+        Lecture Notes in Computer Science, vol 2565, p. 302-314.
+        Berlin, Heidelberg: Springer.
+        DOI:10.1007/3-540-36569-9_19
+
+    Parallel direct tall-skinny QR decomposition (tags:BGD, DGHL):
+        Benson, A. R., Gleich, D. F., & Demmel, J. (2013).
+        Direct QR factorizations for tall-and-skinny matrices in MapReduce
+        architectures.
+        2013 IEEE International Conference on Big Data.
+        DOI:10.1109/bigdata.2013.6691583
+
+        Demmel, J., Grigori, L., Hoemmen, M., & Langou, J. (2012).
+        Communication-optimal parallel and sequential QR and LU
+        factorizations.
+        SIAM Journal on Scientific Computing, 34(1), A206-A239.
+        DOI:10.1137/080731992
+
+    DMD properties:
+        Brunton S. L. (2018).
+        Dynamic mode decomposition overview.
+        Seattle, Washington: University of Washington.
+        youtu.be/sQvrK8AGCAo (Retrieved:24-04-20)
+    \endverbatim
+
+    Operands:
+    \table
+      Operand      | Type              | Location
+      input        | {vol,surface}\<Type\>Field <!--
+               --> | $FOAM_CASE/\<time\>/\<inpField\>
+      output file  | dat                        <!--
+               --> | $FOAM_CASE/postProcessing/\<FO\>/\<time\>/\<file\>(s)
+      output field | volScalarField(s) | $FOAM_CASE/\<time\>/\<outField\>(s)
+    \endtable
+
+    where Type={Scalar,SphericalTensor,SymmTensor,Tensor}.
+
+    Output fields:
+    \verbatim
+      modeReal<modeIndex><field> | Real part of a mode field
+      modeImag<modeIndex><field> | Imaginary part of a mode field
+    \endverbatim
+
+    Output files:
+    \verbatim
+      uSTDMD.dat                 | Unfiltered STDMD output
+      STDMD.dat                  | Filtered STDMD output
+    \endverbatim
+
+    wherein for each mode, the following quantities are output into files:
+    \verbatim
+      Frequency
+      Magnitude
+      Amplitude (real part)
+      Amplitude (imaginary part)
+      Eigenvalue (real part)
+      Eigenvalue (imaginary part)
+    \endverbatim
+
+Note
+    Operations on boundary fields, e.g. \c wallShearStress, are currently not
+    available.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    STDMD1
+    {
+        // Mandatory entries (unmodifiable)
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               <inpField>;
+
+        // Conditional mandatory entries (unmodifiable)
+        // either of the below
+        stdmdInterval       5.5;
+        executeInterval     10;
+
+        // Optional entries (unmodifiable)
+        modeSorter          kiewat;
+        nModes              50;
+        maxRank             50;
+        nGramSchmidt        5;
+        fMin                0;
+        fMax                1000000000;
+
+        // Optional entries (run-time modifiable, yet not recommended)
+        testEigen           false;
+        dumpEigen           false;
+        minBasis            0.00000001;
+        minEVal             0.00000001;
+        absTol              0.001;
+        relTol              0.0001;
+        sortLimiter         500;
+
+        // Optional (inherited) entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Req'd | Dflt
+      type         | Type name: STDMD                   | word |  yes  | -
+      libs         | Library name: fieldFunctionObjects | word |  yes  | -
+      field        | Name of the operand field          | word |  yes  | -
+      stdmdInterval | STDMD time-step size [s]   | scalar| conditional | <!--
+                  --> executeInterval*(current time-step of the simulation)
+      modeSorter   | Mode-sorting algorithm variant     | word | no  | firstSnap
+      nModes | Number of output modes in input freq range | label | no | GREAT
+      maxRank  | Max columns in accumulated matrices      | label | no | GREAT
+      nGramSchmidt | Number of Gram-Schmidt iterations    | label | no | 5
+      fMin         | Min (non-negative) output frequency  | label | no | 0
+      fMax         | Max output frequency                 | label | no | GREAT
+      testEigen    | Flag to verify eigendecompositions   | bool  | no | false
+      dumpEigen    | Flag to log operands of eigendecomps | bool  | no | false
+      minBasis | Orthogonal basis expansion threshold     | scalar| no | 1e-8
+      minEVal  | Min EVal for below EVals are omitted     | scalar| no | 1e-8
+      absTol   | Min abs tol in eigendecomposition tests  | scalar| no | 1e-4
+      relTol   | Relative tol in eigendecomposition tests | scalar| no | 1e-6
+      sortLimiter | Maximum allowable magnitude for <!--
+               -->  mag(eigenvalue)*(number of STDMD steps) to be used in <!--
+               -->  modeSorter={kiewat,kouZhang} to avoid overflow errors <!--
+               -->                                        | scalar | no | 500.0
+    \endtable
+
+    Options for the \c modeSorter entry:
+    \verbatim
+      kiewat    | Modified weighted amplitude scaling method
+      kouZhang  | Weighted amplitude scaling method
+      firstSnap | Method of first snapshot amplitude magnitude
+    \endverbatim
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+     - \link writeFile.H \endlink
+
+    Usage by \c postProcess utility is not available.
+
+Note
+    - To specify the STDMD time-step size (not necessarily equal to the
+    time step of the simulation), entries of either \c stdmdInterval or
+    \c executeInterval must be available (highest to lowest precedence). While
+    \c stdmdInterval allows users to directly specify the STDMD time-step size
+    in seconds, in absence of \c stdmdInterval, for convenience,
+    \c executeInterval allows users to compute the STDMD time-step internally
+    by multiplying itself with the current time-step size of the simulation.
+    - Limitation: Restart is currently not available since intermediate writing
+    of STDMD matrices are not supported.
+    - Limitation: Non-physical input (e.g. full-zero fields) may upset STDMD.
+    - Warning: DMD is an active research area at the time of writing; therefore,
+    there would be cases whereat oddities might be encountered.
+    - Warning: This STDMD release is the \b beta release; therefore,
+    small-to-medium changes in input/output interfaces and internal structures
+    should be expected in the next versions.
+    - Warning: In \c modeSorter algorithms of \c kiewat and \c kouZhang, there
+    exists a function of \f$power(x,y)\f$ where \c x is the magnitude of an
+    eigenvalue, and \c y is the number of STDMD snapshots. This function poses
+    a risk of overflow errors since, for example, if \c x ends up above 1.5 with
+    500 snapshots or more, this function automatically throws the floating point
+    error meaning overflow. Therefore, the domain-range of this function is
+    heuristically constrained by the optional entry \c sortLimiter which skips
+    the evaluation of this function for a given
+    mag(eigenvalue)*(no. of STDMD steps), i.e. x*y, whose magnitude is larger
+    than \c sortLimiter.
+
+See also
+    - Foam::functionObjects::fvMeshFunctionObject
+    - Foam::functionObjects::writeFile
+    - ExtendedCodeGuide::functionObjects::field::STDMD
+
+SourceFiles
+    STDMD.C
+    STDMDTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_STDMD_H
+#define functionObjects_STDMD_H
+
+#include "fvMeshFunctionObject.H"
+#include "writeFile.H"
+#include "EigenMatrix.H"
+#include "QRMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                              Class STDMD Declaration
+\*---------------------------------------------------------------------------*/
+
+class STDMD
+:
+    public fvMeshFunctionObject,
+    public writeFile
+{
+
+    typedef SquareMatrix<scalar> SMatrix;
+    typedef SquareMatrix<complex> SCMatrix;
+    typedef RectangularMatrix<scalar> RMatrix;
+    typedef RectangularMatrix<complex> RCMatrix;
+
+
+    // Private Enumerations
+
+        //- Options for the mode-sorting algorithm
+        enum modeSorterType
+        {
+            KIEWAT,          //!< "Modified weighted amplitude scaling method"
+            KOU_ZHANG,       //!< "Weighted amplitude scaling method"
+            FIRST_SNAPSHOT   //!< "Method of first snapshot amplitude magnitude"
+        };
+
+        //- Names for modeSorterType
+        static const Enum<modeSorterType> modeSorterTypeNames;
+
+
+    // Private Data
+
+        //- Mode-sorting algorithm (default=modeSorterType::FIRST_SNAPSHOT)
+        const enum modeSorterType modeSorter_;
+
+        //- Name of the operand volume or surface field
+        const word fieldName_;
+
+        //- Flag: First execution-step initialisation
+        bool initialised_;
+
+        //- Flag: Verify eigendecompositions by using theoretical relations
+        bool testEigen_;
+
+        //- Flag: Write operands of eigendecompositions to log
+        //  To activate it, testEigen=true
+        bool dumpEigen_;
+
+        //- Number of output modes within input frequency range
+        //- starting from the most energetic mode
+        const label nModes_;
+
+        //- Maximum allowable matrix column for Qz_ and Gz_
+        //  Qz_ is assumed to always have full-rank, thus Qz_.n() = rank
+        const label maxRank_;
+
+        //- Number of maximum iterations of the classical Gram-Schmidt procedure
+        const label nGramSchmidt_;
+
+        //- Min frequency: Output only entries corresponding to freqs >= fMin
+        const label fMin_;
+
+        //- Max frequency: Output only entries corresponding to freqs <= fMax
+        const label fMax_;
+
+        //- Number of base components of input field, e.g. 3 for vector
+        label nComps_;
+
+        //- Number of elements in a snapshot
+        //  A snapshot is an input dataset to be processed per execution-step
+        label nSnap_;
+
+        //- Current execution-step index of STDMD,
+        //- not necessarily that of the simulation
+        label currIndex_;
+
+        //- Maximum allowable magnitude for mag(eigenvalue)*(number of
+        //- STDMD steps) to be used in modeSorter={kiewat,kouZhang} to avoid
+        //- overflow errors
+        scalar sortLimiter_;
+
+        //- Min absolute tolerance being used in eigendecomposition tests
+        scalar absTol_;
+
+        //- Relative tolerance being used in eigendecomposition tests
+        scalar relTol_;
+
+        //- Min value to execute orthogonal basis expansion of Qz_ and Gz_
+        scalar minBasis_;
+
+        //- STDMD time-step size that equals to
+        //- (executeInterval of STDMD)*(deltaT of simulation) [s]
+        scalar dt_;
+
+        //- Min EVal magnitude threshold below where EVals are omitted
+        scalar minMagEVal_;
+
+        //- L2-norm of column vector z_
+        scalar zNorm_;
+
+        //- L2-norm of column vector ez_
+        scalar ezNorm_;
+
+        //- Augmented snapshot matrix (effectively a column vector) (K:Eq. 60)
+        //  Upper half z_ = current-time snapshot slot
+        //  Lower half z_ = previous-time snapshot slot
+        RMatrix z_;
+
+        //- Working copy of the augmented snapshot matrix z_
+        //- being used in the classical Gram-Schmidt process
+        RMatrix ez_;
+
+        //- First-processed snapshot required by the mode-sorting
+        //- algorithms at the final output computations (K:p. 43)
+        RMatrix X1_;
+
+        //- Accumulated-in-time unitary similarity matrix produced by the
+        //- orthonormal decomposition of the augmented snapshot matrix z_
+        //  (K:Eq. 60)
+        RMatrix Qz_;
+
+        //- Accumulated-in-time (squared) upper triangular matrix produced by
+        //- the orthonormal decomposition of the augmented snapshot matrix z_
+        //  (K:Eq. 64)
+        SMatrix Gz_;
+
+        //- Upper half of Qz_
+        RMatrix QzUH_;
+
+        //- Lower half of Qz_
+        RMatrix QzLH_;
+
+        //- Moore-Penrose pseudo-inverse of R produced by
+        //- the QR decomposition of the last time-step QzUH_
+        RMatrix RxInv_;
+
+        //- Projected STDMD operator (K:Eq. 78)
+        SMatrix Ap_;
+
+        //- Output eigenvectors
+        RCMatrix oEVecs_;
+
+        //- Output eigenvalues
+        List<complex> oEVals_;
+
+        //- Output amplitudes
+        List<complex> oAmps_;
+
+        //- Output (non-negative) frequencies
+        List<scalar> oFreqs_;
+
+        //- Indices of oFreqs_ where freqs are
+        //- non-negative and within [fMin_, fMax_]
+        DynamicList<label> iFreqs_;
+
+        //- Output (descending) magnitudes of (complex) amplitudes
+        List<scalar> oMags_;
+
+        //- Indices of oMags_
+        List<label> iMags_;
+
+
+    // Private Member Functions
+
+        // Process
+
+        //- Move previous-time snapshot into previous-time slot in z_
+        //- and copy new current-time snapshot into current-time slot in z_
+        template<class GeoFieldType>
+        bool getSnapshot();
+
+        //- Get the input field type to be processed by snapshot()
+        template<class Type>
+        bool getSnapshotType();
+
+        //- Get the number of base components of input field
+        template<class GeoFieldType>
+        bool getComps();
+
+        //- Return (parallel) L2-norm of a given column vector
+        scalar parnorm(const RMatrix& colVector) const;
+
+        //- Move the current-time snapshot into the previous-time snapshot in z_
+        //- and copy the new field into the current-time snapshot
+        void snapshot();
+
+        //- Initialise all working matrices at the first execution-step
+        void init();
+
+        //- Initialise Qz_, Gz_ (both require the first two snapshots) and X1_
+        void initBasis();
+
+        //- Execute (parallel) classical Gram-Schmidt
+        //- process to orthonormalise ez_ (Ka:Fig. 5)
+        void GramSchmidt();
+
+        //- Expand orthonormal bases Qz_ and Gz_ by stacking a column
+        //- (ez_/ezNorm_) to Qz_, and a row (Zero) and column (Zero)
+        //- to Gz_ if (minBasis_ < ezNorm_/zNorm_)
+        void expandBasis();
+
+        //- Update Gz_ before the potential orthonormal basis compression
+        void updateGz();
+
+        //- Compress orthonormal basis for Qz_ and Gz_ if (maxRank_ < Qz_.n())
+        void compressBasis();
+
+
+        // Postprocess
+
+        //- Return a new List containing elems of List at indices
+        template<class Type>
+        void filterIndexed
+        (
+            List<Type>& lst,
+            const UList<label>& indices
+        );
+
+        //- Return a new Matrix containing columns of Matrix at indices
+        template<class MatrixType>
+        void filterIndexed
+        (
+            MatrixType& lst,
+            const UList<label>& indices
+        );
+
+        //- Compute global Ap (K:Eq. 78)
+        void calcAp();
+
+        //- Compute eigenvalues and eigenvectors
+        void calcEigen();
+
+        //- Weak-type floating-point comparison
+        //  bit.ly/2Trdbgs (Eq. 2), and PEP-485
+        bool close
+        (
+            const scalar s1,
+            const scalar s2,
+            const scalar absTol = 0,   //<! comparisons near zero
+            const scalar relTol = 1e-8 //<! e.g. vals the same within 8 decimals
+        ) const;
+
+        //- Test real/complex eigenvalues by using
+        //- the theoretical relation: (sum(eigenvalues) - trace(A) ~ 0)
+        void testEigenvalues
+        (
+            const SquareMatrix<scalar>& A,
+            const DiagonalMatrix<scalar>& EValsRe
+        ) const;
+
+        //- Test real eigenvectors by using the theoretical relation:
+        //- ((A & EVec - EVal*EVec).norm() ~ 0)
+        void testEigenvectors
+        (
+            const SquareMatrix<scalar>& A,
+            const DiagonalMatrix<scalar>& EValsRe,
+            const SquareMatrix<scalar>& EVecs
+        ) const;
+
+        //- Test complex eigenvectors by using the theoretical relation:
+        //- ((A & EVec - EVal*EVec).norm() ~ 0)
+        void testEigenvectors
+        (
+            const SquareMatrix<scalar>& A,
+            const List<complex>& EVals,
+            const RectangularMatrix<complex>& EVecs
+        ) const;
+
+        //- Remove any eigenvalues whose magnitude is smaller than
+        //- minMagEVal_ while keeping the order of elements the same
+        void filterEVals();
+
+        //- Compute and filter frequencies (K:Eq. 81)
+        void calcFreqs();
+
+        //- Compute frequency indices
+        //  Locate indices where oFreqs_ are
+        //  in [fMin_, fMax_], and store them in iFreqs_ indices
+        void calcFreqI();
+
+        //- Compute amplitudes
+        void calcAmps();
+
+        //- Compute magnitudes
+        //  Includes advanced sorting methods using
+        //  the chosen weighted amplitude scaling method
+        void calcMags();
+
+        //- Compute the ordered magnitude indices
+        void calcMagI();
+
+        //- Compute modes
+        void calcModes();
+
+        //- Eigenvalue weighted amplitude scaling (KZ:Eq. 33)
+        //- Modified eigenvalue weighted amplitude scaling (K)
+        scalar sorter
+        (
+            const List<scalar>& weight,
+            const complex& amplitude,
+            const complex& eval,
+            const scalar modeNorm
+        ) const;
+
+        //- Output file header information
+        virtual void writeFileHeader(Ostream& os) const;
+
+        //- Filter objects according to iFreqs_ and iMags_
+        void filterOutput();
+
+        //- Write unfiltered/filtered data
+        void writeOutput(OFstream& os) const;
+
+        //- Compute STDMD output
+        void calcOutput();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("STDMD");
+
+
+    // Constructors
+
+        //- No default construct
+        STDMD() = delete;
+
+        //- Construct from Time and dictionary
+        STDMD
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+        //- No copy construct
+        STDMD(const STDMD&) = delete;
+
+        //- No copy assignment
+        void operator=(const STDMD&) = delete;
+
+
+    //- Destructor
+    virtual ~STDMD() = default;
+
+
+    // Member Functions
+
+        //- Read STDMD settings
+        virtual bool read(const dictionary&);
+
+        //- Execute STDMD
+        virtual bool execute();
+
+        //- Write STDMD output
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "STDMDTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/functionObjects/field/STDMD/STDMDTemplates.C

@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class GeoFieldType>
+bool Foam::functionObjects::STDMD::getSnapshot()
+{
+    if (!initialised_)
+    {
+        init();
+    }
+
+    // Move previous-time snapshot into previous-time slot in z_
+    // Effectively moves the lower half of z_ to its upper half
+    std::rotate(z_.begin(), z_.begin() + nSnap_, z_.end());
+
+    // Copy new current-time snapshot into current-time slot in z_
+    // Effectively copies the new field elements into the lower half of z_
+    const GeoFieldType& Field = lookupObject<GeoFieldType>(fieldName_);
+    const label nField = Field.size();
+
+    for (direction dir = 0; dir < nComps_; ++dir)
+    {
+        z_.subColumn(0, nSnap_ + dir*nField, nField) = Field.component(dir);
+    }
+
+    return true;
+}
+
+
+template<class Type>
+bool Foam::functionObjects::STDMD::getSnapshotType()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+
+    if (foundObject<VolFieldType>(fieldName_))
+    {
+        return getSnapshot<VolFieldType>();
+    }
+    else if (foundObject<SurfaceFieldType>(fieldName_))
+    {
+        return getSnapshot<SurfaceFieldType>();
+    }
+
+    return false;
+}
+
+
+template<class Type>
+bool Foam::functionObjects::STDMD::getComps()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+
+    if (foundObject<VolFieldType>(fieldName_))
+    {
+        nComps_ = pTraits<typename VolFieldType::value_type>::nComponents;
+        return true;
+    }
+    else if (foundObject<SurfaceFieldType>(fieldName_))
+    {
+        nComps_ = pTraits<typename SurfaceFieldType::value_type>::nComponents;
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class Type>
+void Foam::functionObjects::STDMD::filterIndexed
+(
+    List<Type>& lst,
+    const UList<label>& indices
+)
+{
+    // Elems within [a, b]
+    List<Type> lstWithin(indices.size());
+
+    // Copy if frequency of elem is within [a, b]
+    label j = 0;
+    for (const auto& i : indices)
+    {
+        lstWithin[j] = lst[i];
+        ++j;
+    }
+    lst.transfer(lstWithin);
+}
+
+
+template<class MatrixType>
+void Foam::functionObjects::STDMD::filterIndexed
+(
+    MatrixType& mat,
+    const UList<label>& indices
+)
+{
+    // Elems within [a, b]
+    MatrixType matWithin(labelPair(mat.m(), indices.size()));
+
+    // Copy if frequency of elem is within [a, b]
+    label j = 0;
+    for (const auto& i : indices)
+    {
+        matWithin.subColumn(j) = mat.subColumn(i);
+        ++j;
+    }
+    mat.transfer(matWithin);
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/0.orig/U

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0.012984 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    topAndBottom
+    {
+        type            symmetry;
+    }
+
+    "left|right"
+    {
+        type            empty;
+    }
+
+    cylinder
+    {
+        type            noSlip;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/0.orig/p

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    topAndBottom
+    {
+        type            symmetry;
+    }
+
+    "left|right"
+    {
+        type            empty;
+    }
+
+    cylinder
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allclean

@@ -0,0 +1,12 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase0
+
+rm -rf dynamicCode
+rm -rf constant/coarseMesh
+rm -f system/coarseMesh/decomposeParDict
+
+#------------------------------------------------------------------------------

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allrun

@@ -0,0 +1,25 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+./Allrun.pre
+
+runApplication decomposePar
+
+runParallel renumberMesh -overwrite
+
+cp -f system/decomposeParDict system/coarseMesh
+
+runApplication -s coarseMesh decomposePar -region coarseMesh
+
+runParallel -s coarseMesh renumberMesh -overwrite -region coarseMesh
+
+runParallel $(getApplication)
+
+runParallel -s main redistributePar -reconstruct -latestTime
+
+runParallel -s coarseMesh redistributePar -reconstruct \
+    -region coarseMesh -time '50,100,200'
+
+#------------------------------------------------------------------------------

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/Allrun.pre

@@ -0,0 +1,20 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+runApplication -s coarseMesh blockMesh -dict system/blockMeshDict.coarse
+
+runApplication snappyHexMesh -overwrite
+
+mkdir -p constant/coarseMesh
+
+mv -f constant/polyMesh constant/coarseMesh
+
+runApplication -s main blockMesh -dict system/blockMeshDict.main
+
+runApplication mirrorMesh -overwrite
+
+restore0Dir
+
+#------------------------------------------------------------------------------

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/constant/transportProperties

@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              1.558e-5;
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/constant/turbulenceProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.coarse

@@ -0,0 +1,74 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale    1;
+
+x0       -0.16;
+x1       0.8;
+y0       -0.24;
+y1       0.24;
+z0       -0.00375;
+z1       0.00375;
+
+nx       90;
+ny       60;
+nz       1;
+
+vertices
+(
+    ($x0 $y0 $z0)
+    ($x1 $y0 $z0)
+    ($x1 $y1 $z0)
+    ($x0 $y1 $z0)
+    ($x0 $y0 $z1)
+    ($x1 $y0 $z1)
+    ($x1 $y1 $z1)
+    ($x0 $y1 $z1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    allPatches
+    {
+        type     empty;
+        inGroups (allPatches);
+        faces
+        (
+            (3 7 6 2)
+            (1 5 4 0)
+            (0 4 7 3)
+            (2 6 5 1)
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
+

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.main

@@ -0,0 +1,277 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale     1;
+
+R         0.06; // Cylinder radius
+x0        #eval{ 0.5*$R };
+x2        #eval{ 2.38732*$R };
+x3        #eval{ 10.0*$R };
+xOutlet   #eval{ 18.6667*$R };
+xInlet    #eval{ -10.125*$R };
+
+RsinPi8   #eval{ $R*sin(0.125*pi()) };
+RsinPi8n  #eval{ -$RsinPi8 };
+RcosPi8   #eval{ $R*cos(0.125*pi()) };
+RcosPi8n  #eval{ -$RcosPi8 };
+RsinPi4   #eval{ $R*sin(0.25*pi()) };
+
+x2sinPi8  #eval{ $x2*sin(0.125*pi()) };
+x2sinPi8n #eval{ -$x2sinPi8 };
+x2cosPi8  #eval{ $x2*cos(0.125*pi()) };
+x2cosPi8n #eval{ -$x2cosPi8 };
+x2sinPi4  #eval{ $x2*sin(0.25*pi()) };
+
+z0        -0.0075;
+z1        0.0075;
+nz        1;
+
+
+vertices #codeStream
+{
+    codeInclude
+    #{
+        #include "pointField.H"
+    #};
+
+    code
+    #{
+        pointField points(19);
+        points[0] = point($R, 0, $z0);
+        points[1] = point($x2, 0, $z0);
+        points[2] = point($x3, 0, $z0);
+        points[3] = point($x3, $x2sinPi4, $z0);
+        points[4] = point($x2sinPi4, $x2sinPi4, $z0);
+        points[5] = point($RsinPi4, $RsinPi4, $z0);
+        points[6] = point($x3, $x3, $z0);
+        points[7] = point($x2sinPi4, $x3, $z0);
+
+        // Mirror +x points to -x side
+        points[11] = point(-points[0].x(), points[0].y(), points[0].z());
+        points[12] = point(-points[1].x(), points[1].y(), points[1].z());
+        points[13] = point(-points[2].x(), points[2].y(), points[2].z());
+        points[14] = point(-points[3].x(), points[3].y(), points[3].z());
+        points[15] = point(-points[4].x(), points[4].y(), points[4].z());
+        points[16] = point(-points[5].x(), points[5].y(), points[5].z());
+        points[17] = point(-points[6].x(), points[6].y(), points[6].z());
+        points[18] = point(-points[7].x(), points[7].y(), points[7].z());
+
+        // Points on the y-z-plane
+        points[8] = point(0, $x3, $z0);
+        points[9] = point(0, $x2, $z0);
+        points[10] = point(0, $R, $z0);
+
+        // Mirror -z points to +z side
+        label sz = points.size();
+        points.setSize(2*sz);
+        for (label i = 0; i < sz; ++i)
+        {
+            const point& pt = points[i];
+            points[i + sz] = point(pt.x(), pt.y(), $z1);
+        }
+
+        // Add an inner cylinder
+        sz = points.size();
+        label nAdd = 6;
+        points.setSize(sz + nAdd);
+
+        // Points within the inner cylinder
+        points[sz] = point(0, 0, $z0);
+        points[sz + 1] = point(0, $x0, $z0);
+        points[sz + 2] = point($x0, $x0, $z0);
+        points[sz + 3] = point($x0, 0, $z0);
+
+        // Mirror points from +x side to -x side
+        points[sz + 4] =
+            point(-points[sz + 2].x(), points[sz + 2].y(), points[sz + 2].z());
+        points[sz + 5] =
+            point(-points[sz + 3].x(), points[sz + 3].y(), points[sz + 3].z());
+
+        // Mirror -z points to +z side
+        sz = points.size();
+        points.setSize(sz + nAdd);
+        for (label i = 0; i < nAdd; ++i)
+        {
+            const point& pt = points[i+sz-nAdd];
+            points[i+sz] = point(pt.x(), pt.y(), $z1);
+        }
+
+        // Add downstream and upstream blocks
+        sz = points.size();
+        nAdd = 6;
+        points.setSize(sz + nAdd);
+
+        // Points on outlet
+        points[sz] = point($xOutlet, 0, $z0);
+        points[sz + 1] = point($xOutlet, $x3, $z0);
+        points[sz + 4] = point($xOutlet, $x2sinPi4, $z0);
+
+        // Points on inlet
+        points[sz + 2] = point($xInlet, 0, $z0);
+        points[sz + 3] = point($xInlet, $x3, $z0);
+        points[sz + 5] = point($xInlet, $x2sinPi4, $z0);
+
+        // Mirror -z points to +z side
+        sz = points.size();
+        points.setSize(sz + nAdd);
+        for (label i = 0; i < nAdd; ++i)
+        {
+            const point& pt = points[i + sz - nAdd];
+            points[i + sz] = point(pt.x(), pt.y(), $z1);
+        }
+
+        os  << points;
+    #};
+};
+
+
+blocks
+(
+    hex ( 5  4  9 10 24 23 28 29) (16 15 $nz) simpleGrading (1.6 1 1)
+    hex ( 0  1  4  5 19 20 23 24) (16 15 $nz) simpleGrading (1.6 1 1)
+    hex ( 1  2  3  4 20 21 22 23) (61 15 $nz) simpleGrading (1 1 1)
+    hex ( 4  3  6  7 23 22 25 26) (61 61 $nz) simpleGrading (1 1 1)
+    hex ( 9  4  7  8 28 23 26 27) (15 61 $nz) simpleGrading (1 1 1)
+    hex (15 16 10  9 34 35 29 28) (16 15 $nz) simpleGrading (0.625 1 1)
+    hex (12 11 16 15 31 30 35 34) (16 15 $nz) simpleGrading (0.625 1 1)
+    hex (13 12 15 14 32 31 34 33) (61 15 $nz) simpleGrading (1 1 1)
+    hex (14 15 18 17 33 34 37 36) (61 61 $nz) simpleGrading (1 1 1)
+    hex (15  9  8 18 34 28 27 37) (15 61 $nz) simpleGrading (1 1 1)
+    hex (2 50 54 3 21 56 60 22)  (69 15 $nz) simpleGrading (1 1 1) // downstream
+    hex (3 54 51 6 22 60 57 25)  (69 61 $nz) simpleGrading (1 1 1)
+    hex (52 13 14 55 58 32 33 61)  (1 15 $nz) simpleGrading (1 1 1) // upstream
+    hex (55 14 17 53 61 33 36 59)  (1 61 $nz) simpleGrading (1 1 1)
+);
+
+edges
+(
+    arc  0  5 ($RcosPi8   $RsinPi8  $z0)
+    arc  5 10 ($RsinPi8   $RcosPi8  $z0)
+    arc  1  4 ($x2cosPi8  $x2sinPi8 $z0)
+    arc  4  9 ($x2sinPi8  $x2cosPi8 $z0)
+    arc 19 24 ($RcosPi8   $RsinPi8  $z1)
+    arc 24 29 ($RsinPi8   $RcosPi8  $z1)
+    arc 20 23 ($x2cosPi8  $x2sinPi8 $z1)
+    arc 23 28 ($x2sinPi8  $x2cosPi8 $z1)
+    arc 11 16 ($RcosPi8n  $RsinPi8  $z0)
+    arc 16 10 ($RsinPi8n  $RcosPi8  $z0)
+    arc 12 15 ($x2cosPi8n $x2sinPi8 $z0)
+    arc 15  9 ($x2sinPi8n $x2cosPi8 $z0)
+    arc 30 35 ($RcosPi8n  $RsinPi8  $z1)
+    arc 35 29 ($RsinPi8n  $RcosPi8  $z1)
+    arc 31 34 ($x2cosPi8n $x2sinPi8 $z1)
+    arc 34 28 ($x2sinPi8n $x2cosPi8 $z1)
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (53 59 61 55)
+            (55 61 58 52)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (60 57 51 54)
+            (56 60 54 50)
+        );
+    }
+    topAndBottom
+    {
+        type symmetry;
+        faces
+        (
+            (59 53 17 36)
+            (36 17 18 37)
+            (37 18  8 27)
+            (27  8  7 26)
+            (26  7  6 25)
+            (25  6 51 57)
+            (52 58 32 13)
+            (13 32 31 12)
+            (12 31 30 11)
+            (0 19 20 1)
+            (1 20 21 2)
+            (2 21 56 50)
+        );
+    }
+    left
+    {
+        type empty;
+        faces
+        (
+            (55 52 13 14)
+            (53 55 14 17)
+            (13 12 15 14)
+            (14 15 18 17)
+            (15 9 8 18)
+            (12 11 16 15)
+            (15 16 10 9)
+            (9 10 5 4)
+            (5 0 1 4)
+            (9 4 7 8)
+            (4 3 6 7)
+            (1 2 3 4)
+            (3 2 50 54)
+            (6 3 54 51)
+        );
+    }
+    right
+    {
+        type empty;
+        faces
+        (
+            (59 36 33 61)
+            (61 33 32 58)
+            (36 37 34 33)
+            (33 34 31 32)
+            (37 27 28 34)
+            (27 26 23 28)
+            (26 25 22 23)
+            (23 22 21 20)
+            (25 57 60 22)
+            (22 60 56 21)
+            (24 23 20 19)
+            (28 23 24 29)
+            (34 28 29 35)
+            (31 34 35 30)
+        );
+    }
+    cylinder
+    {
+        type wall;
+        faces
+        (
+            (29 24 5 10)
+            (24 19 0 5)
+            (30 35 16 11)
+            (35 29 10 16)
+        );
+    }
+);
+
+mergePatchPairs
+();
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/coarseMesh/fvSchemes

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+}
+
+gradSchemes
+{
+}
+
+divSchemes
+{
+}
+
+laplacianSchemes
+{
+}
+
+interpolationSchemes
+{
+}
+
+snGradSchemes
+{
+}
+
+fluxRequired
+{
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/coarseMesh/fvSolution

@@ -0,0 +1,27 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+}
+
+PIMPLE
+{
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/controlDict

@@ -0,0 +1,229 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     pimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         200;
+
+deltaT          0.05;
+
+writeControl    runTime;
+
+writeInterval   50;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  16;
+
+timeFormat      general;
+
+timePrecision   8;
+
+runTimeModifiable false;
+
+
+functions
+{
+    mapFields1
+    {
+        type            mapFields;
+        libs            (fieldFunctionObjects);
+        mapRegion       coarseMesh;
+        mapMethod       cellVolumeWeight;
+        consistent      no;
+        patchMap        ();
+        cuttingPatches  ();
+        fields          ( U p );
+
+        timeStart       10;
+        timeEnd         2000;
+
+        executeControl  timeStep;
+        executeInterval 10;
+
+        writeControl    timeStep;
+        writeInterval   50;
+    }
+
+    STDMD1U
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               U;
+
+        // Optional entries
+        modeSorter          firstSnap;
+
+        // Optional (inherited) entries
+        region              coarseMesh;  // mapFields must be executed before.
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD1p
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               p;
+
+        // Optional entries
+        modeSorter          firstSnap;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD2U
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               U;
+
+        // Optional entries
+        modeSorter          firstSnap;
+        nModes              2;
+        maxRank             2;
+        nGramSchmidt        5;
+        fMin                0;
+        fMax                1000000000;
+        testEigen           true;
+        dumpEigen           true;
+        minBasis            0.00000001;
+        minEVal             0.00000001;
+        absTol              0.001;
+        relTol              0.0001;
+
+        // Optional (inherited) entries
+        timeStart           50;
+        timeEnd             2000;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD3U
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               U;
+
+        // Optional entries
+        modeSorter          firstSnap;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        timeEnd             150;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        runTime;
+        writeInterval       50;
+    }
+
+    STDMD4U
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               U;
+
+        // Optional entries
+        modeSorter          kiewat;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD4p
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               p;
+
+        // Optional entries
+        modeSorter          kiewat;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD5U
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               U;
+
+        // Optional entries
+        modeSorter          kouZhang;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+
+    STDMD5p
+    {
+        // Mandatory entries
+        type                STDMD;
+        libs                (fieldFunctionObjects);
+        field               p;
+
+        // Optional entries
+        modeSorter          kouZhang;
+
+        // Optional (inherited) entries
+        region              coarseMesh;
+        timeStart           10;
+        executeControl      timeStep;
+        executeInterval     10;
+        writeControl        onEnd;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/decomposeParDict

@@ -0,0 +1,28 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  8;
+
+method              hierarchical;
+
+coeffs
+{
+    n               (8 1 1);
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/fvSchemes

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/fvSolution

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver           GAMG;
+        smoother         DICGaussSeidel;
+        tolerance        1e-6;
+        relTol           0.01;
+    }
+
+    pFinal
+    {
+        $p;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    UFinal
+    {
+        $U;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    nCorrectors         2;
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/mirrorMeshDict

@@ -0,0 +1,27 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      mirrorMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+planeType           pointAndNormal;
+
+pointAndNormalDict
+{
+    point   (0 0 0);
+    normal  (0 -1 0);
+}
+
+planeTolerance      1e-3;
+
+// ************************************************************************* //

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/snappyHexMeshDict

@@ -0,0 +1,83 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      snappyHexMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+castellatedMesh true;
+snap            false;
+addLayers       false;
+
+geometry
+{
+    cylinder
+    {
+        type    searchableCylinder;
+        point1  (0 0 -0.00375);
+        point2  (0 0 0.00375);
+        radius  0.06;
+    }
+}
+
+
+castellatedMeshControls
+{
+    maxLocalCells       100000;
+    maxGlobalCells      100000000;
+    minRefinementCells  10;
+    maxLoadUnbalance    0.10;
+    nCellsBetweenLevels 1;
+    resolveFeatureAngle 30;
+    allowFreeStandingZoneFaces true;
+    features
+    ();
+    refinementSurfaces
+    {
+        cylinder
+        {
+            level (0 0);
+            patchInfo
+            {
+                type     empty;
+                inGroups (allPatches);
+            }
+        }
+    };
+    refinementRegions{};
+    locationInMesh (-0.1 -0.2 0);
+}
+
+
+snapControls
+{
+    nSmoothPatch    3;
+    tolerance       2.0;
+    nSolveIter      100;
+    nRelaxIter      5;
+}
+
+
+addLayersControls
+{}
+
+
+meshQualityControls
+{
+    #includeEtc "caseDicts/meshQualityDict"
+}
+
+debug 0;
+mergeTolerance 1e-6;
+
+
+// ************************************************************************* //