diff --git a/applications/test/TestTools/TestTools.H b/applications/test/TestTools/TestTools.H
new file mode 100644
index 0000000000000000000000000000000000000000..c24cb9fa051f69335395ec93ebe31049c8f81fab
--- /dev/null
+++ b/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_;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/matrices/DiagonalMatrix/Make/files b/applications/test/matrices/DiagonalMatrix/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..0f5877472a346b75694df94f6ff6b489c99cabdc
--- /dev/null
+++ b/applications/test/matrices/DiagonalMatrix/Make/files
@@ -0,0 +1,3 @@
+Test-DiagonalMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-DiagonalMatrix
diff --git a/applications/test/matrices/DiagonalMatrix/Make/options b/applications/test/matrices/DiagonalMatrix/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..b23590707c8d302fe75ce411d5bf921c029a7ff3
--- /dev/null
+++ b/applications/test/matrices/DiagonalMatrix/Make/options
@@ -0,0 +1,2 @@
+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */
diff --git a/applications/test/matrices/DiagonalMatrix/Test-DiagonalMatrix.C b/applications/test/matrices/DiagonalMatrix/Test-DiagonalMatrix.C
new file mode 100644
index 0000000000000000000000000000000000000000..3d8ae06979f579c8977c1fa3d216b42bc1034f46
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/Matrix/Make/files b/applications/test/matrices/Matrix/Make/files
similarity index 100%
rename from applications/test/Matrix/Make/files
rename to applications/test/matrices/Matrix/Make/files
diff --git a/applications/test/Matrix/Make/options b/applications/test/matrices/Matrix/Make/options
similarity index 100%
rename from applications/test/Matrix/Make/options
rename to applications/test/matrices/Matrix/Make/options
diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/matrices/Matrix/Test-Matrix.C
similarity index 99%
rename from applications/test/Matrix/Test-Matrix.C
rename to applications/test/matrices/Matrix/Test-Matrix.C
index 1cbc60254538df916b33a17d12bb5a579795fe18..43115df871371067f2ce87d5f1cf4a1905511c15 100644
--- a/applications/test/Matrix/Test-Matrix.C
+++ b/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.
diff --git a/applications/test/QRMatrix/Make/files b/applications/test/matrices/QRMatrix/Make/files
similarity index 100%
rename from applications/test/QRMatrix/Make/files
rename to applications/test/matrices/QRMatrix/Make/files
diff --git a/applications/test/QRMatrix/Make/options b/applications/test/matrices/QRMatrix/Make/options
similarity index 100%
rename from applications/test/QRMatrix/Make/options
rename to applications/test/matrices/QRMatrix/Make/options
diff --git a/applications/test/QRMatrix/Test-QRMatrix.C b/applications/test/matrices/QRMatrix/Test-QRMatrix.C
similarity index 99%
rename from applications/test/QRMatrix/Test-QRMatrix.C
rename to applications/test/matrices/QRMatrix/Test-QRMatrix.C
index 281c62dd5b40dcd3d8983541fab42af31fc9608d..679072715257733728458b9f2810891f0fa04053 100644
--- a/applications/test/QRMatrix/Test-QRMatrix.C
+++ b/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()
diff --git a/applications/test/matrices/RectangularMatrix/Make/files b/applications/test/matrices/RectangularMatrix/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..49b520d07986ab159a292f32f6f9455e5e12371e
--- /dev/null
+++ b/applications/test/matrices/RectangularMatrix/Make/files
@@ -0,0 +1,3 @@
+Test-RectangularMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-RectangularMatrix
diff --git a/applications/test/matrices/RectangularMatrix/Make/options b/applications/test/matrices/RectangularMatrix/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..b23590707c8d302fe75ce411d5bf921c029a7ff3
--- /dev/null
+++ b/applications/test/matrices/RectangularMatrix/Make/options
@@ -0,0 +1,2 @@
+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */
diff --git a/applications/test/matrices/RectangularMatrix/Test-RectangularMatrix.C b/applications/test/matrices/RectangularMatrix/Test-RectangularMatrix.C
new file mode 100644
index 0000000000000000000000000000000000000000..fbcf4c6343d72b70b7df298a22c2f41d9d9520c3
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/matrices/SquareMatrix/Make/files b/applications/test/matrices/SquareMatrix/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..e64b858cd6ba23bcc369939af81137b364927bdf
--- /dev/null
+++ b/applications/test/matrices/SquareMatrix/Make/files
@@ -0,0 +1,3 @@
+Test-SquareMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-SquareMatrix
diff --git a/applications/test/matrices/SquareMatrix/Make/options b/applications/test/matrices/SquareMatrix/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..b23590707c8d302fe75ce411d5bf921c029a7ff3
--- /dev/null
+++ b/applications/test/matrices/SquareMatrix/Make/options
@@ -0,0 +1,2 @@
+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */
diff --git a/applications/test/matrices/SquareMatrix/Test-SquareMatrix.C b/applications/test/matrices/SquareMatrix/Test-SquareMatrix.C
new file mode 100644
index 0000000000000000000000000000000000000000..d114dd39639e2ce3cc98a4a018d2ed3c3d6f6cd8
--- /dev/null
+++ b/applications/test/matrices/SquareMatrix/Test-SquareMatrix.C
@@ -0,0 +1,1000 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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-SquareMatrix
+
+Description
+ Tests for \c SquareMatrix 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 "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+#include "IOmanip.H"
+#include "TestTools.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Create each constructor of SquareMatrix<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+ {
+ Info<< "# Construct for given size (rows == cols):" << nl;
+ const SquareMatrix<Type> A(2);
+ Info<< A << endl;
+ }
+ {
+ Info<< "# Construct for given size (rows == cols) "
+ << "initializing all elements to zero:" << nl;
+ const SquareMatrix<Type> A(2, Zero);
+ Info<< A << endl;
+ }
+ {
+ Info<< "# Construct for given size (rows == cols) "
+ << "initializing all elements to a value:" << nl;
+ const SquareMatrix<Type> A(2, Type(3));
+ Info<< A << endl;
+ }
+ {
+ Info<< "# Construct for given size (rows == cols) "
+ << "initializing to the identity matrix:" << nl;
+ const SquareMatrix<Type> A(2, I);
+ Info<< A << endl;
+ }
+ {
+ Info<< "# Construct for given size (rows == cols) by using a labelPair"
+ << "initializing to the identity matrix:" << nl;
+ const SquareMatrix<Type> A(labelPair(2, 2), I);
+ Info<< A << endl;
+ }
+ {
+ Info<< "# Construct given number of rows/columns "
+ << "by using a label pair (checked to be equal):" << nl;
+ const SquareMatrix<Type> A(labelPair(2, 2));
+ Info<< A << endl;
+ }
+ {
+ Info<< "# Construct given number of rows/columns "
+ << "by using a label pair (checked to be equal) "
+ << "and initializing all elements to zero:" << nl;
+ const SquareMatrix<Type> A(labelPair(2, 2), Zero);
+ Info<< A << endl;
+ }
+ {
+ Info<< "# Construct given number of rows/columns "
+ << "by using a label pair (checked to be equal) "
+ << "and initializing all elements to the given value:" << nl;
+ const SquareMatrix<Type> A(labelPair(2, 2), Type(3));
+ Info<< A << endl;
+ }
+ {
+ Info<< "# Construct given number of rows/columns "
+ << "(checked to be equal) initializing all elements to zero" << nl;
+ const SquareMatrix<Type> A(2, 2, Zero);
+ Info<< A << endl;
+ }
+ {
+ Info<< "# Construct from const sub-matrix block:" << nl;
+ const SquareMatrix<Type> B(labelPair(3, 3), Type(3));
+ const SquareMatrix<Type> A(B.subMatrix(1, 1));
+ Info<< A << endl;
+ }
+ {
+ Info<< "# Construct from sub-matrix block:" << nl;
+ SquareMatrix<Type> B(labelPair(3, 3), Type(3));
+ const SquareMatrix<Type> A(B.subMatrix(1, 1));
+ Info<< A << endl;
+ }
+ {
+ Info<< "#: Construct as copy of a RectangularMatrix:" << nl;
+ const RectangularMatrix<Type> B(2, 2, Zero);
+ const SquareMatrix<Type> A(B);
+ Info<< A << endl;
+ }
+
+ // ## Pending ##
+ // Construct from Istream
+ // Clone
+ // #############
+}
+
+
+// Execute each member function of SquareMatrix<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+ SquareMatrix<Type> A(labelPair(3, 3), Zero);
+ assignMatrix
+ (
+ A,
+ {
+ Type(1), Type(-2.2), Type(-3.4),
+ Type(-0.35), Type(1), Type(5),
+ Type(-0.35), Type(1), Type(5)
+ }
+ );
+
+ Info<< "# Operand: " << nl
+ << " SquareMatrix = " << A << endl;
+
+
+ {
+ Info<< "# Access:" << nl;
+
+ cmp(" The number of rows = ", Type(A.m()), Type(3));
+ cmp(" The number of columns = ", Type(A.n()), Type(3));
+ cmp(" The number of elements in Matrix = ", Type(A.size()), Type(9));
+ cmp(" Return row/column sizes = ", A.sizes(), labelPair(3, 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), 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), 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), 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), 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);
+
+ {
+ Random rndGen(1234);
+ const label mRows = 10;
+ SquareMatrix<Type> B
+ (
+ makeRandomMatrix<SquareMatrix<Type>>(labelPair(mRows, mRows), rndGen)
+ );
+ // Symmetrise
+ for (label n = 0; n < B.n() - 1; ++n)
+ {
+ for (label m = B.m() - 1; n < m; --m)
+ {
+ B(n, m) = B(m, n);
+ }
+ }
+ cmp
+ (
+ " Return true if the square matrix is "
+ "effectively symmetric/Hermitian",
+ B.symmetric(),
+ true
+ );
+ }
+ // ## Pending ##
+ // Return true if the square matrix is reduced tridiagonal
+ // #############
+ }
+ {
+ Info<< "# Edit:" << nl;
+
+ SquareMatrix<Type> cpA0(A);
+ cpA0.clear();
+ cmp(" Clear Matrix, i.e. set sizes to zero = ", cpA0.empty(), true);
+
+ SquareMatrix<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),
+ Type(-0.35), Type(1), Type(5)
+ })
+ );
+
+ SquareMatrix<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(2, 2);
+ cmp
+ (
+ " Change Matrix dimensions, preserving the elements = ",
+ flt(cpA2),
+ List<Type>({Type(1), Type(-2.2), Type(-0.35), Type(1)})
+ );
+
+ cpA2.resize(1);
+ cmp
+ (
+ " Resize the matrix preserving the elements = ",
+ flt(cpA2),
+ List<Type>({Type(1)})
+ );
+
+ SquareMatrix<Type> cpA3(A);
+ cpA3.setSize(2);
+ cmp
+ (
+ " Resize the matrix preserving the elements = ",
+ flt(cpA3),
+ List<Type>({Type(1), Type(-2.2), Type(-0.35), Type(1)})
+ );
+
+ SquareMatrix<Type> cpA4(A);
+ cpA4.shallowResize(2);
+ cmp
+ (
+ " Resize the matrix without reallocating storage (unsafe) = ",
+ flt(cpA4),
+ List<Type>({Type(1), Type(-2.2), Type(-3.4), Type(-0.35)})
+ );
+
+ SquareMatrix<Type> smallA(2, Type(VSMALL));
+ smallA.round();
+ cmp
+ (
+ " Round elements with mag smaller than tol (SMALL) to zero = ",
+ flt(smallA),
+ List<Type>(4, Zero)
+ );
+ }
+ {
+ Info<< "# Sort:" << nl;
+
+ SquareMatrix<Type> cpA(A);
+ auto descend = [&](Type a, Type b){ return mag(a) > mag(b); };
+ const List<label> sortPermutation(cpA.sortPermutation(descend));
+ cmp
+ (
+ " Return a sort permutation labelList according to "
+ "a given comparison on the diagonal entries",
+ sortPermutation,
+ List<label>({2, 0, 1})
+ );
+
+ SquareMatrix<Type> sortedA0(3, Zero);
+ assignMatrix
+ (
+ sortedA0,
+ {
+ Type(-3.4), Type(1), Type(-2.2),
+ Type(5), Type(-0.35), Type(1),
+ Type(5), Type(-0.35), Type(1)
+ }
+ );
+ const SquareMatrix<Type> sortedA1(applyPermutation(A, sortPermutation));
+ cmp
+ (
+ " Return Matrix column-reordered according to "
+ "a given permutation labelList",
+ flt(sortedA0),
+ flt(sortedA1)
+ );
+
+ cpA.applyPermutation(sortPermutation);
+ cmp
+ (
+ " Column-reorder this Matrix according to "
+ "a given permutation labelList",
+ flt(sortedA0),
+ flt(cpA)
+ );
+ }
+ {
+ Info<< "# Operations:" << nl;
+
+ cmp(" Transpose = ", flt((A.T()).T()), flt(A));
+
+ SquareMatrix<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)
+ }
+ );
+ SquareMatrix<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));
+
+ SquareMatrix<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 SquareMatrix<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 SquareMatrix<Type>, and print output
+template<class Type>
+void test_member_opers(Type)
+{
+ SquareMatrix<Type> A(3, I);
+ const SquareMatrix<Type> cpA(A);
+
+ Info<< "# Operand: " << nl
+ << " SquareMatrix = " << A << endl;
+
+
+ A = Zero;
+ cmp(" Assign all elements to zero =", flt(A), List<Type>(9, Zero));
+ A = cpA;
+
+ A = Type(5);
+ cmp(" Assign all elements to value =", flt(A), List<Type>(9, Type(5)));
+ A = cpA;
+
+ A = I;
+ cmp(" Set to identity matrix =", flt(A), flt(cpA));
+ 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(SquareMatrix<Type>(3, Zero))
+ );
+ A = cpA;
+
+ A = Type(-1.2);
+ cmp
+ (
+ " Assignment of all elements to the given value = ",
+ flt(A),
+ flt(SquareMatrix<Type>(3, 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(SquareMatrix<Type>(3, Zero))
+ );
+ A = cpA;
+
+ A = Zero;
+ A += Type(5);
+ cmp
+ (
+ " Matrix scalar addition = ",
+ flt(A),
+ flt(SquareMatrix<Type>(3, Type(5)))
+ );
+ A = cpA;
+
+ A = Zero;
+ A -= Type(5);
+ cmp
+ (
+ " Matrix scalar subtraction = ",
+ flt(A),
+ flt(SquareMatrix<Type>(3, Type(-5)))
+ );
+ A = cpA;
+
+ A = Type(1);
+ A *= Type(5);
+ cmp
+ (
+ " Matrix scalar multiplication = ",
+ flt(A),
+ flt(SquareMatrix<Type>(3, Type(5)))
+ );
+ A = cpA;
+
+ A = Type(1);
+ A /= Type(5);
+ cmp
+ (
+ " Matrix scalar division = ",
+ flt(A),
+ flt(SquareMatrix<Type>(3, 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(2,2) = 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 SquareMatrix<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;
+
+ {
+ SquareMatrix<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
+ );
+ cmp
+ (
+ " Return the outer product of Field-Field as SquareMatrix = ",
+ flt(symmOuter(F, F)),
+ flt(A),
+ 1e-6
+ );
+ }
+}
+
+
+// Execute each global function of SquareMatrix<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+ SquareMatrix<Type> A(3, Zero);
+ assignMatrix
+ (
+ A,
+ {
+ Type(1), Type(-2.2), Type(-3.4),
+ Type(-0.35), Type(10.32), Type(5),
+ Type(-0.35), Type(10.32), Type(-5)
+ }
+ );
+ const SquareMatrix<Type> cpA(A);
+
+ Info<< "# Operand: " << nl
+ << " SquareMatrix = " << A << endl;
+
+ // ## Pending ##: 'LUDecompose' does not operate with 'complex' base type
+ // cmp
+ // (
+ // " Return the determinant of SquareMatrix",
+ // det(A),
+ // -95.5
+ // );
+ // #############
+
+ // ## 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 SquareMatrix<Type>, and print output
+template<class Type>
+void test_global_opers(Type)
+{
+ SquareMatrix<Type> A(3, Zero);
+ assignMatrix
+ (
+ A,
+ {
+ Type(1), Type(-2.2), Type(-3.4),
+ Type(-0.35), Type(10.32), Type(5),
+ Type(-0.35), Type(10.32), Type(-5)
+ }
+ );
+ const SquareMatrix<Type> cpA(A);
+
+ Info<< "# Operand: " << nl
+ << " SquareMatrix = " << A << endl;
+
+
+ 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(SquareMatrix<Type>(A.sizes(), Type(1)))
+ );
+ cmp
+ (
+ " Scalar addition = ",
+ flt(Type(1) + Type(0)*A),
+ flt(SquareMatrix<Type>(A.sizes(), Type(1)))
+ );
+ cmp
+ (
+ " Scalar subtraction = ",
+ flt(Type(0)*A - Type(1)),
+ flt(SquareMatrix<Type>(A.sizes(), Type(-1)))
+ );
+ cmp
+ (
+ " Scalar subtraction = ",
+ flt(Type(1) - Type(0)*A),
+ flt(SquareMatrix<Type>(A.sizes(), Type(1)))
+ );
+ cmp
+ (
+ " Scalar division = ",
+ flt((Type(1) + Type(0)*A)/Type(2.5)),
+ flt(SquareMatrix<Type>(A.sizes(), Type(0.4)))
+ );
+
+ SquareMatrix<Type> A1(3, Zero);
+ assignMatrix
+ (
+ A1,
+ {
+ Type(2.9600), Type(-59.9920), Type(2.6000),
+ Type(-5.71199), Type(158.87239), Type(27.789997),
+ Type(-2.211999), Type(55.672393), Type(77.789993)
+ }
+ );
+ cmp
+ (
+ " Matrix-Matrix multiplication using ikj-algorithm = ",
+ flt(A*A),
+ flt(A1),
+ 1e-4
+ );
+ SquareMatrix<Type> A2(3, Zero);
+ assignMatrix
+ (
+ A2,
+ {
+ Type(1.245), Type(-9.424), Type(-3.4),
+ Type(-9.424), Type(217.8448), Type(7.47999),
+ Type(-3.4), Type(7.47999), Type(61.56)
+ }
+ );
+ cmp
+ (
+ " Implicit A.T()*B = ",
+ flt(A&A),
+ flt(A2),
+ 1e-4
+ );
+ SquareMatrix<Type> A3(3, Zero);
+ assignMatrix
+ (
+ A3,
+ {
+ Type(17.4), Type(-40.054), Type(-6.054),
+ Type(-40.054), Type(131.6249), Type(81.6249),
+ Type(-6.054), Type(81.6249), Type(131.6249)
+ }
+ );
+ cmp
+ (
+ " Implicit A*B.T() = ",
+ flt(A^A),
+ flt(A3),
+ 1e-4
+ );
+
+ const List<Type> lst({Type(1), Type(2), Type(-3)});
+ const List<Type> out1(A*lst);
+ const List<Type> out2(lst*A);
+ const List<Type> rslt1({Type(6.7999), Type(5.2900), Type(35.29)});
+ const List<Type> rslt2({Type(1.35), Type(-12.52), Type(21.6)});
+ cmp
+ (
+ " Matrix-vector multiplication (A * x), where x is a column vector = ",
+ out1,
+ rslt1,
+ 1e-2
+ );
+ cmp
+ (
+ " Vector-matrix multiplication (x * A), where x is a row vector = ",
+ out2,
+ rslt2,
+ 1e-4
+ );
+}
+
+
+// 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
+ ({
+ "SquareMatrix<floatScalar>",
+ "SquareMatrix<doubleScalar>",
+ "SquareMatrix<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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/matrices/SymmetricSquareMatrix/Make/files b/applications/test/matrices/SymmetricSquareMatrix/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..da88933c249e3bfafd2a25c4842babc8f5c5c1f5
--- /dev/null
+++ b/applications/test/matrices/SymmetricSquareMatrix/Make/files
@@ -0,0 +1,3 @@
+Test-SymmetricSquareMatrix.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-SymmetricSquareMatrix
diff --git a/applications/test/matrices/SymmetricSquareMatrix/Make/options b/applications/test/matrices/SymmetricSquareMatrix/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..b23590707c8d302fe75ce411d5bf921c029a7ff3
--- /dev/null
+++ b/applications/test/matrices/SymmetricSquareMatrix/Make/options
@@ -0,0 +1,2 @@
+EXE_INC = -I../../TestTools
+/* EXE_LIBS = -lfiniteVolume */
diff --git a/applications/test/matrices/SymmetricSquareMatrix/Test-SymmetricSquareMatrix.C b/applications/test/matrices/SymmetricSquareMatrix/Test-SymmetricSquareMatrix.C
new file mode 100644
index 0000000000000000000000000000000000000000..c3d47e59b08e4794f2d3ec3f7998d61b62102ff2
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
index 26d92fba2ffcfc05071aafc29e56f0995d7a1e86..1d5146b4b7170daead7531ccc0590112a483706e 100644
--- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
+++ b/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"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type>
-inline Foam::DiagonalMatrix<Type>::DiagonalMatrix()
-:
- List<Type>()
-{}
-
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
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>
+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>
-Foam::DiagonalMatrix<Type> Foam::inv(const DiagonalMatrix<Type>& mat)
+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
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H
index b8e29d8270692b4fd1fd1333e02fd27a6e8035af..e08800bd77f8101b4bc59fd02af7cf4e87182ac6 100644
--- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H
+++ b/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
+
+ //- Default construct
+ DiagonalMatrix() = default;
+
+ //- Copy construct
+ DiagonalMatrix(const DiagonalMatrix&) = default;
+
+ //- Copy assignment
+ DiagonalMatrix& operator=(const DiagonalMatrix&) = default;
+
- //- Construct null.
- DiagonalMatrix<Type>();
+ // Constructors
//- Construct empty from size
explicit DiagonalMatrix<Type>(const label n);
- //- Construct from size
- //- initializing all elements to the given value
+ //- Construct from size and initialise all elems to zero
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,17 +94,18 @@ public:
// Member Functions
- //- Invert the diagonal matrix and return itself
- DiagonalMatrix<Type>& invert();
-
-};
-
+ //- Return the matrix inverse into itself if no elem is equal to zero
+ void invert();
-// Global functions
+ //- Return a sort permutation labelList according to
+ //- a given comparison on the diagonal entries
+ template<class CompOp>
+ List<label> sortPermutation(CompOp& compare) const;
-//- Return the diagonal Matrix inverse
-template<class Type>
-DiagonalMatrix<Type> inv(const DiagonalMatrix<Type>& mat);
+ //- Column-reorder this Matrix according to
+ //- a given permutation labelList
+ void applyPermutation(const List<label>& p);
+};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index b5ea2e4519d84f1b3cc700ab9c836ffd28cea19a..03beb2c5026b595c15d52f05f4423d728f701ec5 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/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 operator+(const Matrix<Form, Type>& mat, const Type& s)
+{
+ return s + mat;
+}
+
+
+//- Scalar subtraction 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());
@@ -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)
-{
- 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 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
-Foam::operator*
+typename typeOfInnerProduct<Type, Form1, Form2>::type
+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
-Foam::operator&
+typename typeOfInnerProduct<Type, Form1, Form2>::type
+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
-Foam::operator^
+typename typeOfInnerProduct<Type, Form1, Form2>::type
+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"
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index bb66c5e5b7134c8959b1a1d6436ff9e939f1d4a9..673258a7c737baa135dc42f04cc02a74c008410f 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/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
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index 435b22b7c5c07e42ce264b963c7788c355427d2b..ab2909a6738cfd8e22602f9df24b162ca76cda54 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/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
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C
index 7b6ab90153ffb581eec6d94c7fbc15688ad409da..54b88b39ae7ce7ac36b7a178e08846992ac15433 100644
--- a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C
+++ b/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())
{
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
index 5b09ea7568b0e153706f494c4b40122a2a4b2843..b5394f017aa69929b9562bbdedeae619b81cbd79 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
+++ b/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
+
+ //- Default construct
+ RectangularMatrix() = default;
+
+ //- Copy construct
+ RectangularMatrix(const RectangularMatrix&) = default;
+
+ //- Copy assignment
+ RectangularMatrix& operator=(const RectangularMatrix&) = default;
+
- //- Construct null
- inline RectangularMatrix();
+ // 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
- //- initializing all elements to zero
+ //- Construct given number of rows/columns by using a label pair
+ //- and initializing all elements to zero
inline RectangularMatrix(const labelPair& dims, const zero);
- //- Construct given number of rows/columns
- //- initializing all elements to the given value
+ //- Construct given number of rows/columns by using a label pair
+ //- 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
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
index 7cef3593e77c223e6e24c17022a14b6797ab9374..9b98e72caaa5fcee6b1dcfee4595abf711a89a36 100644
--- a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
+++ b/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
-
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
index 4bc1f0f78f6be55b324467e0ea42e17bd941f885..1506cb997236626c4341d5ccaeb6e27bec5ffa95 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
+++ b/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,10 +30,80 @@ License
#include "RectangularMatrix.H"
#include "labelList.H"
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class Type>
+template<class CompOp>
+Foam::List<Foam::label> Foam::SquareMatrix<Type>::sortPermutation
+(
+ CompOp& compare
+) const
+{
+ 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)); }
+ );
+
+ return p;
+}
+
+
+template<class Type>
+void Foam::SquareMatrix<Type>::applyPermutation(const List<label>& p)
+{
+ #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
+
+ SquareMatrix<Type> reordered(this->sizes());
+
+ label j = 0;
+ for (const label i : p)
+ {
+ reordered.subColumn(j) = this->subColumn(i);
+ ++j;
+ }
+
+ this->transfer(reordered);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+
+template<class Type>
+template<class AnyType>
+void Foam::SquareMatrix<Type>::operator=(const Identity<AnyType>)
+{
+ Matrix<SquareMatrix<Type>, Type>::operator=(Zero);
+
+ for (label i = 0; i < this->n(); ++i)
+ {
+ this->operator()(i, i) = pTraits<Type>::one;
+ }
+}
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
+
+//- Return the determinant of the LU decomposed SquareMatrix
template<class Type>
-Foam::scalar Foam::detDecomposed
+scalar detDecomposed
(
const SquareMatrix<Type>& matrix,
const label sign
@@ -50,8 +120,9 @@ Foam::scalar Foam::detDecomposed
}
+//- Return the determinant of SquareMatrix
template<class Type>
-Foam::scalar Foam::det(const SquareMatrix<Type>& matrix)
+scalar det(const SquareMatrix<Type>& matrix)
{
SquareMatrix<Type> matrixTmp = matrix;
@@ -63,8 +134,9 @@ Foam::scalar Foam::det(const SquareMatrix<Type>& matrix)
}
+//- Return the SquareMatrix det and the LU decomposition in the original matrix
template<class Type>
-Foam::scalar Foam::det(SquareMatrix<Type>& matrix)
+scalar det(SquareMatrix<Type>& matrix)
{
labelList pivotIndices(matrix.m());
label sign;
@@ -74,19 +146,52 @@ Foam::scalar Foam::det(SquareMatrix<Type>& matrix)
}
-// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
-
+//- Return Matrix column-reordered according to
+//- a given permutation labelList
template<class Type>
-template<class AnyType>
-void Foam::SquareMatrix<Type>::operator=(const Identity<AnyType>)
+SquareMatrix<Type> applyPermutation
+(
+ const SquareMatrix<Type>& mat,
+ const List<label>& p
+)
{
- Matrix<SquareMatrix<Type>, Type>::operator=(Zero);
+ #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
- for (label i=0; i < this->n(); ++i)
+ SquareMatrix<Type> reordered(mat.sizes());
+
+ label j = 0;
+ for (const label i : p)
{
- this->operator()(i, i) = pTraits<Type>::one;
+ 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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
index e404910e069a8e4fad048a10d2056e064fcd4ed8..a79f2d9f4230d6a8335adfb5390c6a2610be9ddb 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
+++ b/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
- dimension is known and used for subscript bounds checking, etc.
+ A templated (N x N) square matrix of objects of \<Type\>,
+ 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
+
+ //- Default construct
+ SquareMatrix() = default;
+
+ //- Copy construct
+ SquareMatrix(const SquareMatrix&) = default;
- //- Construct null
- inline SquareMatrix();
+ //- 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)
- //- initializing all elements to zero
+ //- Construct given number of rows/columns
+ //- 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)
- //- initializing all elements to the given value
+ //- Construct given number of rows/columns
+ //- 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
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
index 75e893e068f8c9d91e32d83aaec652364b967455..61933881e1139aced6e9e74a8caa7403e29f917c 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
+++ b/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)
:
@@ -83,7 +76,7 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
:
Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
{
- for (label i=0; i < n; ++i)
+ for (label i = 0; i < n; ++i)
{
this->operator()(i, i) = pTraits<Type>::one;
}
@@ -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
(
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
index de6291e361b6e09160ff5ca3c7305e66238c3ef3..b50a96e7b1b3f009c60147cc15d4d8a99b6a8e30 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
+++ b/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>)
}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
index 49af9fae42f1c897b9aed10cec78ab2919cd4d97..c04ec0c8559f1a6677847b5d5006b5609afe3ed5 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
+++ b/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
- n x n matrix dimension is known and used for subscript bounds checking, etc.
+ A templated (N x N) square matrix of objects of \<Type\>,
+ 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
+
+ //- Default construct
+ SymmetricSquareMatrix() = default;
+
+ //- Copy construct
+ SymmetricSquareMatrix(const SymmetricSquareMatrix&) = default;
- //- Construct null
- inline SymmetricSquareMatrix();
+ //- 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>);
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
index 79c6b2a05e766529fce29885ff0c56d0ab48279f..31eff6076595dfecab80516349f2e9196614daa6 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
+++ b/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)
: