Skip to content

Commit

Permalink
review corrections + update stats for eigen solver tests + add vector…
Browse files Browse the repository at this point in the history 
…s tests
  • Loading branch information
maxence-wz committed Apr 16, 2024
1 parent 9be5fa9 commit 1522c83
Show file tree
Hide file tree
Showing 3 changed files with 80 additions and 68 deletions.
52 changes: 26 additions & 26 deletions docs/web/release-notes-5.0.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,10 +35,10 @@ const auto m3 = m1 * m2;

- New eigen solvers have been introduced.

### New eigen solvers{#sec:eigensolvers}
### New eigen solvers{#sec:tfel-5.0.0:tfel-math:eigensolvers}

New eigen solver based on Harari analytical solution have been introduced for symmetric tensors. The
computation of eigen values is done with the Harai algorithm [@harari_computation_2023] and the
computation of eigen values is done with Harari's algorithm [@harari_computation_2023] and the
computation of eigen vectors is done with the default eigen solver for symmetric tensors of `TFEL`. Such
computations are more efficient and more accurate than the default `TFEL` algorithm.

Expand Down Expand Up @@ -74,38 +74,38 @@ std::tie(vp,m)=s.computeEigenVectors<stensor::HARARIEIGENSOLVER>();

| Algorithm | Failure ratio | \(\Delta_{\infty}\) | Times (ns) | Time ratio |
| :-------------------------: | :-----------: | :-----------------: | :--------: | :--------: |
| `TFELEIGENSOLVER` | 0.000642 | 3.29e-05 | 250174564 | 1 |
| `GTESYMMETRICQREIGENSOLVER` | 0 | 1.10e-06 | 359854550 | 1.44 |
| `FSESJACOBIEIGENSOLVER` | 0 | 5.62e-07 | 473263841 | 1.89 |
| `FSESQLEIGENSOLVER` | 0.000397 | 1.69e-06 | 259080052 | 1.04 |
| `FSESCUPPENEIGENSOLVER` | 0.019469 | 3.49e-06 | 274547371 | 1.10 |
| `FSESHYBRIDEIGENSOLVER` | 0.076451 | 5.56e-03 | 126689850 | 0.51 |
| `FSESANALYTICALEIGENSOLVER` | 0.108877 | 1.58e-01 | 127236908 | 0.51 |
| `HARARIEIGENSOLVER` | 0.000018 | 4.04e-06 | 167626780 | 0.67 |
| `TFELEIGENSOLVER` | 0.000557 | 2.37e-05 | 129452335 | 1 |
| `GTESYMMETRICQREIGENSOLVER` | 0 | 9.57e-07 | 236544828 | 1.83 |
| `FSESJACOBIEIGENSOLVER` | 0 | 4.61e-07 | 241131631 | 1.86 |
| `FSESQLEIGENSOLVER` | 0.000173 | 1.67e-06 | 155435496 | 1.20 |
| `FSESCUPPENEIGENSOLVER` | 0.018223 | 2.87e-06 | 151727321 | 1.17 |
| `FSESHYBRIDEIGENSOLVER` | 0.068411 | 3.90e-03 | 80039266 | 0.62 |
| `FSESANALYTICALEIGENSOLVER` | 0.102626 | 6.21e-02 | 76865961 | 0.59 |
| `HARARIEIGENSOLVER` | 0.000018 | 2.46e-06 | 110028802 | 0.85 |
: Test on \(10^{6}\) random symmetric tensors in single precision (`float`). {#tbl:comp_eigensolvers_float}

| Algorithm | Failure ratio | \(\Delta_{\infty}\) | Times (ns) | Time ratio |
| :-------------------------: | :-----------: | :-----------------: | :--------: | :--------: |
| `TFELEIGENSOLVER` | 0.000632 | 7.75e-14 | 252663338 | 1 |
| `GTESYMMETRICQREIGENSOLVER` | 0 | 2.06e-15 | 525845499 | 2.08 |
| `FSESJACOBIEIGENSOLVER` | 0 | 1.05e-15 | 489507133 | 1.94 |
| `FSESQLEIGENSOLVER` | 0.000422 | 3.30e-15 | 367599140 | 1.45 |
| `FSESCUPPENEIGENSOLVER` | 0.020174 | 5.79e-15 | 374190684 | 1.48 |
| `FSESHYBRIDEIGENSOLVER` | 0.090065 | 3.53e-10 | 154911762 | 0.61 |
| `FSESANALYTICALEIGENSOLVER` | 0.110399 | 1.09e-09 | 157613994 | 0.62 |
| `HARARIEIGENSOLVER` | 0.000022 | 8.90e-15 | 175109610 | 0.69 |
| `TFELEIGENSOLVER` | 0.00058 | 6.94e-14 | 137752068 | 1 |
| `GTESYMMETRICQREIGENSOLVER` | 1e-06 | 2.30e-15 | 315593552 | 2.29 |
| `FSESJACOBIEIGENSOLVER` | 0 | 9.08e-16 | 256285090 | 1.86 |
| `FSESQLEIGENSOLVER` | 0.000202 | 3.04e-15 | 214537012 | 1.56 |
| `FSESCUPPENEIGENSOLVER` | 0.019251 | 5.58e-15 | 219113965 | 1.59 |
| `FSESHYBRIDEIGENSOLVER` | 0.081586 | 1.29e-10 | 81861668 | 0.59 |
| `FSESANALYTICALEIGENSOLVER` | 0.103935 | 4.11e-10 | 79701256 | 0.58 |
| `HARARIEIGENSOLVER` | 0.000037 | 2.27e-14 | 116977683 | 0.85 |
: Test on \(10^{6}\) random symmetric tensors in double precision (`double`). {#tbl:comp_eigensolvers_double}

| Algorithm | Failure ratio | \(\Delta_{\infty}\) | Times (ns) | Time ratio |
| :-------------------------: | :-----------: | :-----------------: | :--------: | :--------: |
| `TFELEIGENSOLVER` | 0.000575 | 2.06e-17 | 428333721 | 1 |
| `GTESYMMETRICQREIGENSOLVER` | 0 | 1.00e-18 | 814990447 | 1.90 |
| `FSESJACOBIEIGENSOLVER` | 0 | 5.11e-19 | 748476926 | 1.75 |
| `FSESQLEIGENSOLVER` | 0.00045 | 1.83e-18 | 548604588 | 1.28 |
| `FSESCUPPENEIGENSOLVER` | 0.009134 | 3.23e-18 | 734707748 | 1.71 |
| `FSESHYBRIDEIGENSOLVER` | 0.99959 | 4.01e-10 | 272701749 | 0.64 |
| `FSESANALYTICALEIGENSOLVER` | 0.999669 | 1.36e-11 | 315243286 | 0.74 |
| `HARARIEIGENSOLVER` | 0.000044 | 9.49e-18 | 459905551 | 1.07 |
| `TFELEIGENSOLVER` | 0.000578 | 1.76e-17 | 304165023 | 1 |
| `GTESYMMETRICQREIGENSOLVER` | 0 | 1.01e-18 | 558605772 | 1.84 |
| `FSESJACOBIEIGENSOLVER` | 0 | 5.11e-19 | 408584703 | 1.34 |
| `FSESQLEIGENSOLVER` | 0.00045 | 1.83e-18 | 311028180 | 1.02 |
| `FSESCUPPENEIGENSOLVER` | 0.009134 | 3.23e-18 | 485590150 | 1.60 |
| `FSESHYBRIDEIGENSOLVER` | 0.99959 | 4.01e-10 | 187308886 | 0.62 |
| `FSESANALYTICALEIGENSOLVER` | 0.999669 | 1.36e-11 | 211710377 | 0.70 |
| `HARARIEIGENSOLVER` | 0.000046 | 4.62e-18 | 338589049 | 1.11 |
: Test on \(10^{6}\) random symmetric tensors in extended precision (`long double`). {#tbl:comp_eigensolvers_long_double}


Expand Down
63 changes: 21 additions & 42 deletions include/TFEL/Math/Stensor/Internals/HarariEigenSolver.ixx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,28 +33,17 @@ namespace tfel::math::internals {
const real E,
const real F) {
constexpr auto one = real{1};
constexpr const auto one_half = one / 2;
constexpr const auto one_third = one / 3;

tmatrix<3u, 3u, real> eye;
eye(0, 0) = one;
eye(1, 1) = one;
eye(2, 2) = one;
constexpr auto eye = tmatrix<3u, 3u, real>::Id();

// compute the trace of A
const real I1 = (A + B + C);
const auto I1 = (A + B + C);
const auto tr = one_third * I1;

// compute deviatoric part of M
tmatrix<3u, 3u, real> S;
S(0, 0) = (A - I1 / 3);
S(1, 1) = (B - I1 / 3);
S(2, 2) = (C - I1 / 3);
S(0, 1) = D;
S(0, 2) = E;
S(1, 2) = F;
S(1, 0) = D;
S(2, 0) = E;
S(2, 1) = F;
const auto S = tmatrix<3u, 3u, real>{A - tr, D, E, //
D, B - tr, F, //
E, F, C - tr};

// compute second invariant : J2
const real J2 =
Expand All @@ -65,23 +54,15 @@ namespace tfel::math::internals {
const real s = std::sqrt(J2 / 3);

if (s < std::numeric_limits<real>::min()) {
vp[0] = 0 + I1 / 3;
vp[1] = 0 + I1 / 3;
vp[2] = 0 + I1 / 3;
vp[0] = 0 + tr;
vp[1] = 0 + tr;
vp[2] = 0 + tr;
return;
}

// compute T
tmatrix<3u, 3u, real> S2;
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
S2(i, j) = 0;
for (int k = 0; k < 3; k++) {
S2(i, j) += S(i, k) * S(k, j);
}
}
}
const tmatrix<3u, 3u, real> T = S2 - (real{2} * J2 / 3) * eye;
const auto S2 = eval(S * S);
const tmatrix<3u, 3u, real> T = S2 - (one_third * real{2} * J2) * eye;

// compute d
real TmsS = 0;
Expand All @@ -96,29 +77,31 @@ namespace tfel::math::internals {
TpsS += (T(i, j) + s * S(i, j)) * (T(i, j) + s * S(i, j));
}
}
const real d = std::sqrt((TmsS) / (TpsS));
const auto d = std::sqrt((TmsS) / (TpsS));

const auto sj = [](const real value) -> int {
return (real(0) < value) - (value < real(0));
}(1 - d);

if (sj * (1 - d) < std::numeric_limits<real>::min()) {
vp[0] = std::sqrt(J2) + I1 / 3;
vp[1] = 0 + I1 / 3;
vp[2] = -vp[0] + I1 / 3;
vp[0] = std::sqrt(J2) + tr;
vp[1] = 0 + tr;
vp[2] = -vp[0] + tr;
return;
}

const auto alpha = real{2} / 3 * std::atan(std::pow(d, sj));
// compute alpha
const auto dsj = sj < 0 ? real{1} / d : d;
const auto alpha = 2 * one_third * std::atan(dsj);

// distinct eigenvalue
const auto cd = sj * s * std::cos(alpha);
vp[0] = 2 * cd + I1 / 3;
vp[0] = 2 * cd + tr;

// other eigenvalues
const auto sd = std::sqrt(J2) * std::sin(alpha);
vp[1] = -cd + sd + I1 / 3;
vp[2] = -cd - sd + I1 / 3;
vp[1] = -cd + sd + tr;
vp[2] = -cd - sd + tr;
}

template <typename real>
Expand All @@ -130,11 +113,7 @@ namespace tfel::math::internals {
const real D,
const real E,
const real F) {
constexpr auto one = real{1};
constexpr const auto one_third = one / 3;
constexpr const auto cste = Cste<real>::sqrt2;
// compute the deviatoric part of M : dev(M)
const auto trM_3 = one_third * (A + B + C);
const real s[6u] = {A, B, C, D * cste, E * cste, F * cste};
// computing eigen values
HarariEigensolver3x3::computeEigenValues(vp, A, B, C, D, E, F);
Expand Down
33 changes: 33 additions & 0 deletions tests/Math/stensor/stensor_eigenvectors3.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,9 +130,31 @@ void test2() {
vp2 = vp(1);
vp3 = vp(2);

vec1[0] = m(0, 0);
vec1[1] = m(1, 0);
vec1[2] = m(2, 0);

vec2[0] = m(0, 1);
vec2[1] = m(1, 1);
vec2[2] = m(2, 1);

vec3[0] = m(0, 2);
vec3[1] = m(1, 2);
vec3[2] = m(2, 2);

assert(abs(vp1 - static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vp2 - static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vp3 - static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());

assert(abs(vec1[0] - static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec1[1] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec1[2] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec2[0] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec2[1] - static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec2[2] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec3[0] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec3[1] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec3[2] - static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());
}

template <typename T>
Expand Down Expand Up @@ -181,6 +203,17 @@ void test3() {
assert(abs(vp3 - static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vp2 - static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vp1 - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());

assert(abs(vec1[0] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec1[1] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec1[2] - static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec2[0] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec2[1] - static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec2[2] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec3[0] + static_cast<T>(1)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec3[1] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());
assert(abs(vec3[2] - static_cast<T>(0)) < 20 * numeric_limits<T>::epsilon());

}

/* coverity [UNCAUGHT_EXCEPT]*/
Expand Down

0 comments on commit 1522c83

Please sign in to comment.