EigenValues: added sorting option to verify decreasing order, added hacky matlab-like sort method in VectorT

core/algebra/EigValues.cpp

@@ -117,7 +117,7 @@ EVArnoldi::getEigenvalues ( const GenericMatrix & data, Vector & eigenvalues,
       cerr << "EVArnoldi: [" << iteration << "] delta=" << delta << endl;
   }
   eigenvectors = rmatrix;
-
+  
   for ( uint i = 0; i < k; i++ )
   {
     NICE::Vector tmp;
@@ -125,4 +125,47 @@ EVArnoldi::getEigenvalues ( const GenericMatrix & data, Vector & eigenvalues,
     data.multiply ( tmp, eigenvectors.getColumn ( i ) );
     eigenvalues[i] = tmp.scalarProduct ( eigenvectors.getColumn ( i ) );
   }
+  
+  
+  
+  // post-processing: ensure that eigenvalues are in correct order!
+  
+  if ( this->b_verifyDecreasingOrder && (k > 1) )
+  {
+    NICE::VectorT< int > ewPermutation;
+    eigenvalues.sortDescend ( ewPermutation );
+    
+    NICE::Vector tmpRow;
+    int tmpIdx (-1);
+    for ( uint i = 0; i < k; i++ )
+    { 
+      if ( i == ewPermutation[i] )
+      {
+        //identity - nothing to do here
+      }
+      else
+      {
+        if ( tmpIdx == -1 )
+        {
+          tmpIdx = i;
+          tmpRow = eigenvectors.getColumn ( i );
+          eigenvectors.getColumnRef ( i ) = eigenvectors.getColumn ( ewPermutation[i] );
+        }
+        else
+        {
+          if ( tmpIdx != ewPermutation[i] )
+          {
+            eigenvectors.getColumnRef ( i ) = eigenvectors.getColumn ( ewPermutation[i] );
+          }
+          else // tmpIdx == ewPermutation[i]
+          {
+            eigenvectors.getColumnRef ( i ) = tmpRow;
+            tmpIdx = -1;
+          }
+        }
+      }
+    }
+    
+  }// sorting is only useful if we compute more then 1 ew
+  
 }

core/algebra/EigValues.h

@@ -1,7 +1,7 @@
 /**
 * @file EigValues.h
 * @brief Computing eigenvalues and eigenvectors
-* @author Michael Koch,Erik Rodner
+* @author Michael Koch,Erik Rodner, Alexander Freytag
 * @date 08/19/2008
 
 */
@@ -44,12 +44,14 @@ class EVArnoldi : public EigValues
     uint maxiterations;
     double mindelta;
     bool verbose;
+    bool b_verifyDecreasingOrder;
 
   public:
-    EVArnoldi ( bool verbose = false, uint _maxiterations = 100, double _mindelta = 0.01 )
+    EVArnoldi ( bool verbose = false, uint _maxiterations = 100, double _mindelta = 0.01, bool b_verifyDecreasingOrder = true)
         : maxiterations ( _maxiterations ), mindelta ( _mindelta )
     {
       this->verbose = verbose;
+      this->b_verifyDecreasingOrder = b_verifyDecreasingOrder;
     };
 
     /**

core/vector/VectorT.h

@@ -390,7 +390,7 @@ public:
     inline void sortDescend();
 	 
     /**
-    * @brief sort elements in an descending order.
+    * @brief sort elements in an descending order. Permutation is only correct if all elements are different.
     */
     inline void sortDescend(VectorT<int> & permutation);
 

core/vector/VectorT.tcc

@@ -740,21 +740,32 @@ void VectorT<ElementType>::sortDescend() {
 }
 
 /**
- * @brief sort elements in an descending order.
+ * @brief sort elements in an descending order. Permutation is only correct if all elements are different.
  */
 template<typename ElementType>
 void VectorT<ElementType>::sortDescend(VectorT<int> & permutation) {
-//	VectorT<ElementType> tmp_cp(*this);
-	ippsSortDescend_I(getDataPointer(), this->dataSize);
-//	permutation.resize(this->size());
-//	for (int i = 0; i < this->size(); i++)
-//	{
-//		ElementTyp entry((*this)[i]);
-//		for (int i = 0; i < this->size(); i++)
-//		{
-//			
-//		}
-//	}
+  //copy elements to extract ordering information lateron
+  VectorT<ElementType> tmp_cp(*this);
+
+  // sort the elements
+  ippsSortDescend_I(getDataPointer(), this->dataSize);
+  
+  //compute permutation
+  permutation.resize(this->size());
+
+  int idxSelf ( 0 );
+  for (VectorT<ElementType>::const_iterator itSelf = (*this).begin(); itSelf != (*this).end(); itSelf++, idxSelf++)
+  {
+    int idxOrig ( 0 );
+    for ( VectorT<ElementType>::const_iterator itOrig = tmp_cp.begin(); itOrig != tmp_cp.end(); itOrig++, idxOrig++)
+    {
+      if ( *itOrig == *itSelf )
+      {
+        permutation[idxOrig] = idxSelf;
+        break;
+      }
+    }
+  }
 }
 
 template<typename ElementType>