diff --git a/include/armadillo_bits/sym_helper.hpp b/include/armadillo_bits/sym_helper.hpp index 221eb3ff..2134885d 100644 --- a/include/armadillo_bits/sym_helper.hpp +++ b/include/armadillo_bits/sym_helper.hpp @@ -50,19 +50,26 @@ guess_sympd_worker(const Mat& A) const eT* A_mem = A.memptr(); const eT* A_col = A_mem; + bool diag_below_tol = true; + eT max_diag = eT(0); for(uword j=0; j < N; ++j) { const eT A_jj = A_col[j]; - if(A_jj <= eT(0)) { return false; } + if( A_jj <= eT(0)) { return false; } + if(arma_isfinite(A_jj) == false) { return false; } + + if(A_jj >= tol) { diag_below_tol = false; } max_diag = (A_jj > max_diag) ? A_jj : max_diag; A_col += N; } + if(diag_below_tol) { return false; } // assume matrix is suspect if all diagonal elements are close to zero + A_col = A_mem; const uword Nm1 = N-1; @@ -128,6 +135,8 @@ guess_sympd_worker(const Mat& A) const eT* A_mem = A.memptr(); const eT* A_col = A_mem; + bool diag_below_tol = true; + T max_diag = T(0); for(uword j=0; j < N; ++j) @@ -138,15 +147,21 @@ guess_sympd_worker(const Mat& A) const T A_jj_rabs = std::abs(A_jj_r); const T A_jj_iabs = std::abs(A_jj_i); - if(A_jj_r <= T(0) ) { return false; } // real should be positive - if(A_jj_iabs > tol ) { return false; } // imag should be approx zero - if(A_jj_iabs > A_jj_rabs) { return false; } // corner case: real and imag are close to zero, and imag is dominant + if( A_jj_r <= T(0) ) { return false; } // real should be positive + if(arma_isfinite(A_jj_r) == false) { return false; } + + if(A_jj_iabs > tol ) { return false; } // imag should be approx zero + if(A_jj_iabs > A_jj_rabs) { return false; } // corner case: real and imag are close to zero, and imag is dominant + + if(A_jj_r >= tol) { diag_below_tol = false; } max_diag = (A_jj_r > max_diag) ? A_jj_r : max_diag; A_col += N; } + if(diag_below_tol) { return false; } // assume matrix is suspect if all diagonal elements are close to zero + const T square_max_diag = max_diag * max_diag; if(arma_isfinite(square_max_diag) == false) { return false; } @@ -268,7 +283,7 @@ is_approx_sym_worker(const Mat& A) for(uword j=0; j < N; ++j) { - const eT& A_jj = A_col[j]; + const eT A_jj = A_col[j]; if(arma_isfinite(A_jj) == false) { return false; }