From 2992635aacf9ee1e270baab124eb2a7c6c019a6b Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 8 Dec 2004 16:44:33 +0000
Subject: [PATCH] Don't call EigSort automatically in EigSolve, as it could
screw up the ordering of complex eigenvectors
---
Numeric/EigSolve.cpp | 11 ++++++-----
Numeric/Numeric.h | 1 +
Plugin/PrincipalStresses.cpp | 3 ++-
3 files changed, 9 insertions(+), 6 deletions(-)
diff --git a/Numeric/EigSolve.cpp b/Numeric/EigSolve.cpp
index 7ebba60fc8..4b73c31c5f 100644
--- a/Numeric/EigSolve.cpp
+++ b/Numeric/EigSolve.cpp
@@ -1,4 +1,4 @@
-// $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
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 54953f998b..4e59c3bb82 100644
--- a/Numeric/Numeric.h
+++ b/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,
diff --git a/Plugin/PrincipalStresses.cpp b/Plugin/PrincipalStresses.cpp
index bbde1d57e5..66eb982c10 100644
--- a/Plugin/PrincipalStresses.cpp
+++ b/Plugin/PrincipalStresses.cpp
@@ -1,4 +1,4 @@
-// $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]);
--
GitLab