Capabilities to compute eigenvalues on non-symmetric matrices
Code contributed by Arne Morten Kvarving.

[[Imported from SVN: r7061]]

COPYING

@@ -13,6 +13,7 @@ Copyright holders:
 2010--2012    Christoph Grüninger
 2006          Bernhard Haasdonk
 2012          Olaf Ippisch
+2012          Arne Morten Kvarving
 2009          Leonard Kern
 2005--2007    Sreejith Pulloor Kuttanikkad
 2003--2012    Robert Klöfkorn

dune/common/Makefile.am

@@ -8,6 +8,7 @@ noinst_LTLIBRARIES = libcommon.la
 libcommon_la_SOURCES =				\
 	debugallocator.cc			\
 	fmatrixev.cc \
+	dynmatrixev.cc                           \
 	ios_state.cc				\
 	parametertree.cc                        \
 	parametertreeparser.cc			\

dune/common/dynmatrixev.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FMATRIXEIGENVALUES_EXT_CC
+#define DUNE_FMATRIXEIGENVALUES_EXT_CC
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <cmath>
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+
+#if HAVE_LAPACK
+
+// nonsymetric matrices
+#define DGEEV_FORTRAN FC_FUNC (dgeev, DGEEV)
+
+// dsyev declaration (in liblapack)
+extern "C" {
+
+  /*
+   *
+   **  purpose
+   **  =======
+   **
+   **  xgeev computes for an N-by-N BASE DATA TYPE nonsymmetric matrix A, the
+   **  eigenvalues and, optionally, the left and/or right eigenvectors.
+   **
+   **  The right eigenvector v(j) of A satisfies
+   **                   A * v(j) = lambda(j) * v(j)
+   **  where lambda(j) is its eigenvalue.
+   **  The left eigenvector u(j) of A satisfies
+   **                u(j)**T * A = lambda(j) * u(j)**T
+   **  where u(j)**T denotes the transpose of u(j).
+   **
+   **  The computed eigenvectors are normalized to have Euclidean norm
+   **  equal to 1 and largest component real.
+   **
+   **  arguments
+   **  =========
+   **
+   **  jobvl   (input) char
+   **          = 'n': left eigenvectors of a are not computed;
+   **          = 'v': left eigenvectors of a are computed.
+   **
+   **  jobvr   (input) char
+   **          = 'n': right eigenvectors of a are not computed;
+   **          = 'v': right eigenvectors of a are computed.
+   **
+   **  n       (input) long int
+   **          the order of the matrix v. v >= 0.
+   **
+   **  a       (input/output) BASE DATA TYPE array, dimension (lda,n)
+   **          on entry, the n-by-n matrix a.
+   **          on exit, a has been overwritten.
+   **
+   **  lda     (input) long int
+   **          the leading dimension of the array a.  lda >= max(1,n).
+   **
+   **  wr      (output) BASE DATA TYPE array, dimension (n)
+   **  wi      (output) BASE DATA TYPE array, dimension (n)
+   **          wr and wi contain the real and imaginary parts,
+   **          respectively, of the computed eigenvalues.  complex
+   **          conjugate pairs of eigenvalues appear consecutively
+   **          with the eigenvalue having the positive imaginary part
+   **          first.
+   **
+   **  vl      (output) COMPLEX DATA TYPE array, dimension (ldvl,n)
+   **          if jobvl = 'v', the left eigenvectors u(j) are stored one
+   **          after another in the columns of vl, in the same order
+   **          as their eigenvalues.
+   **          if jobvl = 'n', vl is not referenced.
+   **          if the j-th eigenvalue is real, then u(j) = vl(:,j),
+   **          the j-th column of vl.
+   **          if the j-th and (j+1)-st eigenvalues form a complex
+   **          conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and
+   **          u(j+1) = vl(:,j) - i*vl(:,j+1).
+   **
+   **  ldvl    (input) long int
+   **          the leading dimension of the array vl.  ldvl >= 1; if
+   **          jobvl = 'v', ldvl >= n.
+   **
+   **  vr      (output) COMPLEX DATA TYPE array, dimension (ldvr,n)
+   **          if jobvr = 'v', the right eigenvectors v(j) are stored one
+   **          after another in the columns of vr, in the same order
+   **          as their eigenvalues.
+   **          if jobvr = 'n', vr is not referenced.
+   **          if the j-th eigenvalue is real, then v(j) = vr(:,j),
+   **          the j-th column of vr.
+   **          if the j-th and (j+1)-st eigenvalues form a complex
+   **          conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and
+   **          v(j+1) = vr(:,j) - i*vr(:,j+1).
+   **
+   **  ldvr    (input) long int
+   **          the leading dimension of the array vr.  ldvr >= 1; if
+   **          jobvr = 'v', ldvr >= n.
+   **
+   **  work    (workspace/output) BASE DATA TYPE array, dimension (max(1,lwork))
+   **          on exit, if info = 0, work(1) returns the optimal lwork.
+   **
+   **  lwork   (input) long int
+   **          the dimension of the array work.  lwork >= max(1,3*n), and
+   **          if jobvl = 'v' or jobvr = 'v', lwork >= 4*n.  for good
+   **          performance, lwork must generally be larger.
+   **
+   **          if lwork = -1, then a workspace query is assumed; the routine
+   **          only calculates the optimal size of the work array, returns
+   **          this value as the first entry of the work array, and no error
+   **          message related to lwork is issued by xerbla.
+   **
+   **  info    (output) long int
+   **          = 0:  successful exit
+   **          < 0:  if info = -i, the i-th argument had an illegal value.
+   **          > 0:  if info = i, the qr algorithm failed to compute all the
+   **                eigenvalues, and no eigenvectors have been computed;
+   **                elements i+1:n of wr and wi contain eigenvalues which
+   **                have converged.
+   **
+   **/
+
+  extern void DGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long
+                            int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+                            const long int* ldvl, double* vr, const long int* ldvr, double* work,
+                            const long int* lwork, const long int* info);
+
+} // end extern C
+#endif
+
+namespace Dune {
+
+  namespace DynamicMatrixHelp {
+
+    void eigenValuesNonsymLapackCall(
+      const char* jobvl, const char* jobvr, const long
+      int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+      const long int* ldvl, double* vr, const long int* ldvr, double* work,
+      const long int* lwork, const long int* info)
+    {
+#if HAVE_LAPACK
+      // call LAPACK dgeev
+      DGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr,
+                    work, lwork, info);
+#else
+      DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!");
+#endif
+    }
+
+  } // end namespace FMatrixHelp
+
+} // end namespace Dune
+#endif

dune/common/dynmatrixev.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_DYNMATRIXEIGENVALUES_HH
+#define DUNE_DYNMATRIXEIGENVALUES_HH
+
+#include "dynmatrix.hh"
+
+namespace Dune {
+
+  namespace DynamicMatrixHelp {
+
+    // defined in fmatrixev_ext.cpp
+    extern void eigenValuesNonsymLapackCall(
+      const char* jobvl, const char* jobvr, const long
+      int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+      const long int* ldvl, double* vr, const long int* ldvr, double* work,
+      const long int* lwork, const long int* info);
+
+    template <typename K, class C>
+    static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
+                                  DynamicVector<C>& eigenValues)
+    {
+      {
+        const long int N = matrix.rows();
+        const char jobvl = 'n';
+        const char jobvr = 'n';
+
+
+        // matrix to put into dgeev
+        double matrixVector[N * N];
+
+        // copy matrix
+        int row = 0;
+        for(int i=0; i<N; ++i)
+        {
+          for(int j=0; j<N; ++j, ++row)
+          {
+            matrixVector[ row ] = matrix[ i ][ j ];
+          }
+        }
+
+        // working memory
+        double eigenR[N];
+        double eigenI[N];
+        double work[3*N];
+
+        // return value information
+        long int info = 0;
+        long int lwork = 3*N;
+
+        // call LAPACK routine (see fmatrixev_ext.cc)
+        eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
+                                    &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
+                                    &lwork, &info);
+
+        if( info != 0 )
+        {
+          std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
+          DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
+        }
+        for (int i=0; i<N; ++i)
+          eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
+      }
+    }
+
+  }
+
+}
+
+#endif

dune/common/fmatrixev.cc

@@ -78,6 +78,111 @@ extern "C" {
                             int* n, double* a, const long int* lda, double* w,
                             double* work, const long int* lwork, long int* info);
 
+  /*
+   *
+   **  purpose
+   **  =======
+   **
+   **  xgeev computes for an N-by-N BASE DATA TYPE nonsymmetric matrix A, the
+   **  eigenvalues and, optionally, the left and/or right eigenvectors.
+   **
+   **  The right eigenvector v(j) of A satisfies
+   **                   A * v(j) = lambda(j) * v(j)
+   **  where lambda(j) is its eigenvalue.
+   **  The left eigenvector u(j) of A satisfies
+   **                u(j)**T * A = lambda(j) * u(j)**T
+   **  where u(j)**T denotes the transpose of u(j).
+   **
+   **  The computed eigenvectors are normalized to have Euclidean norm
+   **  equal to 1 and largest component real.
+   **
+   **  arguments
+   **  =========
+   **
+   **  jobvl   (input) char
+   **          = 'n': left eigenvectors of a are not computed;
+   **          = 'v': left eigenvectors of a are computed.
+   **
+   **  jobvr   (input) char
+   **          = 'n': right eigenvectors of a are not computed;
+   **          = 'v': right eigenvectors of a are computed.
+   **
+   **  n       (input) long int
+   **          the order of the matrix v. v >= 0.
+   **
+   **  a       (input/output) BASE DATA TYPE array, dimension (lda,n)
+   **          on entry, the n-by-n matrix a.
+   **          on exit, a has been overwritten.
+   **
+   **  lda     (input) long int
+   **          the leading dimension of the array a.  lda >= max(1,n).
+   **
+   **  wr      (output) BASE DATA TYPE array, dimension (n)
+   **  wi      (output) BASE DATA TYPE array, dimension (n)
+   **          wr and wi contain the real and imaginary parts,
+   **          respectively, of the computed eigenvalues.  complex
+   **          conjugate pairs of eigenvalues appear consecutively
+   **          with the eigenvalue having the positive imaginary part
+   **          first.
+   **
+   **  vl      (output) COMPLEX DATA TYPE array, dimension (ldvl,n)
+   **          if jobvl = 'v', the left eigenvectors u(j) are stored one
+   **          after another in the columns of vl, in the same order
+   **          as their eigenvalues.
+   **          if jobvl = 'n', vl is not referenced.
+   **          if the j-th eigenvalue is real, then u(j) = vl(:,j),
+   **          the j-th column of vl.
+   **          if the j-th and (j+1)-st eigenvalues form a complex
+   **          conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and
+   **          u(j+1) = vl(:,j) - i*vl(:,j+1).
+   **
+   **  ldvl    (input) long int
+   **          the leading dimension of the array vl.  ldvl >= 1; if
+   **          jobvl = 'v', ldvl >= n.
+   **
+   **  vr      (output) COMPLEX DATA TYPE array, dimension (ldvr,n)
+   **          if jobvr = 'v', the right eigenvectors v(j) are stored one
+   **          after another in the columns of vr, in the same order
+   **          as their eigenvalues.
+   **          if jobvr = 'n', vr is not referenced.
+   **          if the j-th eigenvalue is real, then v(j) = vr(:,j),
+   **          the j-th column of vr.
+   **          if the j-th and (j+1)-st eigenvalues form a complex
+   **          conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and
+   **          v(j+1) = vr(:,j) - i*vr(:,j+1).
+   **
+   **  ldvr    (input) long int
+   **          the leading dimension of the array vr.  ldvr >= 1; if
+   **          jobvr = 'v', ldvr >= n.
+   **
+   **  work    (workspace/output) BASE DATA TYPE array, dimension (max(1,lwork))
+   **          on exit, if info = 0, work(1) returns the optimal lwork.
+   **
+   **  lwork   (input) long int
+   **          the dimension of the array work.  lwork >= max(1,3*n), and
+   **          if jobvl = 'v' or jobvr = 'v', lwork >= 4*n.  for good
+   **          performance, lwork must generally be larger.
+   **
+   **          if lwork = -1, then a workspace query is assumed; the routine
+   **          only calculates the optimal size of the work array, returns
+   **          this value as the first entry of the work array, and no error
+   **          message related to lwork is issued by xerbla.
+   **
+   **  info    (output) long int
+   **          = 0:  successful exit
+   **          < 0:  if info = -i, the i-th argument had an illegal value.
+   **          > 0:  if info = i, the qr algorithm failed to compute all the
+   **                eigenvalues, and no eigenvectors have been computed;
+   **                elements i+1:n of wr and wi contain eigenvalues which
+   **                have converged.
+   **
+   **/
+
+  extern void DGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long
+                            int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+                            const long int* ldvl, double* vr, const long int* ldvr, double* work,
+                            const long int* lwork, const long int* info);
+
 } // end extern C
 #endif
 
@@ -98,6 +203,21 @@ namespace Dune {
 #endif
     }
 
+    void eigenValuesNonsymLapackCall(
+      const char* jobvl, const char* jobvr, const long
+      int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+      const long int* ldvl, double* vr, const long int* ldvr, double* work,
+      const long int* lwork, const long int* info)
+    {
+#if HAVE_LAPACK
+      // call LAPACK dgeev
+      DGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr,
+                    work, lwork, info);
+#else
+      DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!");
+#endif
+    }
+
   } // end namespace FMatrixHelp
 
 } // end namespace Dune

dune/common/fmatrixev.hh

@@ -30,6 +30,12 @@ namespace Dune {
       int* n, double* a, const long int* lda, double* w,
       double* work, const long int* lwork, long int* info);
 
+    extern void eigenValuesNonsymLapackCall(
+      const char* jobvl, const char* jobvr, const long
+      int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
+      const long int* ldvl, double* vr, const long int* ldvr, double* work,
+      const long int* lwork, const long int* info);
+
     /** \brief calculates the eigenvalues of a symetric field matrix
         \param[in]  matrix matrix eigenvalues are calculated for
         \param[out] eigenvalues FieldVector that contains eigenvalues in
@@ -122,6 +128,58 @@ namespace Dune {
       }
     }
 
+    template <int dim, typename K, class C>
+    static void eigenValuesNonSym(const FieldMatrix<K, dim, dim>& matrix,
+                                  FieldVector<C, dim>& eigenValues)
+    {
+      {
+        const long int N = dim ;
+        const char jobvl = 'n';
+        const char jobvr = 'n';
+
+        // length of matrix vector
+        const long int w = N * N ;
+
+        // matrix to put into dgeev
+        double matrixVector[dim * dim];
+
+        // copy matrix
+        int row = 0;
+        for(int i=0; i<dim; ++i)
+        {
+          for(int j=0; j<dim; ++j, ++row)
+          {
+            matrixVector[ row ] = matrix[ i ][ j ];
+          }
+        }
+
+        // working memory
+        double eigenR[dim];
+        double eigenI[dim];
+        double work[3*dim];
+
+        // return value information
+        long int info = 0;
+        long int lwork = 3*dim;
+
+        // call LAPACK routine (see fmatrixev_ext.cc)
+        eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
+                                    &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
+                                    &lwork, &info);
+
+        if( info != 0 )
+        {
+          std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
+          DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
+        }
+        for (int i=0; i<N; ++i) {
+          eigenValues[i].real = eigenR[i];
+          eigenValues[i].imag = eigenI[i];
+        }
+      }
+
+    }
+
   } // end namespace FMatrixHelp
 
   /** @} end documentation */

dune/common/test/Makefile.am

@@ -11,6 +11,7 @@ TESTPROGS = \
     diagonalmatrixtest \
     dynmatrixtest \
     dynvectortest \
+    eigenvaluestest \
     enumsettest \
     fassigntest \
     fmatrixtest \
@@ -161,6 +162,9 @@ dynmatrixtest_SOURCES = dynmatrixtest.cc
 
 dynvectortest_SOURCES = dynvectortest.cc
 
+eigenvaluestest_SOURCES = eigenvaluestest.cc
+eigenvaluestest_LDADD = $(LAPACK_LIBS) $(LDADD) $(BLAS_LIBS) $(LIBS) $(FLIBS)
+
 fmatrixtest_SOURCES = fmatrixtest.cc
 fmatrixtest_LDADD = $(LAPACK_LIBS) $(LDADD) $(BLAS_LIBS) $(LIBS) $(FLIBS)
 

dune/common/test/eigenvaluestest.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+//==============================================================================
+//!
+//! \file shapefunctions.hpp
+//!
+//! \date Nov 9 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Classes for shape functions. Loosely based on code in dune-grid-howto
+//!
+//==============================================================================
+
+#include <dune/common/fvector.hh>
+#include <dune/common/dynmatrixev.hh>
+
+#include <complex>
+
+
+//! \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;
+  }
+
+private:
+  //! \brief The nodes
+  std::vector<rtype> nodes;
+
+  size_t node;
+};
+
+//! \brief Represents a tensor-product of 1D functions
+template<class rtype, class ctype, class ftype, int dim>
+class TensorProductFunction
+{
+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
+  {
+    rtype result = 1;
+    for (int i=0; i < dim; ++i)
+      result *= funcs[i].evaluateFunction(local[i]);
+
+    return result;
+  }
+
+  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;
+    }
+
+    return result;
+  }
+};
+
+int main()
+{
+  // Not really a test but better than nothing.
+  PNShapeFunctionSet<2> lbasis(2, 3);
+  return 0;
+}