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

use dedicated function for checking the diagonal

Browse Source
This commit is contained in:
conrad
2022-03-17 13:34:34 +10:00
parent 6e52b7f0d0
commit 05c05cb67f
6 changed files with 46 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/armadillo_bits/op_expmat_meat.hpp
+4 -18
View File
@@ -195,25 +195,9 @@ op_expmat_sym::apply_direct(Mat<typename T1::elem_type>& out, const Base<typenam
arma_debug_check( (X.is_square() == false), "expmat_sym(): given matrix must be square sized" );
const uword N = X.n_rows;
if(is_cx<eT>::yes)
if((arma_config::debug) && (is_cx<eT>::yes) && (sympd_helper::check_diag_imag(X) == false))
{
arma_extra_debug_print("op_expmat_sym: checking imaginary components of diagonal elements");
const T tol = T(10000) * std::numeric_limits<T>::epsilon(); // allow some leeway
const eT* colmem = X.memptr();
for(uword i=0; i<N; ++i)
{
const eT& X_ii = colmem[i];
const T X_ii_imag = access::tmp_imag(X_ii);
if(std::abs(X_ii_imag) > tol) { return false; }
colmem += N;
}
arma_debug_warn_level(1, "inv_sympd(): imaginary components on diagonal are non-zero");
}
if(is_op_diagmat<T1>::value || X.is_diagmat())
@@ -224,6 +208,8 @@ op_expmat_sym::apply_direct(Mat<typename T1::elem_type>& out, const Base<typenam
eT* colmem = out.memptr();
const uword N = X.n_rows;
for(uword i=0; i<N; ++i)
{
eT& out_ii = colmem[i];
include/armadillo_bits/op_inv_spd_meat.hpp
+2 -16
View File
@@ -126,23 +126,9 @@ op_inv_spd_full::apply_direct(Mat<typename T1::elem_type>& out, const Base<typen
// fallthrough if optimisation failed
}
if(is_cx<eT>::yes)
if((arma_config::debug) && (is_cx<eT>::yes) && (sympd_helper::check_diag_imag(out) == false))
{
arma_extra_debug_print("op_inv_spd_full: checking imaginary components of diagonal elements");
const T tol = T(10000) * std::numeric_limits<T>::epsilon(); // allow some leeway
const eT* colmem = out.memptr();
for(uword i=0; i<N; ++i)
{
const eT& out_ii = colmem[i];
const T out_ii_imag = access::tmp_imag(out_ii);
if(std::abs(out_ii_imag) > tol) { return false; }
colmem += N;
}
arma_debug_warn_level(1, "inv_sympd(): imaginary components on diagonal are non-zero");
}
if(is_op_diagmat<T1>::value || out.is_diagmat())
include/armadillo_bits/op_log_det_meat.hpp
+4 -18
View File
@@ -195,25 +195,9 @@ op_log_det_sympd::apply_direct(typename T1::pod_type& out_val, const Base<typena
arma_debug_check( (A.is_square() == false), "log_det_sympd(): given matrix must be square sized" );
const uword N = A.n_rows;
if(is_cx<eT>::yes)
if((arma_config::debug) && (is_cx<eT>::yes) && (sympd_helper::check_diag_imag(A) == false))
{
arma_extra_debug_print("op_log_det_sympd: checking imaginary components of diagonal elements");
const T tol = T(10000) * std::numeric_limits<T>::epsilon(); // allow some leeway
const eT* colmem = A.memptr();
for(uword i=0; i<N; ++i)
{
const eT& A_ii = colmem[i];
const T A_ii_imag = access::tmp_imag(A_ii);
if(std::abs(A_ii_imag) > tol) { return false; }
colmem += N;
}
arma_debug_warn_level(1, "log_det_sympd(): imaginary components on diagonal are non-zero");
}
if(is_op_diagmat<T1>::value || A.is_diagmat())
@@ -224,6 +208,8 @@ op_log_det_sympd::apply_direct(typename T1::pod_type& out_val, const Base<typena
out_val = T(0);
const uword N = A.n_rows;
for(uword i=0; i<N; ++i)
{
const eT& A_ii = colmem[i];
include/armadillo_bits/op_logmat_meat.hpp
+2 -18
View File
@@ -505,25 +505,9 @@ op_logmat_sympd::apply_direct(Mat<typename T1::elem_type>& out, const Base<typen
arma_debug_check( (X.is_square() == false), "logmat_sympd(): given matrix must be square sized" );
if(is_cx<eT>::yes)
if((arma_config::debug) && (is_cx<eT>::yes) && (sympd_helper::check_diag_imag(X) == false))
{
arma_extra_debug_print("op_logmat_sympd: checking imaginary components of diagonal elements");
const T tol = T(10000) * std::numeric_limits<T>::epsilon(); // allow some leeway
const eT* colmem = X.memptr();
const uword N = X.n_rows;
for(uword i=0; i<N; ++i)
{
const eT& X_ii = colmem[i];
const T X_ii_imag = access::tmp_imag(X_ii);
if(std::abs(X_ii_imag) > tol) { return false; }
colmem += N;
}
arma_debug_warn_level(1, "logmat_sympd(): imaginary components on diagonal are non-zero");
}
if(is_op_diagmat<T1>::value || X.is_diagmat())
include/armadillo_bits/op_sqrtmat_meat.hpp
+2 -18
View File
@@ -482,25 +482,9 @@ op_sqrtmat_sympd::apply_direct(Mat<typename T1::elem_type>& out, const Base<type
arma_debug_check( (X.is_square() == false), "sqrtmat_sympd(): given matrix must be square sized" );
if(is_cx<eT>::yes)
if((arma_config::debug) && (is_cx<eT>::yes) && (sympd_helper::check_diag_imag(X) == false))
{
arma_extra_debug_print("op_sqrtmat_sympd: checking imaginary components of diagonal elements");
const T tol = T(10000) * std::numeric_limits<T>::epsilon(); // allow some leeway
const eT* colmem = X.memptr();
const uword N = X.n_rows;
for(uword i=0; i<N; ++i)
{
const eT& X_ii = colmem[i];
const T X_ii_imag = access::tmp_imag(X_ii);
if(std::abs(X_ii_imag) > tol) { return false; }
colmem += N;
}
arma_debug_warn_level(1, "sqrtmat_sympd(): imaginary components on the diagonal are non-zero");
}
if(is_op_diagmat<T1>::value || X.is_diagmat())
include/armadillo_bits/sympd_helper.hpp
+32
View File
@@ -447,6 +447,38 @@ analyse_matrix(bool& is_approx_sym, bool& is_approx_sympd, const Mat<eT>& A)
template<typename eT>
inline
bool
check_diag_imag(const Mat<eT>& A)
{
arma_extra_debug_sigprint();
// NOTE: assuming matrix A is square-sized
typedef typename get_pod_type<eT>::result T;
const T tol = T(10000) * std::numeric_limits<T>::epsilon(); // allow some leeway
const eT* colmem = A.memptr();
const uword N = A.n_rows;
for(uword i=0; i<N; ++i)
{
const eT& A_ii = colmem[i];
const T A_ii_imag = access::tmp_imag(A_ii);
if(std::abs(A_ii_imag) > tol) { return false; }
colmem += N;
}
return true;
}
} // end of namespace sympd_helper