phynexis-deps/armadillo-code
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

more robust check for sympd failure

Browse Source
This commit is contained in:
conrad
2022-03-01 22:34:47 +10:00
parent 92509c4be9
commit 45577ff762
6 changed files with 43 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/armadillo_bits/auxlib_bones.hpp
+3 -3
View File
@@ -44,16 +44,16 @@ class auxlib
inline static bool inv_tr_rcond(Mat<eT>& A, typename get_pod_type<eT>::result& out_rcond, const uword layout);
template<typename eT>
inline static bool inv_sympd(Mat<eT>& A);
inline static bool inv_sympd(Mat<eT>& A, bool& out_sympd_state);
template<typename eT>
inline static bool inv_sympd(Mat<eT>& out, const Mat<eT>& X);
template<typename eT>
inline static bool inv_sympd_rcond(Mat<eT>& A, eT& out_rcond, const eT rcond_threshold);
inline static bool inv_sympd_rcond(Mat<eT>& A, bool& out_sympd_state, eT& out_rcond, const eT rcond_threshold);
template<typename T>
inline static bool inv_sympd_rcond(Mat< std::complex<T> >& A, T& out_rcond, const T rcond_threshold);
inline static bool inv_sympd_rcond(Mat< std::complex<T> >& A, bool& out_sympd_state, T& out_rcond, const T rcond_threshold);
//
include/armadillo_bits/auxlib_meat.hpp
+22 -4
View File
@@ -245,10 +245,12 @@ auxlib::inv_tr_rcond(Mat<eT>& A, typename get_pod_type<eT>::result& out_rcond, c
template<typename eT>
inline
bool
auxlib::inv_sympd(Mat<eT>& A)
auxlib::inv_sympd(Mat<eT>& 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<eT>& 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<eT>& 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<eT>& out, const Mat<eT>& 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<eT>& out, const Mat<eT>& X)
template<typename eT>
inline
bool
auxlib::inv_sympd_rcond(Mat<eT>& A, eT& out_rcond, const eT rcond_threshold)
auxlib::inv_sympd_rcond(Mat<eT>& 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<eT>& 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<T>(A, norm_val);
if( (rcond_threshold > eT(0)) && (out_rcond < rcond_threshold) ) { return false; }
@@ -347,6 +358,7 @@ auxlib::inv_sympd_rcond(Mat<eT>& 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<eT>& A, eT& out_rcond, const eT rcond_threshold)
template<typename T>
inline
bool
auxlib::inv_sympd_rcond(Mat< std::complex<T> >& A, T& out_rcond, const T rcond_threshold)
auxlib::inv_sympd_rcond(Mat< std::complex<T> >& 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<T> >& 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<T>(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<T> >& 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");
include/armadillo_bits/op_inv_gen_meat.hpp
+5 -1
View File
@@ -166,10 +166,14 @@ op_inv_gen::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T
Mat<eT> 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
include/armadillo_bits/op_inv_rcond_meat.hpp
+5 -1
View File
@@ -97,10 +97,14 @@ op_inv_rcond::apply_direct_gen(Mat<typename T1::elem_type>& out, typename T1::po
Mat<eT> 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
include/armadillo_bits/op_inv_spd_meat.hpp
+3 -1
View File
@@ -166,7 +166,9 @@ op_inv_spd::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T
return true;
}
return auxlib::inv_sympd(out);
bool sympd_state_junk = false;
return auxlib::inv_sympd(out, sympd_state_junk);
}
include/armadillo_bits/op_pinv_meat.hpp
+5 -4
View File
@@ -103,17 +103,18 @@ op_pinv::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1::
out = A;
T rcond_junk;
const T rcond_threshold = T((std::max)(uword(100), uword(A.n_rows))) * std::numeric_limits<T>::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<T>::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");