Commit 18bad8e9 authored by Mathias Goldau's avatar Mathias Goldau
Browse files

[CHANGE] Now jacobi method returns a sorted eigensystem

parent 57140131
......@@ -25,8 +25,10 @@
#ifndef WTENSORFUNCTIONS_H
#define WTENSORFUNCTIONS_H
#include <algorithm>
#include <cmath>
#include <complex>
#include <iostream>
#include <utility>
#include <vector>
......@@ -50,6 +52,33 @@ typedef boost::array< std::pair< double, WVector3D >, 3 > RealEigenSystem;
*/
typedef boost::array< std::pair< std::complex< double >, WVector3D >, 3 > EigenSystem;
std::ostream& operator<<( std::ostream& os, const RealEigenSystem& sys )
{
os << sys[0].first << ", " << sys[0].second << std::endl;
os << sys[1].first << ", " << sys[1].second << std::endl;
os << sys[2].first << ", " << sys[2].second << std::endl;
return os;
}
namespace
{
void sortRealEigenSystem( RealEigenSystem* es )
{
if( ( *es )[0].first > ( *es )[2].first )
{
std::swap( ( *es )[0], ( *es )[2] );
}
if( ( *es )[0].first > ( *es )[1].first )
{
std::swap( ( *es )[0], ( *es )[1] );
}
if( ( *es )[1].first > ( *es )[2].first )
{
std::swap( ( *es )[1], ( *es )[2] );
}
}
}
/**
* Compute all eigenvalues as well as the corresponding eigenvectors of a
* symmetric real Matrix.
......@@ -97,6 +126,7 @@ void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat, RealEigenSystem
result[i].second[j] = static_cast< double >( ev( j, i ) );
}
}
sortRealEigenSystem( es );
return;
}
......
......@@ -58,18 +58,10 @@ public:
t( 2, 2 ) = 0.00015239;
RealEigenSystem sys;
jacobiEigenvector3D( t, &sys );
WAssert( !wlimits::isnan( sys[0].first ), "Damn!" );
// Note: We don't use TS_ASSERT_DELTA here since its output is restricted to 4 digits after the point aka comma.
double delta = 1.0e-8;
if( std::abs( sys[0].first - 0.000196397 ) > delta ||
std::abs( sys[1].first - 0.000155074 ) > delta ||
std::abs( sys[2].first - 0.000150625 ) > delta )
{
TS_FAIL( "The eigenvalues are incorrect: " );
using string_utils::operator<<;
std::cout << std::fixed << std::setprecision( 16 ) << std::endl << WVector3D( sys[0].first, sys[1].first, sys[2].first ) << std::endl;
std::cout << "but got:" << std::endl << 0.000196397 << " " << 0.000155074 << " " << 0.00015062 << std::endl;
}
TS_ASSERT_DELTA( sys[0].first, 1.5062467240725114e-04, 1e-9 );
TS_ASSERT_DELTA( sys[1].first, 1.5507354000104679e-04, 1e-9 );
TS_ASSERT_DELTA( sys[2].first, 1.9639678759170208e-04, 1e-9 );
}
/**
......@@ -107,9 +99,9 @@ public:
jacobiEigenvector3D( t, &sys );
TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
TS_ASSERT_DELTA( sys[2].first, -3.0, 1e-6 );
TS_ASSERT_DELTA( sys[0].first, -3.0, 1e-6 );
TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 );
TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 );
TS_ASSERT_DELTA( sys[0].second.dotProduct( sys[1].second ), 0.0, 1e-9 );
TS_ASSERT_DELTA( sys[1].second.dotProduct( sys[2].second ), 0.0, 1e-9 );
TS_ASSERT_DELTA( sys[2].second.dotProduct( sys[0].second ), 0.0, 1e-9 );
......@@ -158,8 +150,8 @@ public:
jacobiEigenvector3D( t, &sys );
TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
TS_ASSERT_DELTA( sys[1].first, 0.0, 1e-6 );
TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 );
TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 );
TS_ASSERT_DELTA( sys[2].first, 1.0, 1e-6 );
TS_ASSERT_DELTA( sys[0].second.dotProduct( sys[1].second ), 0.0, 1e-9 );
TS_ASSERT_DELTA( sys[1].second.dotProduct( sys[2].second ), 0.0, 1e-9 );
......@@ -193,9 +185,9 @@ public:
jacobiEigenvector3D( t, &sys );
TS_ASSERT_DELTA( sys[0].first, 2.000001, 1e-6 );
TS_ASSERT_DELTA( sys[1].first, 0.0, 1e-6 );
TS_ASSERT_DELTA( sys[2].first, 1.999998, 1e-6 );
TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 );
TS_ASSERT_DELTA( sys[1].first, 1.999998, 1e-6 );
TS_ASSERT_DELTA( sys[2].first, 2.000001, 1e-6 );
TS_ASSERT_DELTA( sys[0].second.dotProduct( sys[1].second ), 0.0, 1e-9 );
TS_ASSERT_DELTA( sys[1].second.dotProduct( sys[2].second ), 0.0, 1e-9 );
TS_ASSERT_DELTA( sys[2].second.dotProduct( sys[0].second ), 0.0, 1e-9 );
......@@ -211,9 +203,9 @@ public:
jacobiEigenvector3D( t, &sys );
TS_ASSERT_DELTA( sys[0].first, 3.824572321236e30, 1e-6 );
TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 );
TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 );
TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
TS_ASSERT_DELTA( sys[2].first, 3.824572321236e30, 1e-6 );
TS_ASSERT_DELTA( sys[0].second.dotProduct( sys[1].second ), 0.0, 1e-9 );
TS_ASSERT_DELTA( sys[1].second.dotProduct( sys[2].second ), 0.0, 1e-9 );
TS_ASSERT_DELTA( sys[2].second.dotProduct( sys[0].second ), 0.0, 1e-9 );
......@@ -274,8 +266,8 @@ public:
jacobiEigenvector3D( t, &sys );
TS_ASSERT_DELTA( sys[0].first, 2.0, 1e-6 );
TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 );
TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 );
TS_ASSERT_DELTA( sys[0].second.dotProduct( sys[1].second ), 0.0, 1e-6 );
TS_ASSERT_DELTA( sys[1].second.dotProduct( sys[2].second ), 0.0, 1e-6 );
......
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