Commit 64614cfc authored by Kutalmis Bercin's avatar Kutalmis Bercin
Browse files

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
parent f9e53921
......@@ -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<SMatrix>({mRows, mRows}, rndGen));
Info<< "# SquareMatrix<scalar>.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<scalar>.symmetric():" << nl
<< A.symmetric() << nl;
}
{
const label mRows = 100;
SCMatrix A(mRows);
for (auto& val : A)
{
val.Re() = rndGen.GaussNormal<scalar>();
val.Im() = rndGen.GaussNormal<scalar>();
}
Info<< "# SquareMatrix<complex>.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<complex>.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);
......
......@@ -85,6 +85,9 @@ public:
template<class AnyType>
inline SquareMatrix(const label n, const Identity<AnyType>);
template<class AnyType>
inline SquareMatrix(const labelPair& dims, const Identity<AnyType>);
//- 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
......
......@@ -90,6 +90,25 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
}
template<class Type>
template<class AnyType>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
const labelPair& dims,
const Identity<AnyType>
)
:
Matrix<SquareMatrix<Type>, Type>(dims, Zero)
{
CHECK_MATRIX_IS_SQUARE(dims.first(), dims.second());
for (label i = 0; i < dims.first(); ++i)
{
this->operator()(i, i) = pTraits<Type>::one;
}
}
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
......@@ -206,6 +225,19 @@ inline void Foam::SquareMatrix<Type>::resize(const label m)
}
template<class Type>
inline void Foam::SquareMatrix<Type>::resize(const label m, const label n)
{
if (m != n)
{
FatalErrorInFunction<< "Disallowed use of resize() for SquareMatrix"
<< abort(FatalError);
}
Matrix<SquareMatrix<Type>, Type>::resize(m, m);
}
template<class Type>
inline void Foam::SquareMatrix<Type>::setSize(const label m)
{
......@@ -220,6 +252,49 @@ inline void Foam::SquareMatrix<Type>::shallowResize(const label m)
}
template<class Type>
inline bool Foam::SquareMatrix<Type>::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<class Type>
inline bool Foam::SquareMatrix<Type>::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<class Type>
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment