optimised handling of diagonal matrices

This commit is contained in:
conrad
2021-11-18 14:30:04 +10:00
parent ea906f7916
commit 3ee77a6ded
+22 -17
View File
@@ -306,7 +306,7 @@ op_inv_sympd::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT;
// typedef typename T1::pod_type T;
typedef typename T1::pod_type T;
out = expr.get_ref();
@@ -329,22 +329,27 @@ op_inv_sympd::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename
// fallthrough if optimisation failed
}
// if((is_cx<eT>::no) && (out.is_diagmat()))
// {
// const uword N = (std::min)(out.n_rows, out.n_cols);
//
// for(uword i=0; i<N; ++i)
// {
// eT& out_ii = out.at(i,i);
// const T real_out_ii = access::tmp_real(out_ii);
//
// if(real_out_ii <= T(0)) { return false; }
//
// out_ii = eT(T(1) / real_out_ii);
// }
//
// return true;
// }
if((is_cx<eT>::no) && (out.is_diagmat()))
{
// specialised handling of real matrices only;
// currently auxlib::inv_sympd() does not enforce that
// imaginary components of diagonal elements must be zero;
// strictly enforcing this constraint may break existing user software.
const uword N = (std::min)(out.n_rows, out.n_cols);
for(uword i=0; i<N; ++i)
{
eT& out_ii = out.at(i,i);
const T real_out_ii = access::tmp_real(out_ii);
if(real_out_ii <= T(0)) { return false; }
out_ii = eT(T(1) / real_out_ii);
}
return true;
}
return auxlib::inv_sympd(out);
}