Fixes regression #179

2 files changed, 15 insertions, 4 deletions

diff --git a/src/lapack.cpp b/src/lapack.cpp
index eb5b16b..31b6add 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -253,8 +253,10 @@ double EigenDecomp(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
   return d;
 }
 
-// Same as EigenDecomp but zeroes eigenvalues close to zero. When
-// negative eigenvalues remain a warning is issued.
+// Does NOT set eigenvalues to be positive. G gets destroyed. Returns
+// eigen trace and values in U and eval (eigenvalues).  Same as
+// EigenDecomp but zeroes eigenvalues close to zero. When negative
+// eigenvalues remain a warning is issued.
 double EigenDecomp_Zeroed(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
                           const size_t flag_largematrix) {
   EigenDecomp(G,U,eval,flag_largematrix);
@@ -262,7 +264,8 @@ double EigenDecomp_Zeroed(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
   int count_zero_eigenvalues = 0;
   int count_negative_eigenvalues = 0;
   for (size_t i = 0; i < eval->size; i++) {
-    if (std::abs(gsl_vector_get(eval, i)) < EIGEN_MINVALUE)
+    // if (std::abs(gsl_vector_get(eval, i)) < EIGEN_MINVALUE)
+    if (gsl_vector_get(eval, i) < 1e-10)
       gsl_vector_set(eval, i, 0.0);
     // checks
     if (gsl_vector_get(eval,i) == 0.0)
diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp
index a3c7c85..4b8db05 100644
--- a/src/mvlmm.cpp
+++ b/src/mvlmm.cpp
@@ -256,7 +256,15 @@ double EigenProc(const gsl_matrix *V_g, const gsl_matrix *V_e, gsl_vector *D_l,
   gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, V_e_hi, VgVehi, 0.0, Lambda);
 
   // Eigen decomposition of Lambda.
-  EigenDecomp_Zeroed(Lambda, U_l, D_l, 0);
+  EigenDecomp(Lambda, U_l, D_l, 0);
+
+  for (size_t i = 0; i < d_size; i++) {
+    d = gsl_vector_get(D_l, i);
+    if (d < 0) {
+      gsl_vector_set(D_l, i, 0);
+    }
+  }
+
 
   // Calculate UltVeh and UltVehi.
   gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, U_l, V_e_h, 0.0, UltVeh);