Fix bug in MatrixFunctions for matrices with multiple eigenvalues.

Store indices, not eigenvalues, in clusters.
Bug was introduced in changeset a3a55357db7394281c872e911f13d69aba510aec
.

Eigen/src/Core/util/XprHelper.h

@@ -307,9 +307,9 @@ struct transfer_constness
   *
   * Example. Suppose that a, b, and c are of type Matrix3d. The user forms the expression a*(b+c).
   * b+c is an expression "sum of matrices", which we will denote by S. In order to determine how to nest it,
-  * the Product expression uses: nested<S, 3>::ret, which turns out to be Matrix3d because the internal logic of
+  * the Product expression uses: nested<S, 3>::type, which turns out to be Matrix3d because the internal logic of
   * nested determined that in this case it was better to evaluate the expression b+c into a temporary. On the other hand,
-  * since a is of type Matrix3d, the Product expression nests it as nested<Matrix3d, 3>::ret, which turns out to be
+  * since a is of type Matrix3d, the Product expression nests it as nested<Matrix3d, 3>::type, which turns out to be
   * const Matrix3d&, because the internal logic of nested determined that since a was already a matrix, there was no point
   * in copying it into another matrix.
   */

unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h

@@ -105,10 +105,10 @@ MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A)
   * \returns Iterator to cluster containing \p key, or \c clusters.end() if no cluster in \p m_clusters
   * contains \p key.
   */
-template <typename Scalar, typename ListOfClusters>
+template <typename Index, typename ListOfClusters>
-typename ListOfClusters::iterator matrix_function_find_cluster(Scalar key, ListOfClusters& clusters)
+typename ListOfClusters::iterator matrix_function_find_cluster(Index key, ListOfClusters& clusters)
 {
-  typename std::list<Scalar>::iterator j;
+  typename std::list<Index>::iterator j;
   for (typename ListOfClusters::iterator i = clusters.begin(); i != clusters.end(); ++i) {
     j = std::find(i->begin(), i->end(), key);
     if (j != i->end())
@@ -133,11 +133,11 @@ void matrix_function_partition_eigenvalues(const EivalsType& eivals, std::list<C
 {
   typedef typename EivalsType::Index Index;
   for (Index i=0; i<eivals.rows(); ++i) {
-    // Find set containing eivals(i), adding a new set if necessary
+    // Find cluster containing i-th ei'val, adding a new cluster if necessary
-    typename std::list<Cluster>::iterator qi = matrix_function_find_cluster(eivals(i), clusters);
+    typename std::list<Cluster>::iterator qi = matrix_function_find_cluster(i, clusters);
     if (qi == clusters.end()) {
       Cluster l;
-      l.push_back(eivals(i));
+      l.push_back(i);
       clusters.push_back(l);
       qi = clusters.end();
       --qi;
@@ -146,10 +146,10 @@ void matrix_function_partition_eigenvalues(const EivalsType& eivals, std::list<C
     // Look for other element to add to the set
     for (Index j=i+1; j<eivals.rows(); ++j) {
       if (abs(eivals(j) - eivals(i)) <= matrix_function_separation
-          && std::find(qi->begin(), qi->end(), eivals(j)) == qi->end()) {
+          && std::find(qi->begin(), qi->end(), j) == qi->end()) {
-        typename std::list<Cluster>::iterator qj = matrix_function_find_cluster(eivals(j), clusters);
+        typename std::list<Cluster>::iterator qj = matrix_function_find_cluster(j, clusters);
         if (qj == clusters.end()) {
-          qi->push_back(eivals(j));
+          qi->push_back(j);
         } else {
           qi->insert(qi->end(), qj->begin(), qj->end());
           clusters.erase(qj);
@@ -192,7 +192,7 @@ void matrix_function_compute_map(const EivalsType& eivals, const ListOfClusters&
   Index clusterIndex = 0;
   for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) {
     for (Index i = 0; i < eivals.rows(); ++i) {
-      if (std::find(cluster->begin(), cluster->end(), eivals(i)) != cluster->end()) {
+      if (std::find(cluster->begin(), cluster->end(), i) != cluster->end()) {
         eivalToCluster[i] = clusterIndex;
       }
     }
@@ -435,7 +435,7 @@ struct matrix_function_compute<MatrixType, 1>
     MatrixType U = schurOfA.matrixU();
 
     // partition eigenvalues into clusters of ei'vals "close" to each other
-    std::list<std::list<Scalar> > clusters; 
+    std::list<std::list<Index> > clusters; 
     matrix_function_partition_eigenvalues(T.diagonal(), clusters);
 
     // compute size of each cluster
@@ -446,7 +446,7 @@ struct matrix_function_compute<MatrixType, 1>
     Matrix<Index, Dynamic, 1> blockStart; 
     matrix_function_compute_block_start(clusterSize, blockStart);
 
-    // compute map so that eivalToCluster[i] = j means that ei'val T(i,i) is in j-th cluster 
+    // compute map so that eivalToCluster[i] = j means that i-th ei'val is in j-th cluster 
     Matrix<Index, Dynamic, 1> eivalToCluster;
     matrix_function_compute_map(T.diagonal(), clusters, eivalToCluster);
 

unsupported/Eigen/src/MatrixFunctions/StemFunction.h

@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed