From 64614cfce88fe668ee728ef1d6de273330cfc5b5 Mon Sep 17 00:00:00 2001 From: Kutalmis Bercin Date: Wed, 13 Nov 2019 15:47:36 +0000 Subject: [PATCH] ENH: add new funcs into `SquareMatrix` - query func `symmetric()` - query func `tridiagonal()` - `resize()` - `labelpair` identity constructor STYLE: add `#if(0 | RUNALL)` to improve test control in Test-Matrix --- applications/test/Matrix/Test-Matrix.C | 79 ++++++++++++++++--- .../matrices/SquareMatrix/SquareMatrix.H | 14 ++++ .../matrices/SquareMatrix/SquareMatrixI.H | 75 ++++++++++++++++++ 3 files changed, 158 insertions(+), 10 deletions(-) diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C index c372cf9b9a..1cbc602545 100644 --- a/applications/test/Matrix/Test-Matrix.C +++ b/applications/test/Matrix/Test-Matrix.C @@ -37,6 +37,7 @@ License using namespace Foam; using namespace Foam::MatrixTools; +#define RUNALL true #define isEqual MatrixTools::equal void horizontalLine() @@ -113,7 +114,7 @@ int main(int argc, char *argv[]) Random rndGen(1234); - #if 1 + #if (0 | RUNALL) { horizontalLine(); @@ -218,7 +219,7 @@ int main(int argc, char *argv[]) #endif - #if 1 + #if (0 | RUNALL) { horizontalLine(); @@ -254,7 +255,7 @@ int main(int argc, char *argv[]) #endif - #if 1 + #if (0 | RUNALL) { horizontalLine(); @@ -299,7 +300,7 @@ int main(int argc, char *argv[]) #endif - #if 1 + #if (0 | RUNALL) { horizontalLine(); @@ -410,7 +411,7 @@ int main(int argc, char *argv[]) #endif - #if 1 + #if (0 | RUNALL) { horizontalLine(); @@ -494,7 +495,7 @@ int main(int argc, char *argv[]) #endif - #if 1 + #if (0 | RUNALL) { horizontalLine(); @@ -683,7 +684,7 @@ int main(int argc, char *argv[]) #endif - #if 1 + #if (0 | RUNALL) { horizontalLine(); @@ -804,7 +805,7 @@ int main(int argc, char *argv[]) << nl; } - #if 1 + #if (0 | RUNALL) { Info<< nl << "# Implicit inner/outer products:" << nl; @@ -860,7 +861,65 @@ int main(int argc, char *argv[]) #endif - #if 1 + #if (0 | RUNALL) + { + horizontalLine(); + + Info<< "## Query to determine if a square matrix is symmetric:" << nl; + + { + const label mRows = 100; + SMatrix A(makeRandomMatrix({mRows, mRows}, rndGen)); + + Info<< "# SquareMatrix.symmetric():" << nl + << A.symmetric() << nl; + + // Symmetrise + for (label n = 0; n < A.n() - 1; ++n) + { + for (label m = A.m() - 1; n < m; --m) + { + A(n, m) = A(m, n); + } + } + + Info<< "# SquareMatrix.symmetric():" << nl + << A.symmetric() << nl; + } + + { + const label mRows = 100; + + SCMatrix A(mRows); + + for (auto& val : A) + { + val.Re() = rndGen.GaussNormal(); + val.Im() = rndGen.GaussNormal(); + } + + Info<< "# SquareMatrix.symmetric():" << nl + << A.symmetric() << nl; + + // Symmetrise + for (label n = 0; n < A.n() - 1; ++n) + { + for (label m = A.m() - 1; n < m; --m) + { + A(n, m) = A(m, n); + } + } + + Info<< "# SquareMatrix.symmetric():" << nl + << A.symmetric() << nl; + } + + horizontalLine(); + } + #endif + + + #if (0 | RUNALL) { horizontalLine(); @@ -898,7 +957,7 @@ int main(int argc, char *argv[]) #endif - #if 1 + #if (0 | RUNALL) { scalarSquareMatrix squareMatrix(3, Zero); diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H index b3ac5eb9f3..e404910e06 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H @@ -85,6 +85,9 @@ public: template inline SquareMatrix(const label n, const Identity); + template + inline SquareMatrix(const labelPair& dims, const Identity); + //- Construct given number of rows/columns (checked to be equal) // For constructor consistency in rectangular matrices inline explicit SquareMatrix(const labelPair& dims); @@ -129,12 +132,23 @@ public: //- Resize the matrix preserving the elements inline void resize(const label m); + //- Resize the matrix preserving the elements (compatibility) + inline void resize(const label m, const label n); + //- Resize the matrix preserving the elements inline void setSize(const label m); //- Resize the matrix without reallocating storage (unsafe) inline void shallowResize(const label m); + // Check + + //- Return true if the square matrix is effectively symmetric/Hermitian + inline bool symmetric() const; + + //- Return true if the square matrix is reduced tridiagonal + inline bool tridiagonal() const; + // Member Operators diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H index b41b5e97be..75e893e068 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H @@ -90,6 +90,25 @@ inline Foam::SquareMatrix::SquareMatrix } +template +template +inline Foam::SquareMatrix::SquareMatrix +( + const labelPair& dims, + const Identity +) +: + Matrix, Type>(dims, Zero) +{ + CHECK_MATRIX_IS_SQUARE(dims.first(), dims.second()); + + for (label i = 0; i < dims.first(); ++i) + { + this->operator()(i, i) = pTraits::one; + } +} + + template inline Foam::SquareMatrix::SquareMatrix ( @@ -206,6 +225,19 @@ inline void Foam::SquareMatrix::resize(const label m) } +template +inline void Foam::SquareMatrix::resize(const label m, const label n) +{ + if (m != n) + { + FatalErrorInFunction<< "Disallowed use of resize() for SquareMatrix" + << abort(FatalError); + } + + Matrix, Type>::resize(m, m); +} + + template inline void Foam::SquareMatrix::setSize(const label m) { @@ -220,6 +252,49 @@ inline void Foam::SquareMatrix::shallowResize(const label m) } +template +inline bool Foam::SquareMatrix::symmetric() const +{ + for (label n = 0; n < this->n() - 1; ++n) + { + for (label m = this->m() - 1; n < m; --m) + { + if (SMALL < mag((*this)(n, m) - (*this)(m, n))) + { + return false; + } + } + } + return true; +} + + +template +inline bool Foam::SquareMatrix::tridiagonal() const +{ + for (label i = 0; i < this->m(); ++i) + { + for (label j = 0; j < this->n(); ++j) + { + const Type& val = (*this)(i, j); + + if ((i == j) || (i - 1 == j) || (i + 1 == j)) + { + if (mag(val) < SMALL) + { + return false; + } + } + else if (SMALL < mag(val)) + { + return false; + } + } + } + return true; +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template