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