From b9e712ac2cda3347aab1ff14bca9b3df6727245f Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 6 Nov 2014 16:53:15 +0100
Subject: [PATCH] Rip out the complete test for eigenValuesNonSym and replace
it by the Rosser matrix test from fieldmatrixtest.cc
The old test was complete nonsense: an implementation of some sort of strange
undocumented shape function set that happened to need eigenvalues for something.
Lots of undocumented code, and no checks at all whether the values obtained
by eigenValuesNonSym actually make sense.
Only, for the particular choice of parameters used, eigenValuesNonSym was never
called. Yes, you read this correctly: these hundreds of lines of code didn't
actually test anything.
The new test is copied form fieldmatrixtest.cc and adapted for DynamicMatrix
types. It computes eigenvalues of a given 8x8 matrix and compares with the known
results. Unfortunately, the test matrix is only symmetric. Tests with more
general matrices are left for another day.
---
dune/common/test/eigenvaluestest.cc | 310 +++++++---------------------
1 file changed, 76 insertions(+), 234 deletions(-)
diff --git a/dune/common/test/eigenvaluestest.cc b/dune/common/test/eigenvaluestest.cc
index 1b89210fc..e0e92c7ae 100644
--- a/dune/common/test/eigenvaluestest.cc
+++ b/dune/common/test/eigenvaluestest.cc
@@ -4,257 +4,99 @@
#include "config.h"
#endif
-//==============================================================================
-//!
-//! \date Nov 9 2011
-//!
-//! \author Arne Morten Kvarving / SINTEF
-//!
-//==============================================================================
-
#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
#include <dune/common/dynmatrixev.hh>
#include <algorithm>
#include <complex>
+using namespace Dune;
-//! \brief Represents a cardinal function on a line
-template<class ctype, class rtype>
-class LagrangeCardinalFunction
-{
-public:
- //! \brief Empty default constructor
- LagrangeCardinalFunction() {}
-
- //! \brief Construct a cardinal function with the given nodes
- //! \param[in] nodes_ The nodes
- //! \param[in] i The node this function is associated with
- LagrangeCardinalFunction(const std::vector<rtype>& nodes_,
- size_t i)
- : nodes(nodes_), node(i) {}
-
- //! \brief Evaluate the shape function
- //! \param[in] local The local coordinates
- rtype evaluateFunction(const ctype& local) const
- {
- rtype result = 1;
- for (size_t i=0; i < nodes.size(); ++i) {
- if (i != node)
- result *= (local-nodes[i])/(nodes[node]-nodes[i]);
- }
-
- return result;
- }
-
- //! \brief Evaluate the derivative of the cardinal function
- //! \param[in] local The local coordinates
- rtype evaluateGradient(const ctype& local) const
- {
- rtype result = 0;
- for (size_t i=0; i < nodes.size(); ++i) {
- rtype f = 1;
- for (int j=0; j < nodes.size(); ++j) {
- if (i != j && j != node)
- f *= (local-nodes[j])/(nodes[node]-nodes[j]);
- }
- result += f/(nodes[node]-nodes[i]);
- }
-
- return result;
- }
+/** \brief Test the eigenvalue code with the Rosser test matrix
-private:
- //! \brief The nodes
- std::vector<rtype> nodes;
+ This matrix was a challenge for many matrix eigenvalue
+ algorithms. But the Francis QR algorithm, as perfected by
+ Wilkinson and implemented in EISPACK, has no trouble with it. The
+ matrix is 8-by-8 with integer elements. It has:
- size_t node;
-};
+ * A double eigenvalue
+ * Three nearly equal eigenvalues
+ * Dominant eigenvalues of opposite sign
+ * A zero eigenvalue
+ * A small, nonzero eigenvalue
-//! \brief Represents a tensor-product of 1D functions
-template<class rtype, class ctype, class ftype, int dim>
-class TensorProductFunction
+*/
+template<typename ft>
+void testRosserMatrix()
{
-public:
- //! \brief The dimension of the function
- enum { dimension = dim };
-
- //! \brief Empty default constructor
- TensorProductFunction() {}
-
- //! \brief Construct a tensor-product function
- //! \param[in] funcs_ The functions
- TensorProductFunction(const Dune::FieldVector<ftype, dim>& funcs_)
- : funcs(funcs_) {}
-
- //! \brief Evaluate the function
- //! \param[in] local The local coordinates
- rtype evaluateFunction(const Dune::FieldVector<ctype,dim>& local) const
+ // Hack: I want this matrix to be a DynamicMatrix, but currently I cannot
+ // initialize such a matrix from initializer lists. Therefore I take the
+ // detour over a FieldMatrix.
+ FieldMatrix<ft,8,8> AField = {
+ { 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}
+ };
+
+ DynamicMatrix<ft> A(8,8);
+ for (int i=0; i<8; i++)
+ for (int j=0; j<8; j++)
+ A[i][j] = AField[i][j];
+
+ // compute eigenvalues
+ DynamicVector<std::complex<ft> > eigenComplex;
+ DynamicMatrixHelp::eigenValuesNonSym(A, eigenComplex);
+
+ // test results
+ /*
+ reference solution computed with octave 3.2
+
+ > format long e
+ > eig(rosser())
+
+ */
+ std::vector<ft> reference = {
+ -1.02004901843000e+03,
+ -4.14362871168386e-14,
+ 9.80486407214362e-02,
+ 1.00000000000000e+03,
+ 1.00000000000000e+03,
+ 1.01990195135928e+03,
+ 1.02000000000000e+03,
+ 1.02004901843000e+03
+ };
+
+ std::vector<ft> eigenRealParts(8);
+ for (int i=0; i<8; i++)
+ eigenRealParts[i] = std::real(eigenComplex[i]);
+
+ std::sort(eigenRealParts.begin(), eigenRealParts.end());
+
+ for (int i=0; i<8; i++)
{
- rtype result = 1;
- for (int i=0; i < dim; ++i)
- result *= funcs[i].evaluateFunction(local[i]);
+ if (std::fabs(std::imag(eigenComplex[i])) > 1e-10)
+ DUNE_THROW(MathError, "Symmetric matrix has complex eigenvalue");
- return result;
+ if( std::fabs(reference[i] - eigenRealParts[i]) > 1e-10 )
+ DUNE_THROW(MathError,"error computing eigenvalues");
}
- Dune::FieldVector<rtype, dim>
- evaluateGradient(const Dune::FieldVector<ctype,dim>& local) const
- {
- Dune::FieldVector<rtype, dim> result;
- for (int i=0; i < dim; ++i)
- result[i] = funcs[i].evaluateGradient(local[i]);
- }
-private:
- Dune::FieldVector<ftype, dim> funcs;
-};
-
-template<int dim>
-class PNShapeFunctionSet
-{
-public:
- typedef LagrangeCardinalFunction<double, double> CardinalFunction;
-
- typedef TensorProductFunction<double, double, CardinalFunction, dim>
- ShapeFunction;
-
- PNShapeFunctionSet(int n1, int n2, int n3=0)
- {
- int dims[3] = {n1, n2, n3};
- cfuncs.resize(dim);
- for (int i=0; i < dim; ++i) {
- std::vector<double> grid;
- grid = gaussLobattoLegendreGrid(dims[i]);
- for (int j=0; j<dims[i]; ++j)
- cfuncs[i].push_back(CardinalFunction(grid,j));
- }
- int l=0;
- Dune::FieldVector<CardinalFunction,dim> fs;
- if (dim == 3) {
- f.resize(n1*n2*n3);
- for (int k=0; k < n3; ++k) {
- for (int j=0; j < n2; ++j)
- for (int i=0; i< n1; ++i) {
- fs[0] = cfuncs[0][i];
- fs[1] = cfuncs[1][j];
- fs[2] = cfuncs[2][k];
- f[l++] = ShapeFunction(fs);
- }
- }
- } else {
- f.resize(n1*n2);
- for (int j=0; j < n2; ++j) {
- for (int i=0; i< n1; ++i) {
- fs[0] = cfuncs[0][i];
- fs[1] = cfuncs[1][j];
- f[l++] = ShapeFunction(fs);
- }
- }
- }
- }
-
- //! \brief Obtain a given shape function
- //! \param[in] i The requested shape function
- const ShapeFunction& operator[](int i) const
- {
- return f[i];
- }
-
- int size()
- {
- return f.size();
- }
-protected:
- std::vector< std::vector<CardinalFunction> > cfuncs;
- std::vector<ShapeFunction> f;
-
- double legendre(double x, int n)
- {
- std::vector<double> Ln;
- Ln.resize(n+1);
- Ln[0] = 1.f;
- Ln[1] = x;
- if( n > 1 ) {
- for( int i=1; i<n; i++ )
- Ln[i+1] = (2*i+1.0)/(i+1.0)*x*Ln[i]-i/(i+1.0)*Ln[i-1];
- }
-
- return Ln[n];
- }
-
- double legendreDerivative(double x, int n)
- {
- std::vector<double> Ln;
- Ln.resize(n+1);
-
- Ln[0] = 1.0; Ln[1] = x;
-
- if( (x == 1.0) || (x == -1.0) )
- return( pow(x,n-1)*n*(n+1.0)/2.0 );
- else {
- for( int i=1; i<n; i++ )
- Ln[i+1] = (2.0*i+1.0)/(i+1.0)*x*Ln[i]-(double)i/(i+1.0)*Ln[i-1];
- return( (double)n/(1.0-x*x)*Ln[n-1]-n*x/(1-x*x)*Ln[n] );
- }
- }
-
-
- std::vector<double> gaussLegendreGrid(int n)
- {
- Dune::DynamicMatrix<double> A(n,n,0.0);
-
- A[0][1] = 1.f;
- for (int i=1; i<n-1; ++i) {
- A[i][i-1] = (double)i/(2.0*(i+1.0)-1.0);
- A[i][i+1] = (double)(i+1.0)/(2*(i+1.0)-1.0);
- }
- A[n-1][n-2] = (n-1.0)/(2.0*n-1.0);
-
- Dune::DynamicVector<std::complex<double> > eigenValues(n);
- Dune::DynamicMatrixHelp::eigenValuesNonSym(A, eigenValues);
-
- std::vector<double> result(n);
- for (int i=0; i < n; ++i)
- result[i] = std::real(eigenValues[i]);
- std::sort(result.begin(),result.begin()+n);
-
- return result;
- }
-
- std::vector<double> gaussLobattoLegendreGrid(int n)
- {
- assert(n > 1);
- const double tolerance = 1.e-15;
-
- std::vector<double> result(n);
- result[0] = 0.0;
- result[n-1] = 1.0;
- if (n == 3)
- result[1] = 0.5;
-
- if (n < 4)
- return result;
-
- std::vector<double> glgrid = gaussLegendreGrid(n-1);
- for (int i=1; i<n-1; ++i) {
- result[i] = (glgrid[i-1]+glgrid[i])/2.0;
- double old = 0.0;
- while (std::abs(old-result[i]) > tolerance) {
- old = result[i];
- double L = legendre(old, n-1);
- double Ld = legendreDerivative(old, n-1);
- result[i] += (1.0-old*old)*Ld/((n-1.0)*n*L);
- }
- result[i] = (result[i]+1.0)/2.0;
- }
+ std::cout << "Eigenvalues of Rosser matrix: " << eigenComplex << std::endl;
+}
- return result;
- }
-};
-int main()
+int main() try
{
- // Not really a test but better than nothing.
- PNShapeFunctionSet<2> lbasis(2, 3);
+ testRosserMatrix<double>();
return 0;
+} catch (Exception exception)
+{
+ std::cerr << exception << std::endl;
+ return 1;
}
--
GitLab