Commit 55e7da67 authored by Kutalmis Bercin's avatar Kutalmis Bercin Committed by Andrew Heather
Browse files

ENH: improve analytical eigendecompositions

  - `tensor` and `tensor2D` returns complex eigenvalues/vectors
  - `symmTensor` and `symmTensor2D` returns real eigenvalues/vectors
  - adds new test routines for eigendecompositions
  - improves numerical stability by:
    - using new robust algorithms,
    - reordering the conditional branches in root-type selection
parent 6a53794e
......@@ -512,6 +512,78 @@ void test_global_opers(Type)
}
// Return false if given eigenvalues fail to satisy eigenvalue relations
// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307)
void test_eigenvalues(const symmTensor& sT, const vector& EVals)
{
{
const scalar determinant = det(sT);
const scalar EValsProd = EVals.x()*EVals.y()*EVals.z();
cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant);
}
{
const scalar trace = tr(sT);
scalar EValsSum = 0.0;
for (const auto& val : EVals)
{
EValsSum += val;
}
cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace);
}
}
// Return false if a given eigenvalue-eigenvector pair
// fails to satisfy the characteristic equation
void test_characteristic_equation
(
const symmTensor& sT,
const vector& EVals,
const tensor& EVecs
)
{
Info<< "# Characteristic equation:" << nl;
for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
{
Info<< "EVal = " << EVals[dir] << nl
<< "EVec = " << EVecs.row(dir) << endl;
const vector leftSide(sT & EVecs.row(dir));
const vector rightSide(EVals[dir]*EVecs.row(dir));
const vector X(leftSide - rightSide);
for (const auto x : X)
{
cmp(" (sT & EVec - EVal*EVec) = 0:", x, 0.0, 1e-8, 1e-6);
}
}
}
// Return false if the eigen functions fail to satisfy relations
void test_eigen_funcs(const symmTensor& sT)
{
Info<< "# Operand:" << nl
<< " symmTensor = " << sT << nl;
Info<< "# Return eigenvalues of a given symmTensor:" << nl;
const vector EVals(eigenValues(sT));
Info<< EVals << endl;
test_eigenvalues(sT, EVals);
Info<< "# Return eigenvectors of a given symmTensor corresponding to"
<< " given eigenvalues:" << nl;
const tensor EVecs0(eigenVectors(sT, EVals));
Info<< EVecs0 << endl;
test_characteristic_equation(sT, EVals, EVecs0);
Info<< "# Return eigenvectors of a given symmTensor by computing"
<< " the eigenvalues of the symmTensor in the background:" << nl;
const tensor EVecs1(eigenVectors(sT));
Info<< EVecs1 << endl;
}
// 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
......@@ -557,6 +629,74 @@ int main()
run_tests(types, typeID);
Info<< nl << " ## Test SymmTensor<scalar> eigen functions: ##" << nl;
const label numberOfTests = 10000;
Random rndGen(1234);
for (label i = 0; i < numberOfTests; ++i)
{
const symmTensor T(makeRandomContainer(rndGen));
test_eigen_funcs(T);
}
{
Info<< nl << " ## Test eigen functions by a zero tensor: ##"<< nl;
const symmTensor zeroT(Zero);
test_eigen_funcs(zeroT);
}
{
Info<< nl
<< " ## Test eigen functions by a tensor consisting of"
<< " 2 repeated eigenvalues: ##"
<< nl;
const symmTensor T
(
1.0, 0.0, Foam::sqrt(2.0),
2.0, 0.0,
0.0
);
test_eigen_funcs(T);
}
{
Info<< nl
<< " ## Test eigen functions by a tensor consisting of"
<< " 3 repeated eigenvalues: ##"
<< nl;
const symmTensor T
(
0.023215, -5.0739e-09, -7.0012e-09,
0.023215, -8.148e-10,
0.023215
);
test_eigen_funcs(T);
}
{
Info<< nl
<< " ## Test eigen functions by a tensor consisting of"
<< " SMALL off-diagonal elements: ##"
<< nl;
for (label i = 0; i < numberOfTests; ++i)
{
symmTensor T(makeRandomContainer(rndGen));
T.xy() = SMALL*rndGen.GaussNormal<scalar>();
T.xz() = SMALL*rndGen.GaussNormal<scalar>();
T.yz() = SMALL*rndGen.GaussNormal<scalar>();
test_eigen_funcs(T);
}
}
{
Info<< nl << " ## Test eigen functions by a stiff tensor: ##" << nl;
const symmTensor stiff
(
pow(10.0, 10), pow(10.0, 8), pow(10.0, -8),
pow(10.0, -8), pow(10.0, 8),
pow(10.0, 7)
);
test_eigen_funcs(stiff);
}
if (nFail_)
{
Info<< nl << " #### "
......@@ -568,8 +708,6 @@ int main()
Info<< nl << " #### Passed all " << nTest_ <<" tests ####\n" << endl;
return 0;
return 0;
}
......
......@@ -476,6 +476,80 @@ void test_global_opers(Type)
}
// Return false if given eigenvalues fail to satisy eigenvalue relations
// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307)
void test_eigenvalues(const symmTensor2D& T, const vector2D& EVals)
{
{
Info<< "# Product of eigenvalues = det(T):" << nl;
const scalar determinant = det(T);
const scalar EValsProd = EVals.x()*EVals.y();
cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant);
}
{
const scalar trace = tr(T);
scalar EValsSum = 0.0;
for (const auto& val : EVals)
{
EValsSum += val;
}
cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace);
}
}
// Return false if a given eigenvalue-eigenvector pair
// fails to satisfy the characteristic equation
void test_characteristic_equation
(
const symmTensor2D& T,
const vector2D& EVals,
const tensor2D& EVecs
)
{
Info<< "# Characteristic equation:" << nl;
for (direction dir = 0; dir < pTraits<vector2D>::nComponents; ++dir)
{
Info<< "EVal = " << EVals[dir] << nl
<< "EVec = " << EVecs.row(dir) << nl;
const vector2D leftSide(T & EVecs.row(dir));
const vector2D rightSide(EVals[dir]*EVecs.row(dir));
const vector2D X(leftSide - rightSide);
for (const auto x : X)
{
cmp(" (sT & EVec - EVal*EVec) = 0:", x, 0.0, 1e-8, 1e-6);
}
}
}
// Return false if the eigen functions fail to satisfy relations
void test_eigen_funcs(const symmTensor2D& T)
{
Info<< "# Operand:" << nl
<< " symmTensor2D = " << T << nl;
Info<< "# Return eigenvalues of a given symmTensor2D:" << nl;
const vector2D EVals(eigenValues(T));
Info<< EVals << endl;
test_eigenvalues(T, EVals);
Info<< "# Return eigenvectors of a given symmTensor2D corresponding to"
<< " given eigenvalues:" << nl;
const tensor2D EVecs0(eigenVectors(T, EVals));
Info<< EVecs0 << endl;
test_characteristic_equation(T, EVals, EVecs0);
Info<< "# Return eigenvectors of a given symmTensor2D by computing"
<< " the eigenvalues of the symmTensor2D in the background:" << nl;
const tensor2D EVecs1(eigenVectors(T));
Info<< EVecs1 << endl;
}
// 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
......@@ -521,6 +595,52 @@ int main(int argc, char *argv[])
run_tests(types, typeID);
Info<< nl << " ## Test symmTensor2D eigen functions: ##" << nl;
const label numberOfTests = 10000;
Random rndGen(1234);
for (label i = 0; i < numberOfTests; ++i)
{
const symmTensor2D T(makeRandomContainer(rndGen));
test_eigen_funcs(T);
}
Info<< nl << " ## Test symmTensor2D eigen functions"
<< " with T.xy = VSMALL: ##" << nl;
for (label i = 0; i < numberOfTests; ++i)
{
symmTensor2D T(makeRandomContainer(rndGen));
T.xy() = VSMALL;
test_eigen_funcs(T);
}
Info<< nl << " ## Test symmTensor2D eigen functions"
<< " with T.xx = T.yy: ##" << nl;
for (label i = 0; i < numberOfTests; ++i)
{
symmTensor2D T(makeRandomContainer(rndGen));
T.xx() = T.yy();
test_eigen_funcs(T);
}
{
Info<< nl
<< " ## Test eigen functions by a zero symmTensor2D: ##"<< nl;
const symmTensor2D zeroT(Zero);
test_eigen_funcs(zeroT);
}
{
Info<< nl << " ## Test eigen functions by a stiff tensor2D: ##"<< nl;
const symmTensor2D stiff
(
pow(10.0, 10), pow(10.0, -8),
pow(10.0, 9)
);
test_eigen_funcs(stiff);
}
if (nFail_)
{
Info<< nl << " #### "
......
......@@ -902,6 +902,104 @@ void test_global_opers(Type)
}
// Return false if given eigenvalues fail to satisy eigenvalue relations
// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307)
void test_eigenvalues
(
const tensor& T,
const Vector<complex>& EVals,
const bool prod = true
)
{
if (prod)
{
const scalar determinant = det(T);
// In case of complex EVals, the production is effectively scalar
// due to the (complex*complex conjugate) results in zero imag part
const scalar EValsProd = ((EVals.x()*EVals.y()*EVals.z()).real());
cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant);
}
{
const scalar trace = tr(T);
scalar EValsSum = 0.0;
// In case of complex EVals, the summation is effectively scalar
// due to the (complex+complex conjugate) results in zero imag part
for (const auto& val : EVals)
{
EValsSum += val.real();
}
cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace);
}
}
// Return false if a given eigenvalue-eigenvector pair
// fails to satisfy the characteristic equation
void test_characteristic_equation
(
const tensor& T,
const Vector<complex>& EVals,
const Tensor<complex>& EVecs
)
{
Info<< "# Characteristic equation:" << nl;
Tensor<complex> Tc(Zero);
forAll(T, i)
{
Tc[i] = complex(T[i], 0);
}
for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
{
const Vector<complex> leftSide(Tc & EVecs.row(dir));
const Vector<complex> rightSide(EVals[dir]*EVecs.row(dir));
const Vector<complex> X(leftSide - rightSide);
for (const auto x : X)
{
cmp(" (sT & EVec - EVal*EVec) = 0:", mag(x), 0.0, 1e-8, 1e-6);
}
}
}
// Return false if the eigen functions fail to satisfy relations
void test_eigen_funcs(const tensor& T, const bool prod = true)
{
Info<< "# Operand:" << nl
<< " tensor = " << T << nl;
Info<< "# Return eigenvalues of a given tensor:" << nl;
const Vector<complex> EVals(eigenValues(T));
Info<< EVals << endl;
test_eigenvalues(T, EVals, prod);
Info<< "# Return an eigenvector of a given tensor in a given direction"
<< " corresponding to a given eigenvalue:" << nl;
const Vector<complex> standardBasis1(Zero, pTraits<complex>::one, Zero);
const Vector<complex> standardBasis2(Zero, Zero, pTraits<complex>::one);
const Vector<complex> EVec
(
eigenVector(T, EVals.x(), standardBasis1, standardBasis2)
);
Info<< EVec << endl;
Info<< "# Return eigenvectors of a given tensor corresponding to"
<< " given eigenvalues:" << nl;
const Tensor<complex> EVecs0(eigenVectors(T, EVals));
Info<< EVecs0 << endl;
test_characteristic_equation(T, EVals, EVecs0);
Info<< "# Return eigenvectors of a given tensor by computing"
<< " the eigenvalues of the tensor in the background:" << nl;
const Tensor<complex> EVecs1(eigenVectors(T));
Info<< EVecs1 << endl;
}
// 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
......@@ -947,6 +1045,53 @@ int main()
run_tests(types, typeID);
Info<< nl << " ## Test tensor eigen functions: ##" << nl;
const label numberOfTests = 10000;
Random rndGen(1234);
for (label i = 0; i < numberOfTests; ++i)
{
const tensor T(makeRandomContainer(rndGen));
test_eigen_funcs(T);
}
{
Info<< nl << " ## Test eigen functions by a zero tensor: ##"<< nl;
const tensor zeroT(Zero);
test_eigen_funcs(zeroT);
}
{
Info<< nl
<< " ## Test eigen functions by a skew-symmetric tensor"
<< " consisting of no-real eigenvalues: ##"
<< nl;
const tensor T
(
0, 1, 1,
-1, 0, 1,
-1, -1, 0
);
test_eigen_funcs(T);
}
{
Info<< nl
<< " ## Test eigen functions by a stiff tensor: ##"
<< nl;
const tensor stiff
(
pow(10.0, 10), pow(10.0, 8), pow(10.0, 7),
pow(10.0, -8), pow(10.0, 9), pow(10.0, -8),
pow(10.0, 10), pow(10.0, 8), pow(10.0, 7)
);
// Although eigendecomposition is successful for the stiff tensors,
// cross-check between prod(eigenvalues) ?= det(stiff) is inherently
// problematic; therefore, eigenvalues of the stiff tensors are
// cross-checked by only sum(eigenvalues) ?= trace(stiff)
const bool testProd = false;
test_eigen_funcs(stiff, testProd);
}
if (nFail_)
{
Info<< nl << " #### "
......
......@@ -716,6 +716,94 @@ void test_global_opers(Type)
}
// Return false if given eigenvalues fail to satisy eigenvalue relations
// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307)
void test_eigenvalues(const tensor2D& T, const Vector2D<complex>& EVals)
{
{
const scalar determinant = det(T);
// In case of complex EVals, the production is effectively scalar
// due to the (complex*complex conjugate) results in zero imag part
const scalar EValsProd = ((EVals.x()*EVals.y()).real());
cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant);
}
{
const scalar trace = tr(T);
scalar EValsSum = 0.0;
// In case of complex EVals, the summation is effectively scalar
// due to the (complex+complex conjugate) results in zero imag part
for (const auto& val : EVals)
{
EValsSum += val.real();
}
cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace, 1e-8, 1e-8);
}
}
// Return false if a given eigenvalue-eigenvector pair
// fails to satisfy the characteristic equation
void test_characteristic_equation
(
const tensor2D& T,
const Vector2D<complex>& EVals,
const Tensor2D<complex>& EVecs
)
{
Info<< "# Characteristic equation:" << nl;
Tensor2D<complex> Tc(Zero);
forAll(T, i)
{
Tc[i] = complex(T[i], 0);
}
for (direction dir = 0; dir < pTraits<vector2D>::nComponents; ++dir)
{
const Vector2D<complex> leftSide(Tc & EVecs.row(dir));
const Vector2D<complex> rightSide(EVals[dir]*EVecs.row(dir));
const Vector2D<complex> X(leftSide - rightSide);
for (const auto x : X)
{
cmp(" (Tc & EVec - EVal*EVec) = 0:", mag(x), 0.0, 1e-8, 1e-6);
}
}
}
// Return false if the eigen functions fail to satisfy relations
void test_eigen_funcs(const tensor2D& T)
{
Info<< "# Operand:" << nl
<< " tensor2D = " << T << nl;
Info<< "# Return eigenvalues of a given tensor2D:" << nl;
const Vector2D<complex> EVals(eigenValues(T));
Info<< EVals << endl;
test_eigenvalues(T, EVals);
Info<< "# Return an eigenvector of a given tensor2D in a given direction"
<< " corresponding to a given eigenvalue:" << nl;
const Vector2D<complex> standardBasis(pTraits<complex>::one, Zero);
const Vector2D<complex> EVec(eigenVector(T, EVals.x(), standardBasis));
Info<< EVec << endl;
Info<< "# Return eigenvectors of a given tensor2D corresponding to"
<< " given eigenvalues:" << nl;
const Tensor2D<complex> EVecs0(eigenVectors(T, EVals));
Info<< EVecs0 << endl;
test_characteristic_equation(T, EVals, EVecs0);
Info<< "# Return eigenvectors of a given tensor2D by computing"
<< " the eigenvalues of the tensor2D in the background:" << nl;
const Tensor2D<complex> EVecs1(eigenVectors(T));
Info<< EVecs1 << endl;
}
// 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
......@@ -762,6 +850,58 @@ int main(int argc, char *argv[])
run_tests(types, typeID);
Info<< nl << " ## Test tensor2D eigen functions: ##" << nl;
const label numberOfTests = 10000;
Random rndGen(1234);
for (label i = 0; i < numberOfTests; ++i)
{
const tensor2D T(makeRandomContainer(rndGen));
test_eigen_funcs(T);
}
{
Info<< nl << " ## Test eigen functions by a zero tensor2D: ##"<< nl;
const tensor2D zeroT(Zero);
test_eigen_funcs(zeroT);
}
{
Info<< nl
<< " ## Test eigen functions by a tensor2D consisting of"
<< " repeated eigenvalues: ##"
<< nl;
const tensor2D T
(
-1.0, 2.0,
0.0, -1.0
);
test_eigen_funcs(T);
}
{
Info<< nl
<< " ## Test eigen functions by a skew-symmetric tensor2D"
<< " consisting of no-real eigenvalues: ##"
<< nl;
const tensor2D T
(
0.0, 1.0,
-1.0, 0.0
);