From 65c42d01cdf11a2d2bf30d308cd5cc5c4de8ce1d Mon Sep 17 00:00:00 2001 From: conrad Date: Sun, 26 Oct 2025 20:07:38 +1000 Subject: [PATCH] faster handling of submatrix rows --- include/armadillo_bits/op_accu_meat.hpp | 33 ++- include/armadillo_bits/subview_bones.hpp | 18 +- include/armadillo_bits/subview_meat.hpp | 258 +++++++++++++++++++---- 3 files changed, 265 insertions(+), 44 deletions(-) diff --git a/include/armadillo_bits/op_accu_meat.hpp b/include/armadillo_bits/op_accu_meat.hpp index b958c31d..3cd71064 100644 --- a/include/armadillo_bits/op_accu_meat.hpp +++ b/include/armadillo_bits/op_accu_meat.hpp @@ -602,7 +602,30 @@ op_accu_mat::apply(const subview& X) const uword X_n_rows = X.n_rows; const uword X_n_cols = X.n_cols; - if(X_n_rows == 1) { return op_accu_mat::apply( static_cast< const subview_row& >(X) ); } + if(X_n_rows == 1) + { + const uword X_m_n_rows = X.m.n_rows; + + const eT* mem_ptr = X.colptr(0); + + eT val1 = eT(0); + eT val2 = eT(0); + + uword j; + + for(j=1; j < X_n_cols; j+=2) + { + val1 += (*mem_ptr); mem_ptr += X_m_n_rows; + val2 += (*mem_ptr); mem_ptr += X_m_n_rows; + } + + if((j-1) < X_n_cols) + { + val1 += (*mem_ptr); + } + + return val1 + val2; + } if(X_n_cols == 1) { return arrayops::accumulate( X.colptr(0), X_n_rows ); } @@ -640,7 +663,7 @@ op_accu_mat::apply(const subview_row& X) const uword X_m_n_rows = X.m.n_rows; const uword X_n_cols = X.n_cols; - const eT* row_mem = &(X.m.at(X.aux_row1,X.aux_col1)); + const eT* mem_ptr = X.rowmem; eT val1 = eT(0); eT val2 = eT(0); @@ -649,13 +672,13 @@ op_accu_mat::apply(const subview_row& X) for(j=1; j < X_n_cols; j+=2) { - val1 += (*row_mem); row_mem += X_m_n_rows; - val2 += (*row_mem); row_mem += X_m_n_rows; + val1 += (*mem_ptr); mem_ptr += X_m_n_rows; + val2 += (*mem_ptr); mem_ptr += X_m_n_rows; } if((j-1) < X_n_cols) { - val1 += (*row_mem); + val1 += (*mem_ptr); } return val1 + val2; diff --git a/include/armadillo_bits/subview_bones.hpp b/include/armadillo_bits/subview_bones.hpp index 06cf9a56..f3837606 100644 --- a/include/armadillo_bits/subview_bones.hpp +++ b/include/armadillo_bits/subview_bones.hpp @@ -85,7 +85,7 @@ class subview : public Base< eT, subview > template inline void operator-= (const SpBase& x); template inline void operator%= (const SpBase& x); template inline void operator/= (const SpBase& x); - + template inline typename enable_if2< is_same_type::value, void>::result operator=(const Gen& x); @@ -396,6 +396,11 @@ class subview_col : public subview inline void zeros(); inline void ones(); + arma_warn_unused inline bool is_finite() const; + + arma_warn_unused inline bool has_inf() const; + arma_warn_unused inline bool has_nan() const; + arma_inline eT at_alt (const uword i) const; arma_inline eT& operator[](const uword i); @@ -528,6 +533,8 @@ class subview_row : public subview static constexpr bool is_col = false; static constexpr bool is_xvec = false; + const eT* rowmem; + inline void operator= (const subview& x); inline void operator= (const subview_row& x); inline void operator= (const eT val); @@ -545,6 +552,15 @@ class subview_row : public subview arma_warn_unused arma_inline const Op,op_strans> as_col() const; + inline void fill(const eT val); + inline void zeros(); + inline void ones(); + + arma_warn_unused inline bool is_finite() const; + + arma_warn_unused inline bool has_inf() const; + arma_warn_unused inline bool has_nan() const; + inline eT at_alt (const uword i) const; inline eT& operator[](const uword i); diff --git a/include/armadillo_bits/subview_meat.hpp b/include/armadillo_bits/subview_meat.hpp index 907ed8c5..1146b7bd 100644 --- a/include/armadillo_bits/subview_meat.hpp +++ b/include/armadillo_bits/subview_meat.hpp @@ -1079,16 +1079,9 @@ subview::fill(const eT val) eT* Aptr = &(A.at(s.aux_row1,s.aux_col1)); - uword jj; - for(jj=1; jj < s_n_cols; jj+=2) + for(uword ii=0; ii < s_n_cols; ++ii) { (*Aptr) = val; Aptr += A_n_rows; - (*Aptr) = val; Aptr += A_n_rows; - } - - if((jj-1) < s_n_cols) - { - (*Aptr) = val; } } else @@ -3380,7 +3373,7 @@ subview_col::operator=(const Base& expr) if(is_Mat::value) { - const unwrap U(expr.get_ref()); + const unwrap U(expr.get_ref()); // deliberately not using quasi_unwrap arma_conform_assert_same_size(subview::n_rows, uword(1), U.M.n_rows, U.M.n_cols, "copy into submatrix"); @@ -3498,6 +3491,48 @@ subview_col::ones() +template +inline +bool +subview_col::is_finite() const + { + arma_debug_sigprint(); + + if(arma_config::fast_math_warn) { arma_warn(1, "is_finite(): detection of non-finite values is not reliable in fast math mode"); } + + return arrayops::is_finite(colmem, subview::n_rows); + } + + + +template +inline +bool +subview_col::has_inf() const + { + arma_debug_sigprint(); + + if(arma_config::fast_math_warn) { arma_warn(1, "has_inf(): detection of non-finite values is not reliable in fast math mode"); } + + return arrayops::has_inf(colmem, subview::n_rows); + } + + + +template +inline +bool +subview_col::has_nan() const + { + arma_debug_sigprint(); + + if(arma_config::fast_math_warn) { arma_warn(1, "has_nan(): detection of non-finite values is not reliable in fast math mode"); } + + return arrayops::has_nan(colmem, subview::n_rows); + } + + + template arma_inline eT @@ -4233,6 +4268,7 @@ template inline subview_row::subview_row(const Mat& in_m, const uword in_row) : subview(in_m, in_row, 0, 1, in_m.n_cols) + , rowmem(subview::colptr(0)) { arma_debug_sigprint(); } @@ -4243,6 +4279,7 @@ template inline subview_row::subview_row(const Mat& in_m, const uword in_row, const uword in_col1, const uword in_n_cols) : subview(in_m, in_row, in_col1, 1, in_n_cols) + , rowmem(subview::colptr(0)) { arma_debug_sigprint(); } @@ -4253,6 +4290,7 @@ template inline subview_row::subview_row(const subview_row& in) : subview(in) // interprets 'subview_row' as 'subview' + , rowmem(in.rowmem) { arma_debug_sigprint(); } @@ -4263,8 +4301,11 @@ template inline subview_row::subview_row(subview_row&& in) : subview(std::move(in)) // interprets 'subview_row' as 'subview' + , rowmem(in.rowmem) { arma_debug_sigprint(); + + access::rw(in.rowmem) = nullptr; } @@ -4300,7 +4341,12 @@ subview_row::operator=(const eT val) { arma_debug_sigprint(); - subview::operator=(val); // interprets 'subview_row' as 'subview' + if(subview::n_elem != 1) + { + arma_conform_assert_same_size(subview::n_rows, subview::n_cols, 1, 1, "copy into submatrix"); + } + + access::rw( rowmem[0] ) = val; } @@ -4335,7 +4381,39 @@ subview_row::operator=(const Base& X) { arma_debug_sigprint(); - subview::operator=(X); + if(is_Mat::value) + { + const unwrap U(X.get_ref()); // deliberately not using quasi_unwrap + + arma_conform_assert_same_size(uword(1), subview::n_cols, U.M.n_rows, U.M.n_cols, "copy into submatrix"); + + const eT* UM_mem = U.M.memptr(); + + eT* mem_ptr = access::rwp(rowmem); + + const uword local_s_n_cols = subview::n_cols; + const uword local_m_n_rows = subview::m.n_rows; + + uword j; + + for(j=1; j < local_s_n_cols; j+=2) + { + const eT val_i = (*UM_mem); UM_mem++; + const eT val_j = (*UM_mem); UM_mem++; + + (*mem_ptr) = val_i; mem_ptr += local_m_n_rows; + (*mem_ptr) = val_j; mem_ptr += local_m_n_rows; + } + + if((j-1) < local_s_n_cols) + { + (*mem_ptr) = (*UM_mem); + } + } + else + { + subview::operator=(X); + } } @@ -4408,14 +4486,134 @@ subview_row::as_col() const +template +inline +void +subview_row::fill(const eT val) + { + arma_debug_sigprint(); + + eT* mem_ptr = access::rwp(rowmem); + + const uword local_s_n_cols = subview::n_cols; + const uword local_m_n_rows = subview::m.n_rows; + + for(uword ii=0; ii < local_s_n_cols; ++ii) + { + (*mem_ptr) = val; mem_ptr += local_m_n_rows; + } + } + + + +template +inline +void +subview_row::zeros() + { + arma_debug_sigprint(); + + (*this).fill(eT(0)); + } + + + +template +inline +void +subview_row::ones() + { + arma_debug_sigprint(); + + (*this).fill(eT(1)); + } + + + +template +inline +bool +subview_row::is_finite() const + { + arma_debug_sigprint(); + + if(arma_config::fast_math_warn) { arma_warn(1, "is_finite(): detection of non-finite values is not reliable in fast math mode"); } + + const eT* mem_ptr = rowmem; + + const uword local_s_n_cols = subview::n_cols; + const uword local_m_n_rows = subview::m.n_rows; + + for(uword ii=0; ii < local_s_n_cols; ++ii) + { + const eT val = (*mem_ptr); mem_ptr += local_m_n_rows; + + if(arma_isnonfinite(val)) { return false; } + } + + return true; + } + + + +template +inline +bool +subview_row::has_inf() const + { + arma_debug_sigprint(); + + if(arma_config::fast_math_warn) { arma_warn(1, "has_inf(): detection of non-finite values is not reliable in fast math mode"); } + + const eT* mem_ptr = rowmem; + + const uword local_s_n_cols = subview::n_cols; + const uword local_m_n_rows = subview::m.n_rows; + + for(uword ii=0; ii < local_s_n_cols; ++ii) + { + const eT val = (*mem_ptr); mem_ptr += local_m_n_rows; + + if(arma_isinf(val)) { return true; } + } + + return false; + } + + + +template +inline +bool +subview_row::has_nan() const + { + arma_debug_sigprint(); + + if(arma_config::fast_math_warn) { arma_warn(1, "has_nan(): detection of non-finite values is not reliable in fast math mode"); } + + const eT* mem_ptr = rowmem; + + const uword local_s_n_cols = subview::n_cols; + const uword local_m_n_rows = subview::m.n_rows; + + for(uword ii=0; ii < local_s_n_cols; ++ii) + { + const eT val = (*mem_ptr); mem_ptr += local_m_n_rows; + + if(arma_isnan(val)) { return true; } + } + + return false; + } + + + template inline eT subview_row::at_alt(const uword ii) const { - const uword index = (ii + (subview::aux_col1))*(subview::m).n_rows + (subview::aux_row1); - - return subview::m.mem[index]; + return rowmem[ii * subview::m.n_rows]; } @@ -4425,9 +4623,7 @@ inline eT& subview_row::operator[](const uword ii) { - const uword index = (ii + (subview::aux_col1))*(subview::m).n_rows + (subview::aux_row1); - - return access::rw( (const_cast< Mat& >(subview::m)).mem[index] ); + return access::rw( rowmem[ii * subview::m.n_rows] ); } @@ -4437,9 +4633,7 @@ inline eT subview_row::operator[](const uword ii) const { - const uword index = (ii + (subview::aux_col1))*(subview::m).n_rows + (subview::aux_row1); - - return subview::m.mem[index]; + return rowmem[ii * subview::m.n_rows]; } @@ -4450,10 +4644,8 @@ eT& subview_row::operator()(const uword ii) { arma_conform_check_bounds( (ii >= subview::n_elem), "subview::operator(): index out of bounds" ); - - const uword index = (ii + (subview::aux_col1))*(subview::m).n_rows + (subview::aux_row1); - return access::rw( (const_cast< Mat& >(subview::m)).mem[index] ); + return access::rw( rowmem[ii * subview::m.n_rows] ); } @@ -4465,9 +4657,7 @@ subview_row::operator()(const uword ii) const { arma_conform_check_bounds( (ii >= subview::n_elem), "subview::operator(): index out of bounds" ); - const uword index = (ii + (subview::aux_col1))*(subview::m).n_rows + (subview::aux_row1); - - return subview::m.mem[index]; + return rowmem[ii * subview::m.n_rows]; } @@ -4479,9 +4669,7 @@ subview_row::operator()(const uword in_row, const uword in_col) { arma_conform_check_bounds( ((in_row > 0) || (in_col >= subview::n_cols)), "subview::operator(): index out of bounds" ); - const uword index = (in_col + (subview::aux_col1))*(subview::m).n_rows + (subview::aux_row1); - - return access::rw( (const_cast< Mat& >(subview::m)).mem[index] ); + return access::rw( rowmem[in_col * subview::m.n_rows] ); } @@ -4493,9 +4681,7 @@ subview_row::operator()(const uword in_row, const uword in_col) const { arma_conform_check_bounds( ((in_row > 0) || (in_col >= subview::n_cols)), "subview::operator(): index out of bounds" ); - const uword index = (in_col + (subview::aux_col1))*(subview::m).n_rows + (subview::aux_row1); - - return subview::m.mem[index]; + return rowmem[in_col * subview::m.n_rows]; } @@ -4505,9 +4691,7 @@ inline eT& subview_row::at(const uword, const uword in_col) { - const uword index = (in_col + (subview::aux_col1))*(subview::m).n_rows + (subview::aux_row1); - - return access::rw( (const_cast< Mat& >(subview::m)).mem[index] ); + return access::rw( rowmem[in_col * subview::m.n_rows] ); } @@ -4517,9 +4701,7 @@ inline eT subview_row::at(const uword, const uword in_col) const { - const uword index = (in_col + (subview::aux_col1))*(subview::m).n_rows + (subview::aux_row1); - - return subview::m.mem[index]; + return rowmem[in_col * subview::m.n_rows]; }