Don't call EigSort automatically in EigSolve, as it could screw up the
ordering of complex eigenvectors

Numeric/EigSolve.cpp

-// $Id: EigSolve.cpp,v 1.1 2004-12-08 05:38:56 geuzaine Exp $
+// $Id: EigSolve.cpp,v 1.2 2004-12-08 16:44:33 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -568,6 +568,10 @@ static void swap_columns(int N, double **v, int i, int k)
 
 void EigSort(int N, double *wr, double *wi, double **B)
 {
+  // Sorts the eigenvalues/vectors in ascending order according to
+  // their real part. Warning: this might screw up the ordering of
+  // complex eigenvectors --- VERIFY THIS
+
   for (int i = 0; i < N - 1; i++){
     int k = i;
     double ek = wr[i];
@@ -593,8 +597,7 @@ int EigSolve(int N, double **A, double *wr, double *wi, double **B,
   // Computes the N eigenvalues of the square, real matrix A. All
   // vectors have dimension N and all matrices have dimension
   // NxN. Warning: the matrix A gets modified during the
-  // computation. The eigenvalues/vectors are sorted in ascending
-  // order according to the real part of the eigenvalues.
+  // computation.
 
   // Balance the matrix to improve accuracy of
   // eigenvalues. (Introduces no rounding errors, since it scales A by
@@ -852,8 +855,6 @@ int EigSolve(int N, double **A, double *wr, double *wi, double **B,
     
   }//End if ierr = -1
 
-  EigSort(N, wr, wi, B);
-  
   return ierr;
 }  //End of Eig3Solve
 

Numeric/Numeric.h

@@ -75,6 +75,7 @@ void FindCubicRoots(const double coeff[4], double re[3], double im[3]);
 void eigsort(double d[3]);
 int EigSolve(int N, double **A, double *wr, double *wi, double **B, 
 	     int *itmp, double *dtmp);
+void EigSort(int N, double *wr, double *wi, double **B);
 
 double InterpolateIso(double *X, double *Y, double *Z, 
 		      double *Val, double V, int I1, int I2, 

Plugin/PrincipalStresses.cpp

-// $Id: PrincipalStresses.cpp,v 1.3 2004-12-08 06:22:23 geuzaine Exp $
+// $Id: PrincipalStresses.cpp,v 1.4 2004-12-08 16:44:33 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -116,6 +116,7 @@ static void principal_stresses(List_T *inList, int inNb, int nbNod, int nbTime,
 	A[1][0] = val[3]; A[1][1] = val[4]; A[1][2] = val[5];
 	A[2][0] = val[6]; A[2][1] = val[7]; A[2][2] = val[8];
 	EigSolve(3, A, wr, wi, B, itmp, dtmp);
+	EigSort(3, wr, wi, B);
 	nbcomplex += nonzero(wi); 
 	//printf("djf=%g %g %g\n", wr[0], wr[1], wr[2]);
 	//printf("vec1=%g %g %g\n", B[0][0], B[1][0], B[2][0]);