Matrices positive definite issues

6 files changed, 21 insertions, 14 deletions

diff --git a/Makefile b/Makefile
index 83a6845..82b1a6b 100644
--- a/Makefile
+++ b/Makefile
@@ -65,7 +65,7 @@ WITH_GSLCBLAS          =                  # Force linking gslcblas (if OpenBlas
 OPENBLAS_LEGACY        =                  # Using older OpenBlas
 FORCE_STATIC           =                  # Static linking of libraries
 # GCC_FLAGS              = -Wall -O3 -std=gnu++11 # extra flags -Wl,--allow-multiple-definition
-GCC_FLAGS              = -Wall -Og -std=gnu++11 # extra flags -Wl,--allow-multiple-definition
+GCC_FLAGS              = -pthread -Wall -std=gnu++11 # extra flags -Wl,--allow-multiple-definition
 TRAVIS_CI              =                  # used by TRAVIS for testing
 
 GSL_INCLUDE_PATH =
@@ -120,10 +120,10 @@ ifdef WITH_OPENBLAS
 endif
 
 ifdef DEBUG
-  CPPFLAGS += -g $(GCC_FLAGS) $(GSL_INCLUDE_PATH) -isystem$(EIGEN_INCLUDE_PATH) -Icontrib/catch-1.9.7 -Isrc $(RPATH)
+  CPPFLAGS += -g -Og $(GCC_FLAGS) $(GSL_INCLUDE_PATH) -isystem$(EIGEN_INCLUDE_PATH) -Icontrib/catch-1.9.7 -Isrc $(RPATH)
 else
   # release mode
-  CPPFLAGS += -DNDEBUG $(GCC_FLAGS) $(GSL_INCLUDE_PATH) -isystem$(EIGEN_INCLUDE_PATH) -Icontrib/catch-1.9.7 -Isrc $(RPATH)
+  CPPFLAGS += -DNDEBUG -O3 $(GCC_FLAGS) $(GSL_INCLUDE_PATH) -isystem$(EIGEN_INCLUDE_PATH) -Icontrib/catch-1.9.7 -Isrc $(RPATH)
 endif
 
 ifeq ($(SYS), WIN)
@@ -135,13 +135,14 @@ ifdef SHOW_COMPILER_WARNINGS
 endif
 
 ifndef FORCE_STATIC
-  LIBS = -lgsl -lopenblas -pthread -lz
+  LIBS = -lgsl -lopenblas -lz
   ifdef WITH_GSLCBLAS
     LIBS += -lgslcblas
   else
     LIBS += -lgfortran -lquadmath
   endif
 else
+  LIBS = -L$(GUIX)/lib -lopenblas -lgsl -lz -lgslcblas -lgfortran -lquadmath
   ifndef TRAVIS_CI # Travis static compile we cheat a little
     CPPFLAGS += -static
   endif
diff --git a/src/debug.cpp b/src/debug.cpp
index 529d603..3a62d2a 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -332,11 +332,11 @@ void do_validate_K(const gsl_matrix *K, const char *__pretty_function, const cha
       warning_at_msg(__file,__line,"K is ill conditioned!");
     if (!isMatrixSymmetric(K))
       warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,"K is not symmetric!" );
-    const bool negative_values = has_negative_values_but_one(eigenvalues);
-    if (negative_values) {
+    const bool negative_eigen_values = has_negative_values_but_one(eigenvalues);
+    if (negative_eigen_values) {
       warning_at_msg(__file,__line,"K has more than one negative eigenvalues!");
     }
-    if (count_small>1 && negative_values && !isMatrixPositiveDefinite(K))
+    if (count_small>1 && negative_eigen_values && !isMatrixPositiveDefinite(K))
       warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,"K is not positive definite!");
     gsl_vector_free(eigenvalues);
   }
diff --git a/src/gemma.cpp b/src/gemma.cpp
index c3caa96..4e3ae88 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1612,6 +1612,7 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) {
     } else if (strcmp(argv[i], "-strict") == 0) {
       // cPar.mode_strict = true;
       debug_set_strict_mode(true);
+      debug_set_debug_mode(true);
     } else if (strcmp(argv[i], "-legacy") == 0) {
       debug_set_legacy_mode(true);
       warning_msg("you are running in legacy mode - support may drop in future versions of gemma");
@@ -3567,5 +3568,7 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
 
   outfile.close();
   outfile.clear();
+
+  info_msg("Done");
   return;
 }
diff --git a/src/lapack.cpp b/src/lapack.cpp
index a66705b..bf73338 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -267,7 +267,7 @@ double EigenDecomp_Zeroed(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
     // checks
     if (gsl_vector_get(eval,i) == 0.0)
       count_zero_eigenvalues += 1;
-    if (gsl_vector_get(eval,i) < 0.0) // count smaller than -EIGEN_MINVALUE
+    if (gsl_vector_get(eval,i) < -EIGEN_MINVALUE) // count smaller than -EIGEN_MINVALUE
       count_negative_eigenvalues += 1;
     d += gsl_vector_get(eval, i);
   }
@@ -279,11 +279,11 @@ double EigenDecomp_Zeroed(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
     msg += " eigenvalues close to zero";
     warning_msg(msg);
   }
-  if (count_negative_eigenvalues > 0) {
+  const bool negative_eigen_values = has_negative_values_but_one(eval);
+  if (negative_eigen_values) {
     write(eval,"eigenvalues");
-    warning_msg("Matrix G has more than one negative eigenvalues!");
+    warning_msg("K has more than one negative eigenvalues!");
   }
-
   return d;
 }
 
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index 542093e..8acbd52 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -104,7 +104,7 @@ bool is_float(const std::string & s){
 }
 
 double safe_log(const double d) {
-  if (!is_legacy_mode() && !is_check_mode())
+  if (!is_legacy_mode() && (is_check_mode() || is_debug_mode()))
     enforce_msg(d > 0.0, (std::string("Trying to take the log of ") + std::to_string(d)).c_str());
   return log(d);
 }
@@ -113,7 +113,7 @@ double safe_sqrt(const double d) {
   double d1 = d;
   if (fabs(d < 0.001))
     d1 = fabs(d);
-  if (!is_legacy_mode() && !is_check_mode())
+  if (!is_legacy_mode() && (is_check_mode() || is_debug_mode()))
     enforce_msg(d1 >= 0.0, (std::string("Trying to take the sqrt of ") + std::to_string(d)).c_str());
   if (d1 < 0.0 )
     return nan("");
@@ -354,7 +354,7 @@ tuple<double, double, double> abs_minmax(const gsl_vector *v) {
 bool has_negative_values_but_one(const gsl_vector *v) {
   bool one_skipped = false;
   for (size_t i=0; i<v->size; i++) {
-    if (v->data[i] < 0.0) {
+    if (v->data[i] < -EIGEN_MINVALUE) {
       if (one_skipped)
         return true;
       one_skipped = true;
diff --git a/test/performance/releases.org b/test/performance/releases.org
index 4497d2f..ca77ad3 100644
--- a/test/performance/releases.org
+++ b/test/performance/releases.org
@@ -10,6 +10,9 @@ core Eigenlib version and 0.97 went multi-core with
 openblas. Unfortunately I linked in lapack and an older BLAS which
 slowed things down. In 0.98 openblas is mostly used and is faster.
 
+Note also that the recent static versions are slower than the
+dynamically linked ones. Not sure why that is.
+
 * GEMMA 0.98-beta1
 
 #+BEGIN_SRC bash