diff --git a/Numeric/EigSolve.cpp b/Numeric/EigSolve.cpp index 7ebba60fc8b85910997b93d412ec233cf0e9b53d..4b73c31c5ffcd5e28e17d633510ef22aa9c4128b 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 54953f998b96ee2fd05e97bb00e4846ce45819db..4e59c3bb82efa23b4e0c9ecdb3e29f5e3d2cb469 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 bbde1d57e55d30aa86d3317c3e4e10793becdbf1..66eb982c10aed19958104f5cbe13cd751bf5644b 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]);