From 5eabafe2b3ece8fb39ed31b8f635dd2ebd424c53 Mon Sep 17 00:00:00 2001 From: conrad Date: Wed, 2 Oct 2024 14:46:00 +1000 Subject: [PATCH] partial rewrite of powext --- include/armadillo_bits/arma_forward.hpp | 12 + include/armadillo_bits/fn_powext.hpp | 20 + include/armadillo_bits/glue_powext_bones.hpp | 13 +- include/armadillo_bits/glue_powext_meat.hpp | 414 +++---------------- 4 files changed, 93 insertions(+), 366 deletions(-) diff --git a/include/armadillo_bits/arma_forward.hpp b/include/armadillo_bits/arma_forward.hpp index b35d64ef..3155f390 100644 --- a/include/armadillo_bits/arma_forward.hpp +++ b/include/armadillo_bits/arma_forward.hpp @@ -243,6 +243,18 @@ struct traits_glue_or }; +struct traits_glue_first_arg + { + template + struct traits + { + static constexpr bool is_row = T1::is_row; + static constexpr bool is_col = T1::is_col; + static constexpr bool is_xvec = T1::is_xvec; + }; + }; + + template class gemm; template class gemv; diff --git a/include/armadillo_bits/fn_powext.hpp b/include/armadillo_bits/fn_powext.hpp index e13a1304..b49c09f7 100644 --- a/include/armadillo_bits/fn_powext.hpp +++ b/include/armadillo_bits/fn_powext.hpp @@ -82,6 +82,26 @@ pow +// TODO +// template +// arma_warn_unused +// inline +// Cube +// pow +// ( +// const BaseCube& X, +// const Base & Y +// ) +// { +// arma_debug_sigprint(); +// +// typedef typename T1::elem_type eT; +// +// ... +// } + + + template arma_warn_unused inline diff --git a/include/armadillo_bits/glue_powext_bones.hpp b/include/armadillo_bits/glue_powext_bones.hpp index d5698c5c..3a6f1a1a 100644 --- a/include/armadillo_bits/glue_powext_bones.hpp +++ b/include/armadillo_bits/glue_powext_bones.hpp @@ -1,4 +1,3 @@ - // SPDX-License-Identifier: Apache-2.0 // // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) @@ -24,7 +23,7 @@ class glue_powext - : public traits_glue_or + : public traits_glue_first_arg { public: @@ -32,21 +31,17 @@ class glue_powext template inline static void apply(Mat& out, const Mat& A, const Mat& B); - template inline static Mat apply(const subview_each1& X, const Base& Y); - // template inline static void apply(Cube& out, const GlueCube& X); template inline static void apply(Cube& out, const Cube& A, const Cube& B); - - template inline static Cube apply(const subview_cube_each1& X, const Base& Y); }; class glue_powext_cx - : public traits_glue_or + : public traits_glue_first_arg { public: @@ -54,15 +49,11 @@ class glue_powext_cx template inline static void apply(Mat< std::complex >& out, const Mat< std::complex >& A, const Mat& B); - template inline static Mat apply(const subview_each1& X, const Base& Y); - // template inline static void apply(Cube& out, const mtGlueCube& X); template inline static void apply(Cube< std::complex >& out, const Cube< std::complex >& A, const Cube& B); - - template inline static Cube< std::complex > apply(const subview_cube_each1< std::complex >& X, const Base& Y); }; diff --git a/include/armadillo_bits/glue_powext_meat.hpp b/include/armadillo_bits/glue_powext_meat.hpp index 927ead6c..0d948c41 100644 --- a/include/armadillo_bits/glue_powext_meat.hpp +++ b/include/armadillo_bits/glue_powext_meat.hpp @@ -36,6 +36,49 @@ glue_powext::apply(Mat& out, const Glue& A = UA.M; const Mat& B = UB.M; + if( (A.is_vec() == false) && ((T2::is_row || B.is_rowvec()) || (T2::is_col || B.is_colvec())) ) + { + // rudimentary handling of broadcasting operations + // mainly for compat with previous ill-designed direct handling of .each_row() and .each_col() + + Mat BB; + + if(T2::is_row || B.is_rowvec()) + { + arma_conform_assert_same_size(A.n_rows, A.n_cols, A.n_rows, B.n_cols, "element-wise pow()"); + + BB.set_size(A.n_rows, B.n_cols); + + BB.each_row() = B; + } + else + if(T2::is_col || B.is_colvec()) + { + arma_conform_assert_same_size(A.n_rows, A.n_cols, B.n_rows, A.n_cols, "element-wise pow()"); + + BB.set_size(B.n_rows, A.n_cols); + + BB.each_col() = B; + } + + const bool UA_bad_alias = UA.is_alias(out) && (UA.has_subview); // allow inplace operation + + if(UA_bad_alias) + { + Mat tmp; + + glue_powext::apply(tmp, A, BB); + + out.steal_mem(tmp); + } + else + { + glue_powext::apply(out, A, BB); + } + + return; + } + arma_conform_assert_same_size(A, B, "element-wise pow()"); const bool UA_bad_alias = UA.is_alias(out) && (UA.has_subview); // allow inplace operation @@ -97,116 +140,6 @@ glue_powext::apply(Mat& out, const Mat& A, const Mat& B) -template -inline -Mat -glue_powext::apply - ( - const subview_each1& X, - const Base& Y - ) - { - arma_debug_sigprint(); - - typedef typename parent::elem_type eT; - - const parent& A = X.P; - - const uword A_n_rows = A.n_rows; - const uword A_n_cols = A.n_cols; - - Mat out(A_n_rows, A_n_cols, arma_nozeros_indicator()); - - const quasi_unwrap tmp(Y.get_ref()); - const Mat& B = tmp.M; - - X.check_size(B); - - const eT* B_mem = B.memptr(); - - if(mode == 0) // each column - { - if( arma_config::openmp && mp_gate::eval(A.n_elem) ) - { - #if defined(ARMA_USE_OPENMP) - { - const int n_threads = int( (std::min)(uword(mp_thread_limit::get()), A_n_cols) ); - - #pragma omp parallel for schedule(static) num_threads(n_threads) - for(uword i=0; i < A_n_cols; ++i) - { - const eT* A_mem = A.colptr(i); - eT* out_mem = out.colptr(i); - - for(uword row=0; row < A_n_rows; ++row) - { - out_mem[row] = eop_aux::pow(A_mem[row], B_mem[row]); - } - } - } - #endif - } - else - { - for(uword i=0; i < A_n_cols; ++i) - { - const eT* A_mem = A.colptr(i); - eT* out_mem = out.colptr(i); - - for(uword row=0; row < A_n_rows; ++row) - { - out_mem[row] = eop_aux::pow(A_mem[row], B_mem[row]); - } - } - } - } - - if(mode == 1) // each row - { - if( arma_config::openmp && mp_gate::eval(A.n_elem) ) - { - #if defined(ARMA_USE_OPENMP) - { - const int n_threads = int( (std::min)(uword(mp_thread_limit::get()), A_n_cols) ); - - #pragma omp parallel for schedule(static) num_threads(n_threads) - for(uword i=0; i < A_n_cols; ++i) - { - const eT* A_mem = A.colptr(i); - eT* out_mem = out.colptr(i); - - const eT B_val = B_mem[i]; - - for(uword row=0; row < A_n_rows; ++row) - { - out_mem[row] = eop_aux::pow(A_mem[row], B_val); - } - } - } - #endif - } - else - { - for(uword i=0; i < A_n_cols; ++i) - { - const eT* A_mem = A.colptr(i); - eT* out_mem = out.colptr(i); - - const eT B_val = B_mem[i]; - - for(uword row=0; row < A_n_rows; ++row) - { - out_mem[row] = eop_aux::pow(A_mem[row], B_val); - } - } - } - } - - return out; - } - - - template inline void @@ -222,6 +155,22 @@ glue_powext::apply(Cube& out, const GlueCube& A = UA.M; const Cube& B = UB.M; + if((A.n_slices != 1) && (B.n_slices == 1)) + { + // rudimentary handling of broadcasting operations + // mainly for compat with previous ill-designed direct handling of .each_slice() + + arma_conform_assert_same_size(A.n_rows, A.n_cols, A.n_slices, B.n_rows, B.n_cols, A.n_slices, "element-wise pow()"); + + Cube BB(B.n_rows, B.n_cols, A.n_slices, arma_nozeros_indicator()); + + BB.each_slice() = B.slice(0); + + glue_powext::apply(out, A, BB); + + return; + } + arma_conform_assert_same_size(A, B, "element-wise pow()"); if(UB.is_alias(out)) @@ -280,72 +229,6 @@ glue_powext::apply(Cube& out, const Cube& A, const Cube& B) -template -inline -Cube -glue_powext::apply - ( - const subview_cube_each1& X, - const Base& Y - ) - { - arma_debug_sigprint(); - - const Cube& A = X.P; - - const uword A_n_rows = A.n_rows; - const uword A_n_cols = A.n_cols; - const uword A_n_slices = A.n_slices; - - Cube out(A_n_rows, A_n_cols, A_n_slices, arma_nozeros_indicator()); - - const quasi_unwrap tmp(Y.get_ref()); - const Mat& B = tmp.M; - - X.check_size(B); - - const eT* B_mem = B.memptr(); - const uword B_n_elem = B.n_elem; - - if( arma_config::openmp && mp_gate::eval(A.n_elem) ) - { - #if defined(ARMA_USE_OPENMP) - { - const int n_threads = int( (std::min)(uword(mp_thread_limit::get()), A_n_slices) ); - - #pragma omp parallel for schedule(static) num_threads(n_threads) - for(uword s=0; s < A_n_slices; ++s) - { - const eT* A_slice_mem = A.slice_memptr(s); - eT* out_slice_mem = out.slice_memptr(s); - - for(uword i=0; i < B_n_elem; ++i) - { - out_slice_mem[i] = eop_aux::pow(A_slice_mem[i], B_mem[i]); - } - } - } - #endif - } - else - { - for(uword s=0; s < A_n_slices; ++s) - { - const eT* A_slice_mem = A.slice_memptr(s); - eT* out_slice_mem = out.slice_memptr(s); - - for(uword i=0; i < B_n_elem; ++i) - { - out_slice_mem[i] = eop_aux::pow(A_slice_mem[i], B_mem[i]); - } - } - } - - return out; - } - - - // @@ -426,117 +309,6 @@ glue_powext_cx::apply(Mat< std::complex >& out, const Mat< std::complex >& -template -inline -Mat -glue_powext_cx::apply - ( - const subview_each1& X, - const Base& Y - ) - { - arma_debug_sigprint(); - - typedef typename parent::elem_type eT; - typedef typename parent::pod_type T; - - const parent& A = X.P; - - const uword A_n_rows = A.n_rows; - const uword A_n_cols = A.n_cols; - - Mat out(A_n_rows, A_n_cols, arma_nozeros_indicator()); - - const quasi_unwrap tmp(Y.get_ref()); - const Mat& B = tmp.M; - - X.check_size(B); - - const T* B_mem = B.memptr(); - - if(mode == 0) // each column - { - if( arma_config::openmp && mp_gate::eval(A.n_elem) ) - { - #if defined(ARMA_USE_OPENMP) - { - const int n_threads = int( (std::min)(uword(mp_thread_limit::get()), A_n_cols) ); - - #pragma omp parallel for schedule(static) num_threads(n_threads) - for(uword i=0; i < A_n_cols; ++i) - { - const eT* A_mem = A.colptr(i); - eT* out_mem = out.colptr(i); - - for(uword row=0; row < A_n_rows; ++row) - { - out_mem[row] = std::pow(A_mem[row], B_mem[row]); - } - } - } - #endif - } - else - { - for(uword i=0; i < A_n_cols; ++i) - { - const eT* A_mem = A.colptr(i); - eT* out_mem = out.colptr(i); - - for(uword row=0; row < A_n_rows; ++row) - { - out_mem[row] = std::pow(A_mem[row], B_mem[row]); - } - } - } - } - - if(mode == 1) // each row - { - if( arma_config::openmp && mp_gate::eval(A.n_elem) ) - { - #if defined(ARMA_USE_OPENMP) - { - const int n_threads = int( (std::min)(uword(mp_thread_limit::get()), A_n_cols) ); - - #pragma omp parallel for schedule(static) num_threads(n_threads) - for(uword i=0; i < A_n_cols; ++i) - { - const eT* A_mem = A.colptr(i); - eT* out_mem = out.colptr(i); - - const eT B_val = B_mem[i]; - - for(uword row=0; row < A_n_rows; ++row) - { - out_mem[row] = std::pow(A_mem[row], B_val); - } - } - } - #endif - } - else - { - for(uword i=0; i < A_n_cols; ++i) - { - const eT* A_mem = A.colptr(i); - eT* out_mem = out.colptr(i); - - const eT B_val = B_mem[i]; - - for(uword row=0; row < A_n_rows; ++row) - { - out_mem[row] = std::pow(A_mem[row], B_val); - } - } - } - } - - return out; - } - - - template inline void @@ -603,72 +375,4 @@ glue_powext_cx::apply(Cube< std::complex >& out, const Cube< std::complex -template -inline -Cube< std::complex > -glue_powext_cx::apply - ( - const subview_cube_each1< std::complex >& X, - const Base& Y - ) - { - arma_debug_sigprint(); - - typedef typename std::complex eT; - - const Cube& A = X.P; - - const uword A_n_rows = A.n_rows; - const uword A_n_cols = A.n_cols; - const uword A_n_slices = A.n_slices; - - Cube out(A_n_rows, A_n_cols, A_n_slices, arma_nozeros_indicator()); - - const quasi_unwrap tmp(Y.get_ref()); - const Mat& B = tmp.M; - - X.check_size(B); - - const T* B_mem = B.memptr(); - const uword B_n_elem = B.n_elem; - - if( arma_config::openmp && mp_gate::eval(A.n_elem) ) - { - #if defined(ARMA_USE_OPENMP) - { - const int n_threads = int( (std::min)(uword(mp_thread_limit::get()), A_n_slices) ); - - #pragma omp parallel for schedule(static) num_threads(n_threads) - for(uword s=0; s < A_n_slices; ++s) - { - const eT* A_slice_mem = A.slice_memptr(s); - eT* out_slice_mem = out.slice_memptr(s); - - for(uword i=0; i < B_n_elem; ++i) - { - out_slice_mem[i] = std::pow(A_slice_mem[i], B_mem[i]); - } - } - } - #endif - } - else - { - for(uword s=0; s < A_n_slices; ++s) - { - const eT* A_slice_mem = A.slice_memptr(s); - eT* out_slice_mem = out.slice_memptr(s); - - for(uword i=0; i < B_n_elem; ++i) - { - out_slice_mem[i] = std::pow(A_slice_mem[i], B_mem[i]); - } - } - } - - return out; - } - - - //! @}