ENH: improve Matrix classes and tests

applications/test/TestTools/TestTools.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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

+Test-DiagonalMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-DiagonalMatrix

applications/test/matrices/DiagonalMatrix/Make/options

+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */

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

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/Matrix/Test-Matrix.C → 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/QRMatrix/Test-QRMatrix.C → 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

+Test-RectangularMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-RectangularMatrix

applications/test/matrices/RectangularMatrix/Make/options

+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */

applications/test/matrices/SquareMatrix/Make/files

+Test-SquareMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-SquareMatrix

applications/test/matrices/SquareMatrix/Make/options

+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */

applications/test/matrices/SymmetricSquareMatrix/Make/files

+Test-SymmetricSquareMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-SymmetricSquareMatrix