diff --git a/include/armadillo_bits/auxlib_bones.hpp b/include/armadillo_bits/auxlib_bones.hpp index f65a8b8c..16fdc935 100644 --- a/include/armadillo_bits/auxlib_bones.hpp +++ b/include/armadillo_bits/auxlib_bones.hpp @@ -44,16 +44,16 @@ class auxlib inline static bool inv_tr_rcond(Mat& A, typename get_pod_type::result& out_rcond, const uword layout); template - inline static bool inv_sympd(Mat& A); + inline static bool inv_sympd(Mat& A, bool& out_sympd_state); template inline static bool inv_sympd(Mat& out, const Mat& X); template - inline static bool inv_sympd_rcond(Mat& A, eT& out_rcond, const eT rcond_threshold); + inline static bool inv_sympd_rcond(Mat& A, bool& out_sympd_state, eT& out_rcond, const eT rcond_threshold); template - inline static bool inv_sympd_rcond(Mat< std::complex >& A, T& out_rcond, const T rcond_threshold); + inline static bool inv_sympd_rcond(Mat< std::complex >& A, bool& out_sympd_state, T& out_rcond, const T rcond_threshold); // diff --git a/include/armadillo_bits/auxlib_meat.hpp b/include/armadillo_bits/auxlib_meat.hpp index 2d1404b5..4efe343d 100644 --- a/include/armadillo_bits/auxlib_meat.hpp +++ b/include/armadillo_bits/auxlib_meat.hpp @@ -245,10 +245,12 @@ auxlib::inv_tr_rcond(Mat& A, typename get_pod_type::result& out_rcond, c template inline bool -auxlib::inv_sympd(Mat& A) +auxlib::inv_sympd(Mat& A, bool& out_sympd_state) { arma_extra_debug_sigprint(); + out_sympd_state = false; + if(A.is_empty()) { return true; } #if defined(ARMA_USE_LAPACK) @@ -266,6 +268,8 @@ auxlib::inv_sympd(Mat& A) if(info != 0) { return false; } + out_sympd_state = true; + arma_extra_debug_print("lapack::potri()"); lapack::potri(&uplo, &n, A.memptr(), &n, &info); @@ -278,6 +282,7 @@ auxlib::inv_sympd(Mat& A) #else { arma_ignore(A); + arma_ignore(out_sympd_state); arma_stop_logic_error("inv_sympd(): use of LAPACK must be enabled"); return false; } @@ -295,7 +300,9 @@ auxlib::inv_sympd(Mat& out, const Mat& X) out = X; - return auxlib::inv_sympd(out); + bool sympd_state_junk = false; + + return auxlib::inv_sympd(out, sympd_state_junk); } @@ -303,10 +310,12 @@ auxlib::inv_sympd(Mat& out, const Mat& X) template inline bool -auxlib::inv_sympd_rcond(Mat& A, eT& out_rcond, const eT rcond_threshold) +auxlib::inv_sympd_rcond(Mat& A, bool& out_sympd_state, eT& out_rcond, const eT rcond_threshold) { arma_extra_debug_sigprint(); + out_sympd_state = false; + if(A.is_empty()) { return true; } #if defined(ARMA_USE_LAPACK) @@ -331,6 +340,8 @@ auxlib::inv_sympd_rcond(Mat& A, eT& out_rcond, const eT rcond_threshold) if(info != 0) { out_rcond = eT(0); return false; } + out_sympd_state = true; + out_rcond = auxlib::lu_rcond_sympd(A, norm_val); if( (rcond_threshold > eT(0)) && (out_rcond < rcond_threshold) ) { return false; } @@ -347,6 +358,7 @@ auxlib::inv_sympd_rcond(Mat& A, eT& out_rcond, const eT rcond_threshold) #else { arma_ignore(A); + arma_ignore(out_sympd_state); arma_ignore(out_rcond); arma_ignore(rcond_threshold); arma_stop_logic_error("inv_sympd_rcond(): use LAPACK must be enabled"); @@ -360,15 +372,18 @@ auxlib::inv_sympd_rcond(Mat& A, eT& out_rcond, const eT rcond_threshold) template inline bool -auxlib::inv_sympd_rcond(Mat< std::complex >& A, T& out_rcond, const T rcond_threshold) +auxlib::inv_sympd_rcond(Mat< std::complex >& A, bool& out_sympd_state, T& out_rcond, const T rcond_threshold) { arma_extra_debug_sigprint(); + out_sympd_state = false; + if(A.is_empty()) { return true; } #if defined(ARMA_CRIPPLED_LAPACK) { arma_ignore(A); + arma_ignore(out_sympd_state); arma_ignore(out_rcond); arma_ignore(rcond_threshold); return false; @@ -393,6 +408,8 @@ auxlib::inv_sympd_rcond(Mat< std::complex >& A, T& out_rcond, const T rcond_t if(info != 0) { out_rcond = T(0); return false; } + out_sympd_state = true; + out_rcond = auxlib::lu_rcond_sympd(A, norm_val); if( (rcond_threshold > T(0)) && (out_rcond < rcond_threshold) ) { return false; } @@ -409,6 +426,7 @@ auxlib::inv_sympd_rcond(Mat< std::complex >& A, T& out_rcond, const T rcond_t #else { arma_ignore(A); + arma_ignore(out_sympd_state); arma_ignore(out_rcond); arma_ignore(rcond_threshold); arma_stop_logic_error("inv_sympd_rcond(): use LAPACK must be enabled"); diff --git a/include/armadillo_bits/op_inv_gen_meat.hpp b/include/armadillo_bits/op_inv_gen_meat.hpp index b33e9c95..fbb08ed5 100644 --- a/include/armadillo_bits/op_inv_gen_meat.hpp +++ b/include/armadillo_bits/op_inv_gen_meat.hpp @@ -166,10 +166,14 @@ op_inv_gen::apply_direct(Mat& out, const Base tmp = out; - const bool status = auxlib::inv_sympd(tmp); + bool sympd_state = false; + + const bool status = auxlib::inv_sympd(tmp, sympd_state); if(status) { out.steal_mem(tmp); return true; } + if((status == false) && (sympd_state == true)) { return false; } + arma_extra_debug_print("op_inv: sympd optimisation failed"); // fallthrough if optimisation failed diff --git a/include/armadillo_bits/op_inv_rcond_meat.hpp b/include/armadillo_bits/op_inv_rcond_meat.hpp index c3e92ee5..1ce87095 100644 --- a/include/armadillo_bits/op_inv_rcond_meat.hpp +++ b/include/armadillo_bits/op_inv_rcond_meat.hpp @@ -97,10 +97,14 @@ op_inv_rcond::apply_direct_gen(Mat& out, typename T1::po Mat tmp = out; - const bool status = auxlib::inv_sympd_rcond(tmp, out_rcond, T(-1)); + bool sympd_state = false; + + const bool status = auxlib::inv_sympd_rcond(tmp, sympd_state, out_rcond, T(-1)); if(status) { out.steal_mem(tmp); return true; } + if((status == false) && (sympd_state == true)) { return false; } + arma_extra_debug_print("op_inv_rcond: sympd optimisation failed"); // fallthrough if optimisation failed diff --git a/include/armadillo_bits/op_inv_spd_meat.hpp b/include/armadillo_bits/op_inv_spd_meat.hpp index 1c76e7f6..abca63e0 100644 --- a/include/armadillo_bits/op_inv_spd_meat.hpp +++ b/include/armadillo_bits/op_inv_spd_meat.hpp @@ -166,7 +166,9 @@ op_inv_spd::apply_direct(Mat& out, const Base& out, const Base::epsilon(); + bool is_sympd_junk = false; + T rcond_junk = T(0); + const T rcond_threshold = T((std::max)(uword(100), uword(A.n_rows))) * std::numeric_limits::epsilon(); - const bool status = auxlib::inv_sympd_rcond(out, rcond_junk, rcond_threshold); + const bool status = auxlib::inv_sympd_rcond(out, is_sympd_junk, rcond_junk, rcond_threshold); if(status) { return true; } arma_extra_debug_print("op_pinv: sympd optimisation failed"); // auxlib::inv_sympd_rcond() will fail if A isn't really positive definite or its rcond is below rcond_threshold } - + ; if(do_sym) { arma_extra_debug_print("op_pinv: symmetric/hermitian optimisation");