From db9079914209d30289c17eced0d99c2130180318 Mon Sep 17 00:00:00 2001 From: Dirk Eddelbuettel Date: Mon, 10 Feb 2025 13:02:49 -0600 Subject: [PATCH 1/5] Use Armadillo 14.3.90 as 14.4.x candidate --- .../armadillo_bits/CubeToMatOp_bones.hpp | 4 +- inst/include/armadillo_bits/Cube_meat.hpp | 48 + inst/include/armadillo_bits/GenCube_bones.hpp | 6 +- inst/include/armadillo_bits/Gen_bones.hpp | 4 +- inst/include/armadillo_bits/Mat_meat.hpp | 48 + inst/include/armadillo_bits/OpCube_bones.hpp | 10 +- inst/include/armadillo_bits/Op_bones.hpp | 8 +- inst/include/armadillo_bits/Proxy.hpp | 106 +- inst/include/armadillo_bits/ProxyCube.hpp | 20 +- inst/include/armadillo_bits/SpOp_bones.hpp | 8 +- inst/include/armadillo_bits/SpProxy.hpp | 68 +- inst/include/armadillo_bits/SpToDOp_bones.hpp | 8 +- inst/include/armadillo_bits/arma_rng.hpp | 30 +- inst/include/armadillo_bits/arma_version.hpp | 6 +- inst/include/armadillo_bits/arrayops_meat.hpp | 11 +- inst/include/armadillo_bits/auxlib_meat.hpp | 4 + .../include/armadillo_bits/diagview_bones.hpp | 2 +- .../armadillo_bits/eGlueCube_bones.hpp | 6 +- inst/include/armadillo_bits/eGlue_bones.hpp | 4 +- inst/include/armadillo_bits/eOpCube_bones.hpp | 11 +- inst/include/armadillo_bits/eOp_bones.hpp | 9 +- inst/include/armadillo_bits/field_bones.hpp | 24 +- inst/include/armadillo_bits/field_meat.hpp | 246 +++- inst/include/armadillo_bits/fn_accu.hpp | 170 ++- inst/include/armadillo_bits/fn_elem.hpp | 184 ++- inst/include/armadillo_bits/fn_reshape.hpp | 49 + inst/include/armadillo_bits/fn_resize.hpp | 52 + .../armadillo_bits/glue_solve_meat.hpp | 7 + .../armadillo_bits/glue_times_meat.hpp | 154 +-- inst/include/armadillo_bits/memory.hpp | 6 +- .../armadillo_bits/mtGlueCube_bones.hpp | 6 +- inst/include/armadillo_bits/mtGlue_bones.hpp | 6 +- .../include/armadillo_bits/mtOpCube_bones.hpp | 14 +- inst/include/armadillo_bits/mtOp_bones.hpp | 12 +- inst/include/armadillo_bits/mtSpOp_bones.hpp | 10 +- .../armadillo_bits/mtSpReduceOp_bones.hpp | 6 +- inst/include/armadillo_bits/op_all_bones.hpp | 6 +- inst/include/armadillo_bits/op_all_meat.hpp | 153 ++- inst/include/armadillo_bits/op_any_bones.hpp | 6 +- inst/include/armadillo_bits/op_any_meat.hpp | 126 +- .../include/armadillo_bits/op_htrans_meat.hpp | 54 +- inst/include/armadillo_bits/op_mean_meat.hpp | 4 +- inst/include/armadillo_bits/op_misc_meat.hpp | 22 +- .../armadillo_bits/op_reshape_meat.hpp | 47 +- .../include/armadillo_bits/op_resize_meat.hpp | 36 +- inst/include/armadillo_bits/op_sort_bones.hpp | 10 +- .../armadillo_bits/op_sort_index_bones.hpp | 23 +- .../armadillo_bits/op_sort_index_meat.hpp | 71 +- inst/include/armadillo_bits/op_sort_meat.hpp | 229 ++-- .../include/armadillo_bits/op_strans_meat.hpp | 54 +- inst/include/armadillo_bits/op_sum_bones.hpp | 27 +- inst/include/armadillo_bits/op_sum_meat.hpp | 290 ++++- .../armadillo_bits/op_trimat_bones.hpp | 4 +- .../include/armadillo_bits/op_trimat_meat.hpp | 112 +- .../armadillo_bits/op_vectorise_meat.hpp | 63 +- .../armadillo_bits/spdiagview_bones.hpp | 2 +- .../include/armadillo_bits/spop_misc_meat.hpp | 16 +- inst/include/armadillo_bits/traits.hpp | 20 +- inst/include/armadillo_bits/unwrap.hpp | 1034 +---------------- 59 files changed, 1786 insertions(+), 1990 deletions(-) diff --git a/inst/include/armadillo_bits/CubeToMatOp_bones.hpp b/inst/include/armadillo_bits/CubeToMatOp_bones.hpp index a53dd2d0..07c88cae 100644 --- a/inst/include/armadillo_bits/CubeToMatOp_bones.hpp +++ b/inst/include/armadillo_bits/CubeToMatOp_bones.hpp @@ -33,8 +33,8 @@ class CubeToMatOp : public Base< typename T1::elem_type, CubeToMatOp constexpr bool is_alias(const Mat&) const { return false; } diff --git a/inst/include/armadillo_bits/Cube_meat.hpp b/inst/include/armadillo_bits/Cube_meat.hpp index 6ef8c271..b42a8613 100644 --- a/inst/include/armadillo_bits/Cube_meat.hpp +++ b/inst/include/armadillo_bits/Cube_meat.hpp @@ -2868,6 +2868,14 @@ Cube::Cube(const eOpCube& X) init_cold(); + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply(*this, reinterpret_cast< const eOpCube& >(X)); return; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply(*this, reinterpret_cast< const eOpCube& >(X)); return; } + } + eop_type::apply(*this, X); } @@ -2890,6 +2898,14 @@ Cube::operator=(const eOpCube& X) init_warm(X.get_n_rows(), X.get_n_cols(), X.get_n_slices()); + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply(*this, reinterpret_cast< const eOpCube& >(X)); return *this; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply(*this, reinterpret_cast< const eOpCube& >(X)); return *this; } + } + eop_type::apply(*this, X); return *this; @@ -2912,6 +2928,14 @@ Cube::operator+=(const eOpCube& X) if(bad_alias) { const Cube tmp(X); return (*this).operator+=(tmp); } + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply_inplace_plus(*this, reinterpret_cast< const eOpCube& >(X)); return *this; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply_inplace_plus(*this, reinterpret_cast< const eOpCube& >(X)); return *this; } + } + eop_type::apply_inplace_plus(*this, X); return *this; @@ -2934,6 +2958,14 @@ Cube::operator-=(const eOpCube& X) if(bad_alias) { const Cube tmp(X); return (*this).operator-=(tmp); } + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply_inplace_minus(*this, reinterpret_cast< const eOpCube& >(X)); return *this; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply_inplace_minus(*this, reinterpret_cast< const eOpCube& >(X)); return *this; } + } + eop_type::apply_inplace_minus(*this, X); return *this; @@ -2956,6 +2988,14 @@ Cube::operator%=(const eOpCube& X) if(bad_alias) { const Cube tmp(X); return (*this).operator%=(tmp); } + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply_inplace_schur(*this, reinterpret_cast< const eOpCube& >(X)); return *this; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply_inplace_schur(*this, reinterpret_cast< const eOpCube& >(X)); return *this; } + } + eop_type::apply_inplace_schur(*this, X); return *this; @@ -2978,6 +3018,14 @@ Cube::operator/=(const eOpCube& X) if(bad_alias) { const Cube tmp(X); return (*this).operator/=(tmp); } + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply_inplace_div(*this, reinterpret_cast< const eOpCube& >(X)); return *this; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply_inplace_div(*this, reinterpret_cast< const eOpCube& >(X)); return *this; } + } + eop_type::apply_inplace_div(*this, X); return *this; diff --git a/inst/include/armadillo_bits/GenCube_bones.hpp b/inst/include/armadillo_bits/GenCube_bones.hpp index 3c6099f3..67f326b5 100644 --- a/inst/include/armadillo_bits/GenCube_bones.hpp +++ b/inst/include/armadillo_bits/GenCube_bones.hpp @@ -33,9 +33,9 @@ class GenCube static constexpr bool use_at = false; static constexpr bool is_simple = (is_same_type::value) || (is_same_type::value); - arma_aligned const uword n_rows; - arma_aligned const uword n_cols; - arma_aligned const uword n_slices; + const uword n_rows; + const uword n_cols; + const uword n_slices; arma_inline GenCube(const uword in_n_rows, const uword in_n_cols, const uword in_n_slices); arma_inline ~GenCube(); diff --git a/inst/include/armadillo_bits/Gen_bones.hpp b/inst/include/armadillo_bits/Gen_bones.hpp index 352bfcdf..f843c91a 100644 --- a/inst/include/armadillo_bits/Gen_bones.hpp +++ b/inst/include/armadillo_bits/Gen_bones.hpp @@ -37,8 +37,8 @@ class Gen static constexpr bool is_col = T1::is_col; static constexpr bool is_xvec = T1::is_xvec; - arma_aligned const uword n_rows; - arma_aligned const uword n_cols; + const uword n_rows; + const uword n_cols; arma_inline Gen(const uword in_n_rows, const uword in_n_cols); arma_inline ~Gen(); diff --git a/inst/include/armadillo_bits/Mat_meat.hpp b/inst/include/armadillo_bits/Mat_meat.hpp index 93fb5c72..e18b1d74 100644 --- a/inst/include/armadillo_bits/Mat_meat.hpp +++ b/inst/include/armadillo_bits/Mat_meat.hpp @@ -5185,6 +5185,14 @@ Mat::Mat(const eOp& X) init_cold(); + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply(*this, reinterpret_cast< const eOp& >(X)); return; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply(*this, reinterpret_cast< const eOp& >(X)); return; } + } + eop_type::apply(*this, X); } @@ -5207,6 +5215,14 @@ Mat::operator=(const eOp& X) init_warm(X.get_n_rows(), X.get_n_cols()); + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply(*this, reinterpret_cast< const eOp& >(X)); return *this; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply(*this, reinterpret_cast< const eOp& >(X)); return *this; } + } + eop_type::apply(*this, X); return *this; @@ -5228,6 +5244,14 @@ Mat::operator+=(const eOp& X) if(bad_alias) { const Mat tmp(X); return (*this).operator+=(tmp); } + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply_inplace_plus(*this, reinterpret_cast< const eOp& >(X)); return *this; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply_inplace_plus(*this, reinterpret_cast< const eOp& >(X)); return *this; } + } + eop_type::apply_inplace_plus(*this, X); return *this; @@ -5249,6 +5273,14 @@ Mat::operator-=(const eOp& X) if(bad_alias) { const Mat tmp(X); return (*this).operator-=(tmp); } + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply_inplace_minus(*this, reinterpret_cast< const eOp& >(X)); return *this; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply_inplace_minus(*this, reinterpret_cast< const eOp& >(X)); return *this; } + } + eop_type::apply_inplace_minus(*this, X); return *this; @@ -5287,6 +5319,14 @@ Mat::operator%=(const eOp& X) if(bad_alias) { const Mat tmp(X); return (*this).operator%=(tmp); } + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply_inplace_schur(*this, reinterpret_cast< const eOp& >(X)); return *this; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply_inplace_schur(*this, reinterpret_cast< const eOp& >(X)); return *this; } + } + eop_type::apply_inplace_schur(*this, X); return *this; @@ -5308,6 +5348,14 @@ Mat::operator/=(const eOp& X) if(bad_alias) { const Mat tmp(X); return (*this).operator/=(tmp); } + if(is_same_type::value) + { + constexpr bool eT_non_int = is_non_integral::value; + + if( X.aux == eT(2) ) { eop_square::apply_inplace_div(*this, reinterpret_cast< const eOp& >(X)); return *this; } + if(eT_non_int && (X.aux == eT(0.5))) { eop_sqrt::apply_inplace_div(*this, reinterpret_cast< const eOp& >(X)); return *this; } + } + eop_type::apply_inplace_div(*this, X); return *this; diff --git a/inst/include/armadillo_bits/OpCube_bones.hpp b/inst/include/armadillo_bits/OpCube_bones.hpp index be4f618f..4b19aa4f 100644 --- a/inst/include/armadillo_bits/OpCube_bones.hpp +++ b/inst/include/armadillo_bits/OpCube_bones.hpp @@ -35,11 +35,11 @@ class OpCube : public BaseCube< typename T1::elem_type, OpCube > inline OpCube(const BaseCube& in_m, const uword in_aux_uword_a, const uword in_aux_uword_b, const uword in_aux_uword_c); inline ~OpCube(); - arma_aligned const T1& m; //!< the operand; must be derived from BaseCube - arma_aligned elem_type aux; //!< auxiliary data, using the element type as used by T1 - arma_aligned uword aux_uword_a; //!< auxiliary data, uword format - arma_aligned uword aux_uword_b; //!< auxiliary data, uword format - arma_aligned uword aux_uword_c; //!< auxiliary data, uword format + const T1& m; //!< the operand; must be derived from BaseCube + elem_type aux; //!< auxiliary data, using the element type as used by T1 + uword aux_uword_a; //!< auxiliary data, uword format + uword aux_uword_b; //!< auxiliary data, uword format + uword aux_uword_c; //!< auxiliary data, uword format }; diff --git a/inst/include/armadillo_bits/Op_bones.hpp b/inst/include/armadillo_bits/Op_bones.hpp index 7fc4088f..fe509bf7 100644 --- a/inst/include/armadillo_bits/Op_bones.hpp +++ b/inst/include/armadillo_bits/Op_bones.hpp @@ -61,10 +61,10 @@ class Op template inline bool is_alias(const Mat& X) const; - arma_aligned const T1& m; //!< the operand; must be derived from Base - arma_aligned elem_type aux; //!< auxiliary data, using the element type as used by T1 - arma_aligned uword aux_uword_a; //!< auxiliary data, uword format - arma_aligned uword aux_uword_b; //!< auxiliary data, uword format + const T1& m; //!< the operand; must be derived from Base + elem_type aux; //!< auxiliary data, using the element type as used by T1 + uword aux_uword_a; //!< auxiliary data, uword format + uword aux_uword_b; //!< auxiliary data, uword format }; diff --git a/inst/include/armadillo_bits/Proxy.hpp b/inst/include/armadillo_bits/Proxy.hpp index 441f68ab..0ddb7562 100644 --- a/inst/include/armadillo_bits/Proxy.hpp +++ b/inst/include/armadillo_bits/Proxy.hpp @@ -168,7 +168,7 @@ struct Proxy< Mat > static constexpr bool is_col = false; static constexpr bool is_xvec = false; - arma_aligned const Mat& Q; + const Mat& Q; inline explicit Proxy(const Mat& A) : Q(A) @@ -215,7 +215,7 @@ struct Proxy< Col > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const Col& Q; + const Col& Q; inline explicit Proxy(const Col& A) : Q(A) @@ -262,7 +262,7 @@ struct Proxy< Row > static constexpr bool is_col = false; static constexpr bool is_xvec = false; - arma_aligned const Row& Q; + const Row& Q; inline explicit Proxy(const Row& A) : Q(A) @@ -309,7 +309,7 @@ struct Proxy< Gen > static constexpr bool is_col = Gen::is_col; static constexpr bool is_xvec = Gen::is_xvec; - arma_aligned const Gen& Q; + const Gen& Q; inline explicit Proxy(const Gen& A) : Q(A) @@ -356,7 +356,7 @@ struct Proxy< eOp > static constexpr bool is_col = eOp::is_col; static constexpr bool is_xvec = eOp::is_xvec; - arma_aligned const eOp& Q; + const eOp& Q; inline explicit Proxy(const eOp& A) : Q(A) @@ -403,7 +403,7 @@ struct Proxy< eGlue > static constexpr bool is_col = eGlue::is_col; static constexpr bool is_xvec = eGlue::is_xvec; - arma_aligned const eGlue& Q; + const eGlue& Q; inline explicit Proxy(const eGlue& A) : Q(A) @@ -450,7 +450,7 @@ struct Proxy< Op > static constexpr bool is_col = Op::is_col; static constexpr bool is_xvec = Op::is_xvec; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const Op& A) : Q(A) @@ -497,7 +497,7 @@ struct Proxy< Glue > static constexpr bool is_col = Glue::is_col; static constexpr bool is_xvec = Glue::is_xvec; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const Glue& A) : Q(A) @@ -547,9 +547,9 @@ struct Proxy< Glue > static constexpr bool is_col = this_Glue_type::is_col; static constexpr bool is_xvec = this_Glue_type::is_xvec; - arma_aligned const this_Glue_type& Q; - arma_aligned const Proxy P1; - arma_aligned const Proxy P2; + const this_Glue_type& Q; + const Proxy P1; + const Proxy P2; arma_lt_comparator comparator; @@ -605,9 +605,9 @@ struct Proxy< Glue > static constexpr bool is_col = this_Glue_type::is_col; static constexpr bool is_xvec = this_Glue_type::is_xvec; - arma_aligned const this_Glue_type& Q; - arma_aligned const Proxy P1; - arma_aligned const Proxy P2; + const this_Glue_type& Q; + const Proxy P1; + const Proxy P2; arma_gt_comparator comparator; @@ -660,7 +660,7 @@ struct Proxy< mtOp > static constexpr bool is_col = mtOp::is_col; static constexpr bool is_xvec = mtOp::is_xvec; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const mtOp& A) : Q(A) @@ -707,7 +707,7 @@ struct Proxy< mtGlue > static constexpr bool is_col = mtGlue::is_col; static constexpr bool is_xvec = mtGlue::is_xvec; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const mtGlue& A) : Q(A) @@ -754,7 +754,7 @@ struct Proxy< CubeToMatOp > static constexpr bool is_col = CubeToMatOp::is_col; static constexpr bool is_xvec = CubeToMatOp::is_xvec; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const CubeToMatOp& A) : Q(A) @@ -801,8 +801,8 @@ struct Proxy< CubeToMatOp > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const unwrap_cube U; - arma_aligned const Mat Q; + const unwrap_cube U; + const Mat Q; inline explicit Proxy(const CubeToMatOp& A) : U(A.m) @@ -850,7 +850,7 @@ struct Proxy< SpToDOp > static constexpr bool is_col = SpToDOp::is_col; static constexpr bool is_xvec = SpToDOp::is_xvec; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const SpToDOp& A) : Q(A) @@ -897,8 +897,8 @@ struct Proxy< SpToDOp, op_sp_nonzeros> > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const SpMat& R; - arma_aligned const Mat Q; + const SpMat& R; + const Mat Q; inline explicit Proxy(const SpToDOp, op_sp_nonzeros>& A) : R(A.m) @@ -946,7 +946,7 @@ struct Proxy< SpToDGlue > static constexpr bool is_col = SpToDGlue::is_col; static constexpr bool is_xvec = SpToDGlue::is_xvec; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const SpToDGlue& A) : Q(A) @@ -993,7 +993,7 @@ struct Proxy< subview > static constexpr bool is_col = false; static constexpr bool is_xvec = false; - arma_aligned const subview& Q; + const subview& Q; inline explicit Proxy(const subview& A) : Q(A) @@ -1040,7 +1040,7 @@ struct Proxy< subview_col > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const subview_col& Q; + const subview_col& Q; inline explicit Proxy(const subview_col& A) : Q(A) @@ -1087,8 +1087,8 @@ struct Proxy< subview_cols > static constexpr bool is_col = false; static constexpr bool is_xvec = false; - arma_aligned const subview_cols& sv; - arma_aligned const Mat Q; + const subview_cols& sv; + const Mat Q; inline explicit Proxy(const subview_cols& A) : sv(A) @@ -1136,7 +1136,7 @@ struct Proxy< subview_row > static constexpr bool is_col = false; static constexpr bool is_xvec = false; - arma_aligned const subview_row& Q; + const subview_row& Q; inline explicit Proxy(const subview_row& A) : Q(A) @@ -1183,8 +1183,8 @@ struct Proxy< subview_elem1 > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const subview_elem1& Q; - arma_aligned const Proxy R; + const subview_elem1& Q; + const Proxy R; inline explicit Proxy(const subview_elem1& A) : Q(A) @@ -1237,7 +1237,7 @@ struct Proxy< subview_elem2 > static constexpr bool is_col = false; static constexpr bool is_xvec = false; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const subview_elem2& A) : Q(A) @@ -1284,7 +1284,7 @@ struct Proxy< diagview > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const diagview& Q; + const diagview& Q; inline explicit Proxy(const diagview& A) : Q(A) @@ -1339,8 +1339,8 @@ struct Proxy_diagvec_mat< Op > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const Mat& R; - arma_aligned const diagview Q; + const Mat& R; + const diagview Q; inline explicit Proxy_diagvec_mat(const Op& A) : R(A.m), Q( R.diag() ) @@ -1395,7 +1395,7 @@ struct Proxy_diagvec_expr< Op > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy_diagvec_expr(const Op& A) : Q(A) @@ -1468,7 +1468,7 @@ struct Proxy< Op > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const Op& A) : Q(A) @@ -1614,8 +1614,8 @@ struct Proxy_xtrans_vector< Op > static constexpr bool is_col = Op::is_col; static constexpr bool is_xvec = Op::is_xvec; - arma_aligned const quasi_unwrap U; // avoid copy if T1 is a Row, Col or subview_col - arma_aligned const Mat Q; + const quasi_unwrap U; // avoid copy if T1 is a Row, Col or subview_col + const Mat Q; inline Proxy_xtrans_vector(const Op& A) : U(A.m) @@ -1656,8 +1656,8 @@ struct Proxy_xtrans_vector< Op > static constexpr bool is_col = Op::is_col; static constexpr bool is_xvec = Op::is_xvec; - arma_aligned const quasi_unwrap U; // avoid copy if T1 is a Row, Col or subview_col - arma_aligned const Mat Q; + const quasi_unwrap U; // avoid copy if T1 is a Row, Col or subview_col + const Mat Q; inline Proxy_xtrans_vector(const Op& A) : U(A.m) @@ -1832,7 +1832,7 @@ struct Proxy_subview_row_htrans_cx static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const subview_row_htrans Q; + const subview_row_htrans Q; inline explicit Proxy_subview_row_htrans_cx(const Op, op_htrans>& A) : Q(A.m) @@ -1866,7 +1866,7 @@ struct Proxy_subview_row_htrans_non_cx static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const subview_row_strans Q; + const subview_row_strans Q; inline explicit Proxy_subview_row_htrans_non_cx(const Op, op_htrans>& A) : Q(A.m) @@ -1973,7 +1973,7 @@ struct Proxy< Op, op_strans> > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const subview_row_strans Q; + const subview_row_strans Q; inline explicit Proxy(const Op, op_strans>& A) : Q(A.m) @@ -2174,8 +2174,8 @@ struct Proxy< Op > static constexpr bool is_col = eOp< Op, eop_scalar_times>::is_col; static constexpr bool is_xvec = eOp< Op, eop_scalar_times>::is_xvec; - arma_aligned const Op R; - arma_aligned const eOp< Op, eop_scalar_times > Q; + const Op R; + const eOp< Op, eop_scalar_times > Q; inline explicit Proxy(const Op& A) : R(A.m) @@ -2223,7 +2223,7 @@ struct Proxy< subview_row_strans > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const subview_row_strans& Q; + const subview_row_strans& Q; inline explicit Proxy(const subview_row_strans& A) : Q(A) @@ -2270,7 +2270,7 @@ struct Proxy< subview_row_htrans > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const subview_row_htrans& Q; + const subview_row_htrans& Q; inline explicit Proxy(const subview_row_htrans& A) : Q(A) @@ -2317,7 +2317,7 @@ struct Proxy< xtrans_mat > static constexpr bool is_col = false; static constexpr bool is_xvec = false; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const xtrans_mat& A) : Q(A) @@ -2364,7 +2364,7 @@ struct Proxy< xvec_htrans > static constexpr bool is_col = false; static constexpr bool is_xvec = true; - arma_aligned const Mat Q; + const Mat Q; inline explicit Proxy(const xvec_htrans& A) : Q(A) @@ -2419,8 +2419,8 @@ struct Proxy_vectorise_col_mat< Op > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const unwrap U; - arma_aligned const Mat Q; + const unwrap U; + const Mat Q; inline explicit Proxy_vectorise_col_mat(const Op& A) : U(A.m) @@ -2476,8 +2476,8 @@ struct Proxy_vectorise_col_expr< Op > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const Op& Q; - arma_aligned const Proxy R; + const Op& Q; + const Proxy R; inline explicit Proxy_vectorise_col_expr(const Op& A) : Q(A) diff --git a/inst/include/armadillo_bits/ProxyCube.hpp b/inst/include/armadillo_bits/ProxyCube.hpp index 7f7cc0e0..3a34b5c9 100644 --- a/inst/include/armadillo_bits/ProxyCube.hpp +++ b/inst/include/armadillo_bits/ProxyCube.hpp @@ -48,7 +48,7 @@ struct ProxyCube< Cube > static constexpr bool use_mp = false; static constexpr bool has_subview = false; - arma_aligned const Cube& Q; + const Cube& Q; inline explicit ProxyCube(const Cube& A) : Q(A) @@ -93,7 +93,7 @@ struct ProxyCube< GenCube > static constexpr bool use_mp = false; static constexpr bool has_subview = false; - arma_aligned const GenCube& Q; + const GenCube& Q; inline explicit ProxyCube(const GenCube& A) : Q(A) @@ -138,7 +138,7 @@ struct ProxyCube< OpCube > static constexpr bool use_mp = false; static constexpr bool has_subview = false; - arma_aligned const Cube Q; + const Cube Q; inline explicit ProxyCube(const OpCube& A) : Q(A) @@ -183,7 +183,7 @@ struct ProxyCube< GlueCube > static constexpr bool use_mp = false; static constexpr bool has_subview = false; - arma_aligned const Cube Q; + const Cube Q; inline explicit ProxyCube(const GlueCube& A) : Q(A) @@ -228,7 +228,7 @@ struct ProxyCube< subview_cube > static constexpr bool use_mp = false; static constexpr bool has_subview = true; - arma_aligned const subview_cube& Q; + const subview_cube& Q; inline explicit ProxyCube(const subview_cube& A) : Q(A) @@ -273,7 +273,7 @@ struct ProxyCube< subview_cube_slices > static constexpr bool use_mp = false; static constexpr bool has_subview = false; - arma_aligned const Cube Q; + const Cube Q; inline explicit ProxyCube(const subview_cube_slices& A) : Q(A) @@ -318,7 +318,7 @@ struct ProxyCube< eOpCube > static constexpr bool use_mp = eOpCube::use_mp; static constexpr bool has_subview = eOpCube::has_subview; - arma_aligned const eOpCube& Q; + const eOpCube& Q; inline explicit ProxyCube(const eOpCube& A) : Q(A) @@ -363,7 +363,7 @@ struct ProxyCube< eGlueCube > static constexpr bool use_mp = eGlueCube::use_mp; static constexpr bool has_subview = eGlueCube::has_subview; - arma_aligned const eGlueCube& Q; + const eGlueCube& Q; inline explicit ProxyCube(const eGlueCube& A) : Q(A) @@ -408,7 +408,7 @@ struct ProxyCube< mtOpCube > static constexpr bool use_mp = false; static constexpr bool has_subview = false; - arma_aligned const Cube Q; + const Cube Q; inline explicit ProxyCube(const mtOpCube& A) : Q(A) @@ -453,7 +453,7 @@ struct ProxyCube< mtGlueCube > static constexpr bool use_mp = false; static constexpr bool has_subview = false; - arma_aligned const Cube Q; + const Cube Q; inline explicit ProxyCube(const mtGlueCube& A) : Q(A) diff --git a/inst/include/armadillo_bits/SpOp_bones.hpp b/inst/include/armadillo_bits/SpOp_bones.hpp index af8a229b..156584e9 100644 --- a/inst/include/armadillo_bits/SpOp_bones.hpp +++ b/inst/include/armadillo_bits/SpOp_bones.hpp @@ -40,10 +40,10 @@ class SpOp : public SpBase< typename T1::elem_type, SpOp > arma_inline bool is_alias(const SpMat& X) const; - arma_aligned const T1& m; //!< the operand; must be derived from SpBase - arma_aligned elem_type aux; //!< auxiliary data, using the element type as used by T1 - arma_aligned uword aux_uword_a; //!< auxiliary data, uword format - arma_aligned uword aux_uword_b; //!< auxiliary data, uword format + const T1& m; //!< the operand; must be derived from SpBase + elem_type aux; //!< auxiliary data, using the element type as used by T1 + uword aux_uword_a; //!< auxiliary data, uword format + uword aux_uword_b; //!< auxiliary data, uword format }; diff --git a/inst/include/armadillo_bits/SpProxy.hpp b/inst/include/armadillo_bits/SpProxy.hpp index 82dd7b6c..5b9f93cf 100644 --- a/inst/include/armadillo_bits/SpProxy.hpp +++ b/inst/include/armadillo_bits/SpProxy.hpp @@ -71,42 +71,42 @@ struct SpProxy< SpMat > typedef eT elem_type; typedef typename get_pod_type::result pod_type; typedef SpMat stored_type; - + typedef typename SpMat::const_iterator const_iterator_type; typedef typename SpMat::const_row_iterator const_row_iterator_type; - + static constexpr bool use_iterator = false; static constexpr bool Q_is_generated = false; - + static constexpr bool is_row = false; static constexpr bool is_col = false; static constexpr bool is_xvec = false; - - arma_aligned const SpMat& Q; - + + const SpMat& Q; + inline explicit SpProxy(const SpMat& A) : Q(A) { arma_debug_sigprint(); Q.sync(); } - + arma_inline uword get_n_rows() const { return Q.n_rows; } arma_inline uword get_n_cols() const { return Q.n_cols; } arma_inline uword get_n_elem() const { return Q.n_elem; } arma_inline uword get_n_nonzero() const { return Q.n_nonzero; } - + arma_inline elem_type operator[](const uword i) const { return Q[i]; } arma_inline elem_type at (const uword row, const uword col) const { return Q.at(row, col); } - + arma_inline const eT* get_values() const { return Q.values; } arma_inline const uword* get_row_indices() const { return Q.row_indices; } arma_inline const uword* get_col_ptrs() const { return Q.col_ptrs; } - + arma_inline const_iterator_type begin() const { return Q.begin(); } arma_inline const_iterator_type begin_col(const uword col_num) const { return Q.begin_col(col_num); } arma_inline const_row_iterator_type begin_row(const uword row_num = 0) const { return Q.begin_row(row_num); } - + arma_inline const_iterator_type end() const { return Q.end(); } arma_inline const_row_iterator_type end_row() const { return Q.end_row(); } arma_inline const_row_iterator_type end_row(const uword row_num) const { return Q.end_row(row_num); } @@ -134,7 +134,7 @@ struct SpProxy< SpCol > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const SpCol& Q; + const SpCol& Q; inline explicit SpProxy(const SpCol& A) : Q(A) @@ -186,7 +186,7 @@ struct SpProxy< SpRow > static constexpr bool is_col = false; static constexpr bool is_xvec = false; - arma_aligned const SpRow& Q; + const SpRow& Q; inline explicit SpProxy(const SpRow& A) : Q(A) @@ -227,42 +227,42 @@ struct SpProxy< SpSubview > typedef eT elem_type; typedef typename get_pod_type::result pod_type; typedef SpSubview stored_type; - + typedef typename SpSubview::const_iterator const_iterator_type; typedef typename SpSubview::const_row_iterator const_row_iterator_type; - + static constexpr bool use_iterator = true; static constexpr bool Q_is_generated = false; - + static constexpr bool is_row = false; static constexpr bool is_col = false; static constexpr bool is_xvec = false; - - arma_aligned const SpSubview& Q; - + + const SpSubview& Q; + inline explicit SpProxy(const SpSubview& A) : Q(A) { arma_debug_sigprint(); Q.m.sync(); } - + arma_inline uword get_n_rows() const { return Q.n_rows; } arma_inline uword get_n_cols() const { return Q.n_cols; } arma_inline uword get_n_elem() const { return Q.n_elem; } arma_inline uword get_n_nonzero() const { return Q.n_nonzero; } - + arma_inline elem_type operator[](const uword i) const { return Q[i]; } arma_inline elem_type at (const uword row, const uword col) const { return Q.at(row, col); } - + arma_inline const eT* get_values() const { return Q.m.values; } arma_inline const uword* get_row_indices() const { return Q.m.row_indices; } arma_inline const uword* get_col_ptrs() const { return Q.m.col_ptrs; } - + arma_inline const_iterator_type begin() const { return Q.begin(); } arma_inline const_iterator_type begin_col(const uword col_num) const { return Q.begin_col(col_num); } arma_inline const_row_iterator_type begin_row(const uword row_num = 0) const { return Q.begin_row(row_num); } - + arma_inline const_iterator_type end() const { return Q.end(); } arma_inline const_row_iterator_type end_row() const { return Q.end_row(); } arma_inline const_row_iterator_type end_row(const uword row_num) const { return Q.end_row(row_num); } @@ -290,7 +290,7 @@ struct SpProxy< SpSubview_col > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const SpSubview_col& Q; + const SpSubview_col& Q; inline explicit SpProxy(const SpSubview_col& A) : Q(A) @@ -342,7 +342,7 @@ struct SpProxy< SpSubview_col_list > static constexpr bool is_col = false; static constexpr bool is_xvec = false; - arma_aligned const SpMat Q; + const SpMat Q; inline explicit SpProxy(const SpSubview_col_list& A) : Q(A) @@ -385,7 +385,7 @@ struct SpProxy< SpSubview_row > typedef typename SpSubview::const_iterator const_iterator_type; typedef typename SpSubview::const_row_iterator const_row_iterator_type; - + static constexpr bool use_iterator = true; static constexpr bool Q_is_generated = false; @@ -393,7 +393,7 @@ struct SpProxy< SpSubview_row > static constexpr bool is_col = false; static constexpr bool is_xvec = false; - arma_aligned const SpSubview_row& Q; + const SpSubview_row& Q; inline explicit SpProxy(const SpSubview_row& A) : Q(A) @@ -445,7 +445,7 @@ struct SpProxy< spdiagview > static constexpr bool is_col = true; static constexpr bool is_xvec = false; - arma_aligned const SpMat Q; + const SpMat Q; inline explicit SpProxy(const spdiagview& A) : Q(A) @@ -497,7 +497,7 @@ struct SpProxy< SpOp > static constexpr bool is_col = SpOp::is_col; static constexpr bool is_xvec = SpOp::is_xvec; - arma_aligned const SpMat Q; + const SpMat Q; inline explicit SpProxy(const SpOp& A) : Q(A) @@ -549,7 +549,7 @@ struct SpProxy< SpGlue > static constexpr bool is_col = SpGlue::is_col; static constexpr bool is_xvec = SpGlue::is_xvec; - arma_aligned const SpMat Q; + const SpMat Q; inline explicit SpProxy(const SpGlue& A) : Q(A) @@ -600,7 +600,7 @@ struct SpProxy< mtSpOp > static constexpr bool is_col = mtSpOp::is_col; static constexpr bool is_xvec = mtSpOp::is_xvec; - arma_aligned const SpMat Q; + const SpMat Q; inline explicit SpProxy(const mtSpOp& A) : Q(A) @@ -651,7 +651,7 @@ struct SpProxy< mtSpGlue > static constexpr bool is_col = mtSpGlue::is_col; static constexpr bool is_xvec = mtSpGlue::is_xvec; - arma_aligned const SpMat Q; + const SpMat Q; inline explicit SpProxy(const mtSpGlue& A) : Q(A) @@ -702,7 +702,7 @@ struct SpProxy< mtSpReduceOp > static constexpr bool is_col = mtSpReduceOp::is_col; static constexpr bool is_xvec = mtSpReduceOp::is_xvec; - arma_aligned const SpMat Q; + const SpMat Q; inline explicit SpProxy(const mtSpReduceOp& A) : Q(A) diff --git a/inst/include/armadillo_bits/SpToDOp_bones.hpp b/inst/include/armadillo_bits/SpToDOp_bones.hpp index 2215c97a..957c062d 100644 --- a/inst/include/armadillo_bits/SpToDOp_bones.hpp +++ b/inst/include/armadillo_bits/SpToDOp_bones.hpp @@ -42,10 +42,10 @@ class SpToDOp : public Base< typename T1::elem_type, SpToDOp > template constexpr bool is_alias(const Mat&) const { return false; } - arma_aligned const T1& m; //!< the operand; must be derived from SpBase - arma_aligned elem_type aux; //!< auxiliary data, using the element type as used by T1 - arma_aligned uword aux_uword_a; //!< auxiliary data, uword format - arma_aligned uword aux_uword_b; //!< auxiliary data, uword format + const T1& m; //!< the operand; must be derived from SpBase + elem_type aux; //!< auxiliary data, using the element type as used by T1 + uword aux_uword_a; //!< auxiliary data, uword format + uword aux_uword_b; //!< auxiliary data, uword format }; diff --git a/inst/include/armadillo_bits/arma_rng.hpp b/inst/include/armadillo_bits/arma_rng.hpp index f1048810..21857f7c 100644 --- a/inst/include/armadillo_bits/arma_rng.hpp +++ b/inst/include/armadillo_bits/arma_rng.hpp @@ -267,20 +267,14 @@ arma_rng::set_seed_random() { try { - union - { - seed_type a; - unsigned char b[sizeof(seed_type)]; - } tmp; - - tmp.a = seed_type(0); + char tmp[sizeof(seed_type)] = {}; std::ifstream f("/dev/urandom", std::ifstream::binary); - if(f.good()) { f.read((char*)(&(tmp.b[0])), sizeof(seed_type)); } + if(f.good()) { f.read(&(tmp[0]), sizeof(seed_type)); } + + if(f.good()) { std::memcpy(&seed2, &(tmp[0]), sizeof(seed_type)); } - if(f.good()) { seed2 = tmp.a; } - have_seed = (seed2 != seed_type(0)); } catch(...) {} @@ -297,19 +291,17 @@ arma_rng::set_seed_random() seed3 = static_cast( since_epoch_usec & 0xFFFF ); - union - { - uword* a; - unsigned char b[sizeof(uword*)]; - } tmp; + unsigned char* a = (unsigned char*)std::malloc(std::size_t(4096)); - tmp.a = (uword*)malloc(sizeof(uword)); + unsigned char b[sizeof(unsigned char*)] = {}; - if(tmp.a != nullptr) + if(a != nullptr) { - for(size_t i=0; i::value) + { + const out_eT* src2 = (const out_eT*)src; + + arrayops::copy(dest, src2, n_elem); + + return; + } + uword j; for(j=1; j*, const std::complex*)) { + // TODO: investigate replacement of union-based conversion + union converter { blas_int (*fn)(const std::complex*, const std::complex*); diff --git a/inst/include/armadillo_bits/diagview_bones.hpp b/inst/include/armadillo_bits/diagview_bones.hpp index 4117aa45..0d5d40d5 100644 --- a/inst/include/armadillo_bits/diagview_bones.hpp +++ b/inst/include/armadillo_bits/diagview_bones.hpp @@ -29,7 +29,7 @@ class diagview : public Base< eT, diagview > typedef eT elem_type; typedef typename get_pod_type::result pod_type; - arma_aligned const Mat& m; + const Mat& m; static constexpr bool is_row = false; static constexpr bool is_col = true; diff --git a/inst/include/armadillo_bits/eGlueCube_bones.hpp b/inst/include/armadillo_bits/eGlueCube_bones.hpp index 8c157d8f..91733c72 100644 --- a/inst/include/armadillo_bits/eGlueCube_bones.hpp +++ b/inst/include/armadillo_bits/eGlueCube_bones.hpp @@ -27,13 +27,15 @@ class eGlueCube : public BaseCube< typename T1::elem_type, eGlueCube::result pod_type; + typedef ProxyCube proxy1_type; + typedef ProxyCube proxy2_type; static constexpr bool use_at = (ProxyCube::use_at || ProxyCube::use_at ); static constexpr bool use_mp = (ProxyCube::use_mp || ProxyCube::use_mp ); static constexpr bool has_subview = (ProxyCube::has_subview || ProxyCube::has_subview); - arma_aligned const ProxyCube P1; - arma_aligned const ProxyCube P2; + const ProxyCube P1; + const ProxyCube P2; arma_inline ~eGlueCube(); arma_inline eGlueCube(const T1& in_A, const T2& in_B); diff --git a/inst/include/armadillo_bits/eGlue_bones.hpp b/inst/include/armadillo_bits/eGlue_bones.hpp index a86377a9..412c6642 100644 --- a/inst/include/armadillo_bits/eGlue_bones.hpp +++ b/inst/include/armadillo_bits/eGlue_bones.hpp @@ -38,8 +38,8 @@ class eGlue : public Base< typename T1::elem_type, eGlue > static constexpr bool is_row = (Proxy::is_row || Proxy::is_row ); static constexpr bool is_xvec = (Proxy::is_xvec || Proxy::is_xvec); - arma_aligned const Proxy P1; - arma_aligned const Proxy P2; + const Proxy P1; + const Proxy P2; arma_inline ~eGlue(); arma_inline eGlue(const T1& in_A, const T2& in_B); diff --git a/inst/include/armadillo_bits/eOpCube_bones.hpp b/inst/include/armadillo_bits/eOpCube_bones.hpp index b6bcaba4..7732b6a0 100644 --- a/inst/include/armadillo_bits/eOpCube_bones.hpp +++ b/inst/include/armadillo_bits/eOpCube_bones.hpp @@ -28,16 +28,17 @@ class eOpCube : public BaseCube< typename T1::elem_type, eOpCube > typedef typename T1::elem_type elem_type; typedef typename get_pod_type::result pod_type; + typedef ProxyCube proxy_type; static constexpr bool use_at = ProxyCube::use_at; static constexpr bool use_mp = ProxyCube::use_mp || eop_type::use_mp; static constexpr bool has_subview = ProxyCube::has_subview; - arma_aligned const ProxyCube P; - arma_aligned elem_type aux; //!< storage of auxiliary data, user defined format - arma_aligned uword aux_uword_a; //!< storage of auxiliary data, uword format - arma_aligned uword aux_uword_b; //!< storage of auxiliary data, uword format - arma_aligned uword aux_uword_c; //!< storage of auxiliary data, uword format + const ProxyCube P; + elem_type aux; //!< storage of auxiliary data, user defined format + uword aux_uword_a; //!< storage of auxiliary data, uword format + uword aux_uword_b; //!< storage of auxiliary data, uword format + uword aux_uword_c; //!< storage of auxiliary data, uword format inline ~eOpCube(); inline explicit eOpCube(const BaseCube& in_m); diff --git a/inst/include/armadillo_bits/eOp_bones.hpp b/inst/include/armadillo_bits/eOp_bones.hpp index a200a6bb..3e18fa2b 100644 --- a/inst/include/armadillo_bits/eOp_bones.hpp +++ b/inst/include/armadillo_bits/eOp_bones.hpp @@ -38,11 +38,10 @@ class eOp : public Base< typename T1::elem_type, eOp > static constexpr bool is_col = Proxy::is_col; static constexpr bool is_xvec = Proxy::is_xvec; - arma_aligned const Proxy P; - - arma_aligned elem_type aux; //!< storage of auxiliary data, user defined format - arma_aligned uword aux_uword_a; //!< storage of auxiliary data, uword format - arma_aligned uword aux_uword_b; //!< storage of auxiliary data, uword format + const Proxy P; + elem_type aux; //!< storage of auxiliary data, user defined format + uword aux_uword_a; //!< storage of auxiliary data, uword format + uword aux_uword_b; //!< storage of auxiliary data, uword format inline ~eOp(); inline explicit eOp(const T1& in_m); diff --git a/inst/include/armadillo_bits/field_bones.hpp b/inst/include/armadillo_bits/field_bones.hpp index d3d2b053..c092a04b 100644 --- a/inst/include/armadillo_bits/field_bones.hpp +++ b/inst/include/armadillo_bits/field_bones.hpp @@ -21,13 +21,6 @@ -struct field_prealloc_n_elem - { - static constexpr uword val = 16; - }; - - - //! A lightweight 1D/2D/3D container for arbitrary objects //! (the objects must have a copy constructor) @@ -46,8 +39,7 @@ class field private: - arma_aligned oT** mem; //!< pointers to stored objects - arma_aligned oT* mem_local[ field_prealloc_n_elem::val ]; //!< local storage, for small fields + arma_aligned oT** mem; //!< pointers to stored objects public: @@ -67,12 +59,24 @@ class field inline explicit field(const SizeMat& s); inline explicit field(const SizeCube& s); - inline field& set_size(const uword n_obj_in); + inline field& set_size(const uword n_elem_in); inline field& set_size(const uword n_rows_in, const uword n_cols_in); inline field& set_size(const uword n_rows_in, const uword n_cols_in, const uword n_slices_in); inline field& set_size(const SizeMat& s); inline field& set_size(const SizeCube& s); + inline field& reshape(const uword n_elem_in); + inline field& reshape(const uword n_rows_in, const uword n_cols_in); + inline field& reshape(const uword n_rows_in, const uword n_cols_in, const uword n_slices_in); + inline field& reshape(const SizeMat& s); + inline field& reshape(const SizeCube& s); + + inline field& resize(const uword n_elem_in); + inline field& resize(const uword n_rows_in, const uword n_cols_in); + inline field& resize(const uword n_rows_in, const uword n_cols_in, const uword n_slices_in); + inline field& resize(const SizeMat& s); + inline field& resize(const SizeCube& s); + inline field(const std::vector& x); inline field& operator=(const std::vector& x); diff --git a/inst/include/armadillo_bits/field_meat.hpp b/inst/include/armadillo_bits/field_meat.hpp index fd1e9468..9b48add4 100644 --- a/inst/include/armadillo_bits/field_meat.hpp +++ b/inst/include/armadillo_bits/field_meat.hpp @@ -28,7 +28,7 @@ field::~field() delete_objects(); - if(n_elem > field_prealloc_n_elem::val) { delete [] mem; } + if(n_elem > 0) { delete [] mem; } // try to expose buggy user code that accesses deleted objects mem = nullptr; @@ -272,6 +272,173 @@ field::set_size(const SizeCube& s) +template +inline +field& +field::reshape(const uword n_elem_in) + { + arma_debug_sigprint(); + + return (*this).reshape(n_elem_in, 1, 1); + } + + + +template +inline +field& +field::reshape(const uword n_rows_in, const uword n_cols_in) + { + arma_debug_sigprint(); + + return (*this).reshape(n_rows_in, n_cols_in, 1); + } + + + +template +inline +field& +field::reshape(const uword n_rows_in, const uword n_cols_in, const uword n_slices_in) + { + arma_debug_sigprint(arma_str::format("n_rows_in: %u; n_cols_in: %u; n_slices_in: %u") % n_rows_in % n_cols_in % n_slices_in); + + if((n_rows == n_rows_in) && (n_cols == n_cols_in) && (n_slices == n_slices_in)) + { + // do nothing + } + else + if((n_elem == 0) || ((n_rows == n_cols_in) && (n_cols == n_rows_in) && (n_slices == n_slices_in))) + { + init(n_rows_in, n_cols_in, n_slices_in); + } + else + { + field tmp(n_rows_in, n_cols_in, n_slices_in); + + const uword n_elem_to_copy = (std::min)((*this).n_elem, tmp.n_elem); + + for(uword i=0; i < n_elem_to_copy; ++i) { tmp.at(i) = std::move((*this).at(i)); } + + (*this) = std::move(tmp); + } + + return *this; + } + + + +template +inline +field& +field::reshape(const SizeMat& s) + { + arma_debug_sigprint(); + + return (*this).reshape(s.n_rows, s.n_cols, 1); + } + + + +template +inline +field& +field::reshape(const SizeCube& s) + { + arma_debug_sigprint(); + + return (*this).reshape(s.n_rows, s.n_cols, s.n_slices); + } + + + +template +inline +field& +field::resize(const uword n_elem_in) + { + arma_debug_sigprint(); + + return (*this).resize(n_elem_in, 1, 1); + } + + + +template +inline +field& +field::resize(const uword n_rows_in, const uword n_cols_in) + { + arma_debug_sigprint(); + + return (*this).resize(n_rows_in, n_cols_in, 1); + } + + + +template +inline +field& +field::resize(const uword n_rows_in, const uword n_cols_in, const uword n_slices_in) + { + arma_debug_sigprint(arma_str::format("n_rows_in: %u; n_cols_in: %u; n_slices_in: %u") % n_rows_in % n_cols_in % n_slices_in); + + if((n_rows == n_rows_in) && (n_cols == n_cols_in) && (n_slices == n_slices_in)) + { + // do nothing + } + else + if(n_elem == 0) + { + (*this).set_size(n_rows_in, n_cols_in, n_slices_in); + } + else + { + // better-than-nothing implementation + + field tmp(n_rows_in, n_cols_in, n_slices_in); + + if(tmp.n_elem > 0) + { + const uword end_row = (std::min)(n_rows_in, n_rows ) - 1; + const uword end_col = (std::min)(n_cols_in, n_cols ) - 1; + const uword end_slice = (std::min)(n_slices_in, n_slices) - 1; + + tmp.subfield(0, 0, 0, end_row, end_col, end_slice) = (*this).subfield(0, 0, 0, end_row, end_col, end_slice); + } + + (*this) = std::move(tmp); + } + + return *this; + } + + + +template +inline +field& +field::resize(const SizeMat& s) + { + arma_debug_sigprint(); + + return (*this).resize(s.n_rows, s.n_cols, 1); + } + + + +template +inline +field& +field::resize(const SizeCube& s) + { + arma_debug_sigprint(); + + return (*this).resize(s.n_rows, s.n_cols, s.n_slices); + } + + + template inline field::field(const std::vector& x) @@ -279,6 +446,7 @@ field::field(const std::vector& x) , n_cols (0) , n_slices(0) , n_elem (0) + , mem (nullptr) { arma_debug_sigprint_this(this); @@ -312,6 +480,7 @@ field::field(const std::initializer_list& list) , n_cols (0) , n_slices(0) , n_elem (0) + , mem (nullptr) { arma_debug_sigprint_this(this); @@ -347,6 +516,7 @@ field::field(const std::initializer_list< std::initializer_list >& list) , n_cols (0) , n_slices(0) , n_elem (0) + , mem (nullptr) { arma_debug_sigprint_this(this); @@ -412,19 +582,10 @@ field::field(field&& X) , n_cols (X.n_cols ) , n_slices(X.n_slices) , n_elem (X.n_elem ) + , mem (X.mem ) { arma_debug_sigprint(arma_str::format("this: %x; X: %x") % this % &X); - if(n_elem > field_prealloc_n_elem::val) - { - mem = X.mem; - } - else - { - arrayops::copy(&mem_local[0], &X.mem_local[0], n_elem); - mem = mem_local; - } - access::rw(X.n_rows ) = 0; access::rw(X.n_cols ) = 0; access::rw(X.n_slices) = 0; @@ -450,15 +611,7 @@ field::operator=(field&& X) access::rw(n_slices) = X.n_slices; access::rw(n_elem ) = X.n_elem; - if(n_elem > field_prealloc_n_elem::val) - { - mem = X.mem; - } - else - { - arrayops::copy(&mem_local[0], &X.mem_local[0], n_elem); - mem = mem_local; - } + mem = X.mem; access::rw(X.n_rows ) = 0; access::rw(X.n_cols ) = 0; @@ -1967,34 +2120,15 @@ field::init(const field& x) { arma_debug_sigprint(); - if(this != &x) - { - const uword x_n_rows = x.n_rows; - const uword x_n_cols = x.n_cols; - const uword x_n_slices = x.n_slices; - - init(x_n_rows, x_n_cols, x_n_slices); - - field& t = *this; - - if(x_n_slices == 1) - { - for(uword ucol=0; ucol < x_n_cols; ++ucol) - for(uword urow=0; urow < x_n_rows; ++urow) - { - t.at(urow,ucol) = x.at(urow,ucol); - } - } - else - { - for(uword uslice=0; uslice < x_n_slices; ++uslice) - for(uword ucol=0; ucol < x_n_cols; ++ucol ) - for(uword urow=0; urow < x_n_rows; ++urow ) - { - t.at(urow,ucol,uslice) = x.at(urow,ucol,uslice); - } - } - } + if(this == &x) { return; } + + field& t = (*this); + + t.init(x.n_rows, x.n_cols, x.n_slices); + + const uword t_n_elem = t.n_elem; + + for(uword i=0; i < t_n_elem; ++i) { t.at(i) = x.at(i); } } @@ -2046,13 +2180,11 @@ field::init(const uword n_rows_in, const uword n_cols_in, const uword n_slic { delete_objects(); - if(n_elem > field_prealloc_n_elem::val) { delete [] mem; } + if(n_elem > 0) { delete [] mem; } - if(n_elem_new <= field_prealloc_n_elem::val) - { - mem = (n_elem_new == 0) ? nullptr : mem_local; - } - else + mem = nullptr; + + if(n_elem_new > 0) { mem = new(std::nothrow) oT* [n_elem_new]; @@ -2079,11 +2211,7 @@ field::delete_objects() for(uword i=0; i::value) || (is_subview_col::value) || (is_Mat::stored_type>::value)) + { + const quasi_unwrap U(X); + + return arrayops::accumulate(U.M.memptr(), U.M.n_elem); + } + const Proxy P(X); - if(is_Mat::stored_type>::value || is_subview_col::stored_type>::value) + return (Proxy::use_at) ? accu_proxy_at(P) : accu_proxy_linear(P); + } + + + +template +arma_warn_unused +inline +typename T1::elem_type +accu(const eOp& expr) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + typedef eOp expr_type; + + typedef typename expr_type::proxy_type::stored_type expr_P_stored_type; + + if((is_Mat::value) || (is_subview_col::value)) { - const quasi_unwrap::stored_type> tmp(P.Q); + const quasi_unwrap U(expr.P.Q); + + const eT* X_mem = U.M.memptr(); - return arrayops::accumulate(tmp.M.memptr(), tmp.M.n_elem); + return op_dot::direct_dot(U.M.n_elem, X_mem, X_mem); } - return (Proxy::use_at) ? accu_proxy_at(P) : accu_proxy_linear(P); + const Proxy P(expr); + + return (Proxy::use_at) ? accu_proxy_at(P) : accu_proxy_linear(P); + } + + + +template +arma_warn_unused +inline +typename T1::elem_type +accu(const eOp& expr) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + typedef eOp expr_type; + + if(expr.aux == eT(2)) + { + typedef eOp modified_expr_type; + + return accu( reinterpret_cast< const modified_expr_type& >(expr) ); + } + + if((expr.aux == eT(0.5)) && is_non_integral::value) + { + typedef eOp modified_expr_type; + + return accu( reinterpret_cast< const modified_expr_type& >(expr) ); + } + + const Proxy P(expr); + + return (Proxy::use_at) ? accu_proxy_at(P) : accu_proxy_linear(P); } @@ -836,16 +899,79 @@ accu(const BaseCube& X) { arma_debug_sigprint(); + if((is_Cube::value) || (is_Cube::stored_type>::value)) + { + const unwrap_cube U(X.get_ref()); + + return arrayops::accumulate(U.M.memptr(), U.M.n_elem); + } + const ProxyCube P(X.get_ref()); - if(is_Cube::stored_type>::value) + return (ProxyCube::use_at) ? accu_cube_proxy_at(P) : accu_cube_proxy_linear(P); + } + + + +template +arma_warn_unused +inline +typename T1::elem_type +accu(const eOpCube& expr) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + typedef eOpCube expr_type; + + typedef typename expr_type::proxy_type::stored_type expr_P_stored_type; + + if(is_Cube::value) { - const unwrap_cube::stored_type> tmp(P.Q); + const unwrap_cube U(expr.P.Q); + + const eT* X_mem = U.M.memptr(); - return arrayops::accumulate(tmp.M.memptr(), tmp.M.n_elem); + return op_dot::direct_dot(U.M.n_elem, X_mem, X_mem); } - return (ProxyCube::use_at) ? accu_cube_proxy_at(P) : accu_cube_proxy_linear(P); + const ProxyCube P(expr); + + return (ProxyCube::use_at) ? accu_cube_proxy_at(P) : accu_cube_proxy_linear(P); + } + + + +template +arma_warn_unused +inline +typename T1::elem_type +accu(const eOpCube& expr) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + typedef eOpCube expr_type; + + if(expr.aux == eT(2)) + { + typedef eOpCube modified_expr_type; + + return accu( reinterpret_cast< const modified_expr_type& >(expr) ); + } + + if((expr.aux == eT(0.5)) && is_non_integral::value) + { + typedef eOpCube modified_expr_type; + + return accu( reinterpret_cast< const modified_expr_type& >(expr) ); + } + + const ProxyCube P(expr); + + return (ProxyCube::use_at) ? accu_cube_proxy_at(P) : accu_cube_proxy_linear(P); } @@ -861,8 +987,8 @@ accu(const eGlueCube& expr) typedef eGlueCube expr_type; - typedef typename ProxyCube::stored_type P1_stored_type; - typedef typename ProxyCube::stored_type P2_stored_type; + typedef typename expr_type::proxy_type::stored_type P1_stored_type; + typedef typename expr_type::proxy_type::stored_type P2_stored_type; if(is_Cube::value && is_Cube::value) { @@ -1051,6 +1177,30 @@ accu(const SpOp& expr) if(is_vectorise) { return accu(expr.m); } + if(is_same_type::yes) + { + const SpProxy P(expr.m); + + const uword N = P.get_n_nonzero(); + + if(N == 0) { return eT(0); } + + if(SpProxy::use_iterator == false) + { + return op_dot::direct_dot(N, P.get_values(), P.get_values()); + } + else + { + typename SpProxy::const_iterator_type it = P.begin(); + + eT val = eT(0); + + for(uword i=0; i < N; ++i) { const eT tmp = (*it); ++it; val += tmp*tmp; } + + return val; + } + } + const SpMat tmp = expr; return accu(tmp); diff --git a/inst/include/armadillo_bits/fn_elem.hpp b/inst/include/armadillo_bits/fn_elem.hpp index 0a9649f6..ad4b0d11 100644 --- a/inst/include/armadillo_bits/fn_elem.hpp +++ b/inst/include/armadillo_bits/fn_elem.hpp @@ -52,12 +52,12 @@ real(const BaseCube& X) template arma_warn_unused arma_inline -const T1& -real(const SpBase& A) +typename enable_if2< (is_arma_sparse_type::value && is_cx::no), const T1& >::result +real(const T1& X) { arma_debug_sigprint(); - return A.get_ref(); + return X; } @@ -91,12 +91,12 @@ real(const BaseCube, T1>& X) template arma_warn_unused arma_inline -const mtSpOp -real(const SpBase,T1>& A) +typename enable_if2< (is_arma_sparse_type::value && is_cx::yes), const mtSpOp >::result +real(const T1& X) { arma_debug_sigprint(); - return mtSpOp(A.get_ref()); + return mtSpOp(X); } @@ -107,52 +107,7 @@ real(const SpBase,T1>& A) template arma_warn_unused inline -const Gen< Mat, gen_zeros > -imag(const Base& X) - { - arma_debug_sigprint(); - - const Proxy A(X.get_ref()); - - return Gen< Mat, gen_zeros>(A.get_n_rows(), A.get_n_cols()); - } - - - -template -arma_warn_unused -inline -const GenCube -imag(const BaseCube& X) - { - arma_debug_sigprint(); - - const ProxyCube A(X.get_ref()); - - return GenCube(A.get_n_rows(), A.get_n_cols(), A.get_n_slices()); - } - - - -template -arma_warn_unused -inline -SpMat -imag(const SpBase& A) - { - arma_debug_sigprint(); - - const SpProxy P(A.get_ref()); - - return SpMat(P.get_n_rows(), P.get_n_cols()); - } - - - -template -arma_warn_unused -inline -typename enable_if2< (is_arma_type::value && is_cx::yes), const mtOp >::result +typename enable_if2< is_arma_type::value, const mtOp >::result imag(const T1& X) { arma_debug_sigprint(); @@ -166,7 +121,7 @@ template arma_warn_unused inline const mtOpCube -imag(const BaseCube,T1>& X) +imag(const BaseCube& X) { arma_debug_sigprint(); @@ -178,12 +133,12 @@ imag(const BaseCube,T1>& X) template arma_warn_unused arma_inline -const mtSpOp -imag(const SpBase,T1>& A) +typename enable_if2< is_arma_sparse_type::value, const mtSpOp >::result +imag(const T1& X) { arma_debug_sigprint(); - return mtSpOp(A.get_ref()); + return mtSpOp(X); } @@ -480,27 +435,26 @@ abs(const BaseCube< std::complex,T1>& X, const typename a template arma_warn_unused arma_inline -const SpOp -abs(const SpBase& X, const typename arma_not_cx::result* junk = nullptr) +typename enable_if2< (is_arma_sparse_type::value && is_cx::no), const SpOp >::result +abs(const T1& X) { arma_debug_sigprint(); - arma_ignore(junk); - return SpOp(X.get_ref()); + return SpOp(X); } + template arma_warn_unused arma_inline -const mtSpOp -abs(const SpBase< std::complex, T1>& X, const typename arma_cx_only::result* junk = nullptr) +typename enable_if2< (is_arma_sparse_type::value && is_cx::yes), const mtSpOp >::result +abs(const T1& X) { arma_debug_sigprint(); - arma_ignore(junk); - return mtSpOp(X.get_ref()); + return mtSpOp(X); } @@ -568,13 +522,12 @@ arg(const BaseCube< std::complex,T1>& X, const typename a template arma_warn_unused arma_inline -const SpOp -arg(const SpBase& X, const typename arma_not_cx::result* junk = nullptr) +typename enable_if2< (is_arma_sparse_type::value && is_cx::no), const SpOp >::result +arg(const T1& X) { arma_debug_sigprint(); - arma_ignore(junk); - return SpOp(X.get_ref()); + return SpOp(X); } @@ -582,13 +535,12 @@ arg(const SpBase& X, const typename arma_not_cx arma_warn_unused arma_inline -const mtSpOp -arg(const SpBase< std::complex, T1>& X, const typename arma_cx_only::result* junk = nullptr) +typename enable_if2< (is_arma_sparse_type::value && is_cx::yes), const mtSpOp >::result +arg(const T1& X) { arma_debug_sigprint(); - arma_ignore(junk); - return mtSpOp(X.get_ref()); + return mtSpOp(X); } @@ -625,12 +577,12 @@ square(const BaseCube& A) template arma_warn_unused arma_inline -const SpOp -square(const SpBase& A) +typename enable_if2< is_arma_sparse_type::value, const SpOp >::result +square(const T1& X) { arma_debug_sigprint(); - return SpOp(A.get_ref()); + return SpOp(X); } @@ -667,12 +619,12 @@ sqrt(const BaseCube& A) template arma_warn_unused arma_inline -const SpOp -sqrt(const SpBase& A) +typename enable_if2< is_arma_sparse_type::value, const SpOp >::result +sqrt(const T1& X) { arma_debug_sigprint(); - return SpOp(A.get_ref()); + return SpOp(X); } @@ -709,12 +661,12 @@ cbrt(const BaseCube& A) template arma_warn_unused arma_inline -typename enable_if2< is_cx::no, const SpOp >::result -cbrt(const SpBase& A) +typename enable_if2< (is_arma_sparse_type::value && is_cx::no), const SpOp >::result +cbrt(const T1& X) { arma_debug_sigprint(); - return SpOp(A.get_ref()); + return SpOp(X); } @@ -725,12 +677,12 @@ cbrt(const SpBase& A) template arma_warn_unused arma_inline -const T1& -conj(const Base& A) +typename enable_if2< (is_arma_type::value && is_cx::no), const T1& >::result +conj(const T1& X) { arma_debug_sigprint(); - return A.get_ref(); + return X; } @@ -751,12 +703,12 @@ conj(const BaseCube& A) template arma_warn_unused arma_inline -const T1& -conj(const SpBase& A) +typename enable_if2< (is_arma_sparse_type::value && is_cx::no), const T1& >::result +conj(const T1& X) { arma_debug_sigprint(); - return A.get_ref(); + return X; } @@ -764,12 +716,12 @@ conj(const SpBase& A) template arma_warn_unused arma_inline -const eOp -conj(const Base,T1>& A) +typename enable_if2< (is_arma_type::value && is_cx::yes), const eOp >::result +conj(const T1& A) { arma_debug_sigprint(); - return eOp(A.get_ref()); + return eOp(A); } @@ -790,12 +742,12 @@ conj(const BaseCube,T1>& A) template arma_warn_unused arma_inline -const SpOp -conj(const SpBase,T1>& A) +typename enable_if2< (is_arma_sparse_type::value && is_cx::yes), const SpOp >::result +conj(const T1& X) { arma_debug_sigprint(); - return SpOp(A.get_ref()); + return SpOp(X); } @@ -805,12 +757,12 @@ conj(const SpBase,T1>& A) template arma_warn_unused arma_inline -const eOp -pow(const Base& A, const typename T1::elem_type exponent) +typename enable_if2< is_arma_type::value, const eOp >::result +pow(const T1& A, const typename T1::elem_type exponent) { arma_debug_sigprint(); - return eOp(A.get_ref(), exponent); + return eOp(A, exponent); } @@ -833,14 +785,14 @@ pow(const BaseCube& A, const typename T1::elem_type e template arma_warn_unused arma_inline -const eOp -pow(const Base& A, const typename T1::elem_type::value_type exponent) +typename enable_if2< (is_arma_type::value && is_cx::yes), const eOp >::result +pow(const T1& A, const typename T1::elem_type::value_type exponent) { arma_debug_sigprint(); typedef typename T1::elem_type eT; - return eOp(A.get_ref(), eT(exponent)); + return eOp(A, eT(exponent)); } @@ -849,11 +801,11 @@ template arma_warn_unused arma_inline const eOpCube -pow(const BaseCube& A, const typename T1::elem_type::value_type exponent) +pow(const BaseCube,T1>& A, const typename T1::elem_type::value_type exponent) { arma_debug_sigprint(); - typedef typename T1::elem_type eT; + typedef std::complex eT; return eOpCube(A.get_ref(), eT(exponent)); } @@ -892,12 +844,12 @@ floor(const BaseCube& A) template arma_warn_unused arma_inline -const SpOp -floor(const SpBase& X) +typename enable_if2< is_arma_sparse_type::value, const SpOp >::result +floor(const T1& X) { arma_debug_sigprint(); - return SpOp(X.get_ref()); + return SpOp(X); } @@ -934,12 +886,12 @@ ceil(const BaseCube& A) template arma_warn_unused arma_inline -const SpOp -ceil(const SpBase& X) +typename enable_if2< is_arma_sparse_type::value, const SpOp >::result +ceil(const T1& X) { arma_debug_sigprint(); - return SpOp(X.get_ref()); + return SpOp(X); } @@ -976,12 +928,12 @@ round(const BaseCube& A) template arma_warn_unused arma_inline -const SpOp -round(const SpBase& X) +typename enable_if2< is_arma_sparse_type::value, const SpOp >::result +round(const T1& X) { arma_debug_sigprint(); - return SpOp(X.get_ref()); + return SpOp(X); } @@ -1018,12 +970,12 @@ trunc(const BaseCube& A) template arma_warn_unused arma_inline -const SpOp -trunc(const SpBase& X) +typename enable_if2< is_arma_sparse_type::value, const SpOp >::result +trunc(const T1& X) { arma_debug_sigprint(); - return SpOp(X.get_ref()); + return SpOp(X); } @@ -1073,12 +1025,12 @@ sign(const BaseCube& A) template arma_warn_unused arma_inline -const SpOp -sign(const SpBase& X) +typename enable_if2< is_arma_sparse_type::value, const SpOp >::result +sign(const T1& X) { arma_debug_sigprint(); - return SpOp(X.get_ref()); + return SpOp(X); } diff --git a/inst/include/armadillo_bits/fn_reshape.hpp b/inst/include/armadillo_bits/fn_reshape.hpp index d64603f3..53fbdf38 100644 --- a/inst/include/armadillo_bits/fn_reshape.hpp +++ b/inst/include/armadillo_bits/fn_reshape.hpp @@ -135,4 +135,53 @@ reshape(const SpBase& X, const SizeMat& s) +// + + + +template +arma_warn_unused +inline +field +reshape(const field& A, const uword new_n_rows, const uword new_n_cols, const uword new_n_slices = uword(1)) + { + arma_debug_sigprint(); + + field B(new_n_rows, new_n_cols, new_n_slices); + + const uword n_elem_to_copy = (std::min)(A.n_elem, B.n_elem); + + for(uword i=0; i < n_elem_to_copy; ++i) { B.at(i) = A.at(i); } + + return B; + } + + + +template +arma_warn_unused +inline +field +reshape(const field& A, const SizeMat& s) + { + arma_debug_sigprint(); + + return reshape(A, s.n_rows, s.n_cols); + } + + + +template +arma_warn_unused +inline +field +reshape(const field& A, const SizeCube& s) + { + arma_debug_sigprint(); + + return reshape(A, s.n_rows, s.n_cols, s.n_slices); + } + + + //! @} diff --git a/inst/include/armadillo_bits/fn_resize.hpp b/inst/include/armadillo_bits/fn_resize.hpp index c61bd268..d1d0c060 100644 --- a/inst/include/armadillo_bits/fn_resize.hpp +++ b/inst/include/armadillo_bits/fn_resize.hpp @@ -99,4 +99,56 @@ resize(const SpBase& X, const SizeMat& s) +template +arma_warn_unused +inline +field +resize(const field& A, const uword in_n_rows, const uword in_n_cols, const uword in_n_slices = uword(1)) + { + arma_debug_sigprint(); + + // better-than-nothing implementation + + field B(in_n_rows, in_n_cols, in_n_slices); + + if((B.n_elem > 0) && (A.n_elem > 0)) + { + const uword end_row = (std::min)(in_n_rows, A.n_rows ) - 1; + const uword end_col = (std::min)(in_n_cols, A.n_cols ) - 1; + const uword end_slice = (std::min)(in_n_slices, A.n_slices) - 1; + + B.subfield(0, 0, 0, end_row, end_col, end_slice) = A.subfield(0, 0, 0, end_row, end_col, end_slice); + } + + return B; + } + + + +template +arma_warn_unused +inline +field +resize(const field& A, const SizeMat& s) + { + arma_debug_sigprint(); + + return resize(A, s.n_rows, s.n_cols); + } + + + +template +arma_warn_unused +inline +field +resize(const field& A, const SizeCube& s) + { + arma_debug_sigprint(); + + return resize(A, s.n_rows, s.n_cols, s.n_slices); + } + + + //! @} diff --git a/inst/include/armadillo_bits/glue_solve_meat.hpp b/inst/include/armadillo_bits/glue_solve_meat.hpp index 6aeb6729..d8031248 100644 --- a/inst/include/armadillo_bits/glue_solve_meat.hpp +++ b/inst/include/armadillo_bits/glue_solve_meat.hpp @@ -180,6 +180,13 @@ glue_solve_gen_full::apply(Mat& actual_out, const Base& A_expr, const const bool is_sym = arma_config::optimise_sym && ( (refine || equilibrate || likely_sympd || force_sym || is_band || is_triu || is_tril || auxlib::crippled_lapack(A)) ? false : is_sym_expr::eval(A_expr.get_ref()) ); const bool try_sympd = arma_config::optimise_sym && ( ( no_sympd || is_sym || force_sym || is_band || is_triu || is_tril || auxlib::crippled_lapack(A)) ? false : (likely_sympd ? true : sym_helper::guess_sympd(A, uword(16))) ); + arma_debug_print("glue_solve_gen_full::apply(): internal flags:"); + arma_debug_print("is_band: ", is_band ); + arma_debug_print("is_triu: ", is_triu ); + arma_debug_print("is_tril: ", is_tril ); + arma_debug_print("is_sym: ", is_sym ); + arma_debug_print("try_sympd: ", try_sympd); + if(fast) { // fast mode: solvers without refinement and without rcond estimate diff --git a/inst/include/armadillo_bits/glue_times_meat.hpp b/inst/include/armadillo_bits/glue_times_meat.hpp index 92d69599..28f5ec4b 100644 --- a/inst/include/armadillo_bits/glue_times_meat.hpp +++ b/inst/include/armadillo_bits/glue_times_meat.hpp @@ -31,20 +31,20 @@ glue_times_redirect2_helper::apply(Mat& o typedef typename T1::elem_type eT; - const partial_unwrap tmp1(X.A); - const partial_unwrap tmp2(X.B); + const partial_unwrap U1(X.A); + const partial_unwrap U2(X.B); - const typename partial_unwrap::stored_type& A = tmp1.M; - const typename partial_unwrap::stored_type& B = tmp2.M; + const typename partial_unwrap::stored_type& A = U1.M; + const typename partial_unwrap::stored_type& B = U2.M; constexpr bool use_alpha = partial_unwrap::do_times || partial_unwrap::do_times; - const eT alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val()) : eT(0); + const eT alpha = use_alpha ? (U1.get_val() * U2.get_val()) : eT(0); if( (is_cx::no) && (resolves_to_rowvector::value && resolves_to_colvector::value) ) { arma_debug_print("glue_times: dot product optimisation"); - arma_conform_assert_mul_size(A, B, tmp1.do_trans, tmp2.do_trans, "matrix multiplication"); + arma_conform_assert_mul_size(A, B, U1.do_trans, U2.do_trans, "matrix multiplication"); const eT val = op_dot::direct_dot(A.n_elem, A.memptr(), B.memptr()); @@ -55,7 +55,7 @@ glue_times_redirect2_helper::apply(Mat& o return; } - const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out); + const bool alias = U1.is_alias(out) || U2.is_alias(out); if(alias == false) { @@ -64,7 +64,7 @@ glue_times_redirect2_helper::apply(Mat& o eT, partial_unwrap::do_trans, partial_unwrap::do_trans, - (partial_unwrap::do_times || partial_unwrap::do_times) + use_alpha > (out, A, B, alpha); } @@ -77,7 +77,7 @@ glue_times_redirect2_helper::apply(Mat& o eT, partial_unwrap::do_trans, partial_unwrap::do_trans, - (partial_unwrap::do_times || partial_unwrap::do_times) + use_alpha > (tmp, A, B, alpha); @@ -186,18 +186,18 @@ glue_times_redirect3_helper::apply(Mat& o // we have exactly 3 objects // hence we can safely expand X as X.A.A, X.A.B and X.B - const partial_unwrap tmp1(X.A.A); - const partial_unwrap tmp2(X.A.B); - const partial_unwrap tmp3(X.B ); + const partial_unwrap U1(X.A.A); + const partial_unwrap U2(X.A.B); + const partial_unwrap U3(X.B ); - const typename partial_unwrap::stored_type& A = tmp1.M; - const typename partial_unwrap::stored_type& B = tmp2.M; - const typename partial_unwrap::stored_type& C = tmp3.M; + const typename partial_unwrap::stored_type& A = U1.M; + const typename partial_unwrap::stored_type& B = U2.M; + const typename partial_unwrap::stored_type& C = U3.M; constexpr bool use_alpha = partial_unwrap::do_times || partial_unwrap::do_times || partial_unwrap::do_times; - const eT alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val() * tmp3.get_val()) : eT(0); + const eT alpha = use_alpha ? (U1.get_val() * U2.get_val() * U3.get_val()) : eT(0); - const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out) || tmp3.is_alias(out); + const bool alias = U1.is_alias(out) || U2.is_alias(out) || U3.is_alias(out); if(alias == false) { @@ -207,7 +207,7 @@ glue_times_redirect3_helper::apply(Mat& o partial_unwrap::do_trans, partial_unwrap::do_trans, partial_unwrap::do_trans, - (partial_unwrap::do_times || partial_unwrap::do_times || partial_unwrap::do_times) + use_alpha > (out, A, B, C, alpha); } @@ -221,7 +221,7 @@ glue_times_redirect3_helper::apply(Mat& o partial_unwrap::do_trans, partial_unwrap::do_trans, partial_unwrap::do_trans, - (partial_unwrap::do_times || partial_unwrap::do_times || partial_unwrap::do_times) + use_alpha > (tmp, A, B, C, alpha); @@ -252,14 +252,14 @@ glue_times_redirect3_helper::apply(Mat& out, const arma_conform_check( (A.is_square() == false), "inv(): given matrix must be square sized" ); - const partial_unwrap tmp2(X.A.B); - const partial_unwrap tmp3(X.B ); + const partial_unwrap U2(X.A.B); + const partial_unwrap U3(X.B ); - const typename partial_unwrap::stored_type& B = tmp2.M; - const typename partial_unwrap::stored_type& C = tmp3.M; + const typename partial_unwrap::stored_type& B = U2.M; + const typename partial_unwrap::stored_type& C = U3.M; constexpr bool use_alpha = partial_unwrap::do_times || partial_unwrap::do_times; - const eT alpha = use_alpha ? (tmp2.get_val() * tmp3.get_val()) : eT(0); + const eT alpha = use_alpha ? (U2.get_val() * U3.get_val()) : eT(0); Mat BC; @@ -268,7 +268,7 @@ glue_times_redirect3_helper::apply(Mat& out, const eT, partial_unwrap::do_trans, partial_unwrap::do_trans, - (partial_unwrap::do_times || partial_unwrap::do_times) + use_alpha > (BC, B, C, alpha); @@ -306,8 +306,8 @@ glue_times_redirect3_helper::apply(Mat& out, const arma_conform_check( (B.is_square() == false), "inv(): given matrix must be square sized" ); - const unwrap C_tmp(X.B); - const Mat& C = C_tmp.M; + const quasi_unwrap U3(X.B); + const Mat& C = U3.M; arma_conform_assert_mul_size(B, C, "matrix multiplication"); @@ -330,21 +330,25 @@ glue_times_redirect3_helper::apply(Mat& out, const return; } - const partial_unwrap_check tmp1(X.A.A, out); + const partial_unwrap U1(X.A.A); - const typename partial_unwrap_check::stored_type& A = tmp1.M; + const typename partial_unwrap::stored_type& A = U1.M; - constexpr bool use_alpha = partial_unwrap_check::do_times; - const eT alpha = use_alpha ? tmp1.get_val() : eT(0); + constexpr bool use_alpha = partial_unwrap::do_times; + const eT alpha = use_alpha ? U1.get_val() : eT(0); - glue_times::apply - < - eT, - partial_unwrap_check::do_trans, - false, - partial_unwrap_check::do_times - > - (out, A, solve_result, alpha); + if(U1.is_alias(out)) + { + Mat tmp; + + glue_times::apply::do_trans, false, use_alpha>(tmp, A, solve_result, alpha); + + out.steal_mem(tmp); + } + else + { + glue_times::apply::do_trans, false, use_alpha>(out, A, solve_result, alpha); + } return; } @@ -365,16 +369,16 @@ glue_times_redirect::apply(Mat& out, const Glue tmp1(X.A); - const partial_unwrap tmp2(X.B); + const partial_unwrap U1(X.A); + const partial_unwrap U2(X.B); - const typename partial_unwrap::stored_type& A = tmp1.M; - const typename partial_unwrap::stored_type& B = tmp2.M; + const typename partial_unwrap::stored_type& A = U1.M; + const typename partial_unwrap::stored_type& B = U2.M; constexpr bool use_alpha = partial_unwrap::do_times || partial_unwrap::do_times; - const eT alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val()) : eT(0); + const eT alpha = use_alpha ? (U1.get_val() * U2.get_val()) : eT(0); - const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out); + const bool alias = U1.is_alias(out) || U2.is_alias(out); if(alias == false) { @@ -383,7 +387,7 @@ glue_times_redirect::apply(Mat& out, const Glue::do_trans, partial_unwrap::do_trans, - (partial_unwrap::do_times || partial_unwrap::do_times) + use_alpha > (out, A, B, alpha); } @@ -396,7 +400,7 @@ glue_times_redirect::apply(Mat& out, const Glue::do_trans, partial_unwrap::do_trans, - (partial_unwrap::do_times || partial_unwrap::do_times) + use_alpha > (tmp, A, B, alpha); @@ -446,20 +450,20 @@ glue_times_redirect<4>::apply(Mat& out, const Glue< Glue // there is exactly 4 objects // hence we can safely expand X as X.A.A.A, X.A.A.B, X.A.B and X.B - const partial_unwrap tmp1(X.A.A.A); - const partial_unwrap tmp2(X.A.A.B); - const partial_unwrap tmp3(X.A.B ); - const partial_unwrap tmp4(X.B ); + const partial_unwrap U1(X.A.A.A); + const partial_unwrap U2(X.A.A.B); + const partial_unwrap U3(X.A.B ); + const partial_unwrap U4(X.B ); - const typename partial_unwrap::stored_type& A = tmp1.M; - const typename partial_unwrap::stored_type& B = tmp2.M; - const typename partial_unwrap::stored_type& C = tmp3.M; - const typename partial_unwrap::stored_type& D = tmp4.M; + const typename partial_unwrap::stored_type& A = U1.M; + const typename partial_unwrap::stored_type& B = U2.M; + const typename partial_unwrap::stored_type& C = U3.M; + const typename partial_unwrap::stored_type& D = U4.M; constexpr bool use_alpha = partial_unwrap::do_times || partial_unwrap::do_times || partial_unwrap::do_times || partial_unwrap::do_times; - const eT alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val() * tmp3.get_val() * tmp4.get_val()) : eT(0); + const eT alpha = use_alpha ? (U1.get_val() * U2.get_val() * U3.get_val() * U4.get_val()) : eT(0); - const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out) || tmp3.is_alias(out) || tmp4.is_alias(out); + const bool alias = U1.is_alias(out) || U2.is_alias(out) || U3.is_alias(out) || U4.is_alias(out); if(alias == false) { @@ -470,7 +474,7 @@ glue_times_redirect<4>::apply(Mat& out, const Glue< Glue partial_unwrap::do_trans, partial_unwrap::do_trans, partial_unwrap::do_trans, - (partial_unwrap::do_times || partial_unwrap::do_times || partial_unwrap::do_times || partial_unwrap::do_times) + use_alpha > (out, A, B, C, D, alpha); } @@ -485,7 +489,7 @@ glue_times_redirect<4>::apply(Mat& out, const Glue< Glue partial_unwrap::do_trans, partial_unwrap::do_trans, partial_unwrap::do_trans, - (partial_unwrap::do_times || partial_unwrap::do_times || partial_unwrap::do_times || partial_unwrap::do_times) + use_alpha > (tmp, A, B, C, D, alpha); @@ -504,7 +508,7 @@ glue_times::apply(Mat& out, const Glue constexpr uword N_mat = 1 + depth_lhs< glue_times, Glue >::num; - arma_debug_print(arma_str::format("N_mat: %u") % N_mat); + arma_debug_print(arma_str::format("glue_times::apply(): N_mat: %u") % N_mat); glue_times_redirect::apply(out, X); } @@ -518,7 +522,11 @@ glue_times::apply_inplace(Mat& out, const T1& X) { arma_debug_sigprint(); - out = out * X; + typedef typename T1::elem_type eT; + + Mat tmp = out * X; + + out.steal_mem(tmp); } @@ -533,9 +541,9 @@ glue_times::apply_inplace_plus(Mat& out, const Glue::result T; - if( (is_outer_product::value) || (has_op_inv_any::value) || (has_op_inv_any::value) ) + if( X.is_alias(out) || (is_outer_product::value) || (has_op_inv_any::value) || (has_op_inv_any::value) ) { - // partial workaround for corner cases + // handle aliasing and partial workaround for corner cases const Mat tmp(X); @@ -544,21 +552,21 @@ glue_times::apply_inplace_plus(Mat& out, const Glue tmp1(X.A, out); - const partial_unwrap_check tmp2(X.B, out); + const partial_unwrap U1(X.A); + const partial_unwrap U2(X.B); - typedef typename partial_unwrap_check::stored_type TA; - typedef typename partial_unwrap_check::stored_type TB; + typedef typename partial_unwrap::stored_type TA; + typedef typename partial_unwrap::stored_type TB; - const TA& A = tmp1.M; - const TB& B = tmp2.M; + const TA& A = U1.M; + const TB& B = U2.M; - constexpr bool do_trans_A = partial_unwrap_check::do_trans; - constexpr bool do_trans_B = partial_unwrap_check::do_trans; + constexpr bool do_trans_A = partial_unwrap::do_trans; + constexpr bool do_trans_B = partial_unwrap::do_trans; - const bool use_alpha = partial_unwrap_check::do_times || partial_unwrap_check::do_times || (sign < sword(0)); + const bool use_alpha = partial_unwrap::do_times || partial_unwrap::do_times || (sign < sword(0)); - const eT alpha = use_alpha ? ( tmp1.get_val() * tmp2.get_val() * ( (sign > sword(0)) ? eT(1) : eT(-1) ) ) : eT(0); + const eT alpha = use_alpha ? ( U1.get_val() * U2.get_val() * ( (sign > sword(0)) ? eT(1) : eT(-1) ) ) : eT(0); arma_conform_assert_mul_size(A, B, do_trans_A, do_trans_B, "matrix multiplication"); diff --git a/inst/include/armadillo_bits/memory.hpp b/inst/include/armadillo_bits/memory.hpp index 34547e61..affbfb57 100644 --- a/inst/include/armadillo_bits/memory.hpp +++ b/inst/include/armadillo_bits/memory.hpp @@ -90,7 +90,7 @@ memory::acquire(const uword n_elem) #else { //return ( new(std::nothrow) eT[n_elem] ); - out_memptr = (eT *) malloc(sizeof(eT)*n_elem); + out_memptr = (eT *) std::malloc(sizeof(eT)*n_elem); } #endif @@ -124,7 +124,7 @@ memory::release(eT* mem) } #elif defined(ARMA_HAVE_POSIX_MEMALIGN) { - free( (void *)(mem) ); + std::free( (void *)(mem) ); } #elif defined(_MSC_VER) { @@ -134,7 +134,7 @@ memory::release(eT* mem) #else { //delete [] mem; - free( (void *)(mem) ); + std::free( (void *)(mem) ); } #endif diff --git a/inst/include/armadillo_bits/mtGlueCube_bones.hpp b/inst/include/armadillo_bits/mtGlueCube_bones.hpp index 846d8050..5512ec55 100644 --- a/inst/include/armadillo_bits/mtGlueCube_bones.hpp +++ b/inst/include/armadillo_bits/mtGlueCube_bones.hpp @@ -33,9 +33,9 @@ class mtGlueCube : public BaseCube< out_eT, mtGlueCube > template inline bool is_alias(const Mat& X) const; - arma_aligned const T1& A; //!< first operand; must be derived from Base - arma_aligned const T2& B; //!< second operand; must be derived from Base - arma_aligned uword aux_uword; //!< storage of auxiliary data, uword format + const T1& A; //!< first operand; must be derived from Base + const T2& B; //!< second operand; must be derived from Base + uword aux_uword; //!< storage of auxiliary data, uword format }; diff --git a/inst/include/armadillo_bits/mtOpCube_bones.hpp b/inst/include/armadillo_bits/mtOpCube_bones.hpp index ea9addfd..e793072c 100644 --- a/inst/include/armadillo_bits/mtOpCube_bones.hpp +++ b/inst/include/armadillo_bits/mtOpCube_bones.hpp @@ -45,14 +45,12 @@ class mtOpCube : public BaseCube< out_eT, mtOpCube > inline ~mtOpCube(); - - arma_aligned const T1& m; //!< the operand; must be derived from BaseCube - arma_aligned in_eT aux; //!< auxiliary data, using the element type as used by T1 - arma_aligned out_eT aux_out_eT; //!< auxiliary data, using the element type as specified by the out_eT template parameter - arma_aligned uword aux_uword_a; //!< auxiliary data, uword format - arma_aligned uword aux_uword_b; //!< auxiliary data, uword format - arma_aligned uword aux_uword_c; //!< auxiliary data, uword format - + const T1& m; //!< the operand; must be derived from BaseCube + in_eT aux; //!< auxiliary data, using the element type as used by T1 + out_eT aux_out_eT; //!< auxiliary data, using the element type as specified by the out_eT template parameter + uword aux_uword_a; //!< auxiliary data, uword format + uword aux_uword_b; //!< auxiliary data, uword format + uword aux_uword_c; //!< auxiliary data, uword format }; diff --git a/inst/include/armadillo_bits/mtOp_bones.hpp b/inst/include/armadillo_bits/mtOp_bones.hpp index 87e74573..cd73308a 100644 --- a/inst/include/armadillo_bits/mtOp_bones.hpp +++ b/inst/include/armadillo_bits/mtOp_bones.hpp @@ -51,13 +51,11 @@ class mtOp : public Base< out_eT, mtOp > template inline bool is_alias(const Mat& X) const; - - arma_aligned const T1& m; //!< the operand; must be derived from Base - arma_aligned in_eT aux; //!< auxiliary data, using the element type as used by T1 - arma_aligned out_eT aux_out_eT; //!< auxiliary data, using the element type as specified by the out_eT template parameter - arma_aligned uword aux_uword_a; //!< auxiliary data, uword format - arma_aligned uword aux_uword_b; //!< auxiliary data, uword format - + const T1& m; //!< the operand; must be derived from Base + in_eT aux; //!< auxiliary data, using the element type as used by T1 + out_eT aux_out_eT; //!< auxiliary data, using the element type as specified by the out_eT template parameter + uword aux_uword_a; //!< auxiliary data, uword format + uword aux_uword_b; //!< auxiliary data, uword format }; diff --git a/inst/include/armadillo_bits/mtSpOp_bones.hpp b/inst/include/armadillo_bits/mtSpOp_bones.hpp index 78facfaa..335b25d8 100644 --- a/inst/include/armadillo_bits/mtSpOp_bones.hpp +++ b/inst/include/armadillo_bits/mtSpOp_bones.hpp @@ -44,11 +44,11 @@ class mtSpOp : public SpBase< out_eT, mtSpOp > template arma_inline bool is_alias(const SpMat& X) const; - arma_aligned const T1& m; //!< the operand; must be derived from SpBase - arma_aligned in_eT aux; //!< auxiliary data, using the element type as used by T1 - arma_aligned out_eT aux_out_eT; //!< auxiliary data, using the element type as specified by the out_eT template parameter - arma_aligned uword aux_uword_a; - arma_aligned uword aux_uword_b; + const T1& m; //!< the operand; must be derived from SpBase + in_eT aux; //!< auxiliary data, using the element type as used by T1 + out_eT aux_out_eT; //!< auxiliary data, using the element type as specified by the out_eT template parameter + uword aux_uword_a; + uword aux_uword_b; }; diff --git a/inst/include/armadillo_bits/mtSpReduceOp_bones.hpp b/inst/include/armadillo_bits/mtSpReduceOp_bones.hpp index 51be0a3f..072c9976 100644 --- a/inst/include/armadillo_bits/mtSpReduceOp_bones.hpp +++ b/inst/include/armadillo_bits/mtSpReduceOp_bones.hpp @@ -52,9 +52,9 @@ class mtSpReduceOp : public SpBase< out_eT, mtSpReduceOp > inline mtSpReduceOp(const T1& in_m, const uword in_aux_uword_a, const uword in_aux_uword_b); inline ~mtSpReduceOp(); - arma_aligned const T1& m; //!< the operand; must be derived from SpBase - arma_aligned uword aux_uword_a; //!< auxiliary data, uword format - arma_aligned uword aux_uword_b; //!< auxiliary data, uword format + const T1& m; //!< the operand; must be derived from SpBase + uword aux_uword_a; //!< auxiliary data, uword format + uword aux_uword_b; //!< auxiliary data, uword format }; diff --git a/inst/include/armadillo_bits/op_all_bones.hpp b/inst/include/armadillo_bits/op_all_bones.hpp index b8faf9a2..6ebea6ba 100644 --- a/inst/include/armadillo_bits/op_all_bones.hpp +++ b/inst/include/armadillo_bits/op_all_bones.hpp @@ -68,9 +68,11 @@ class op_all static inline bool all_vec(T1& X); - template - static inline void apply_helper(Mat& out, const Proxy& P, const uword dim); + template + static inline void apply_mat_noalias(Mat& out, const Mat& X, const uword dim); + template + static inline void apply_proxy_noalias(Mat& out, const Proxy& P, const uword dim); template static inline void apply(Mat& out, const mtOp& X); diff --git a/inst/include/armadillo_bits/op_all_meat.hpp b/inst/include/armadillo_bits/op_all_meat.hpp index 185d4eec..a7503e46 100644 --- a/inst/include/armadillo_bits/op_all_meat.hpp +++ b/inst/include/armadillo_bits/op_all_meat.hpp @@ -277,17 +277,15 @@ op_all::all_vec(T1& X) -template +template inline void -op_all::apply_helper(Mat& out, const Proxy& P, const uword dim) +op_all::apply_mat_noalias(Mat& out, const Mat& X, const uword dim) { arma_debug_sigprint(); - const uword n_rows = P.get_n_rows(); - const uword n_cols = P.get_n_cols(); - - typedef typename Proxy::elem_type eT; + const uword n_rows = X.n_rows; + const uword n_cols = X.n_cols; if(dim == 0) // traverse rows (ie. process each column) { @@ -297,37 +295,18 @@ op_all::apply_helper(Mat& out, const Proxy& P, const uword dim) uword* out_mem = out.memptr(); - if(is_Mat::stored_type>::value) + for(uword col=0; col < n_cols; ++col) { - const unwrap::stored_type> U(P.Q); + const eT* colmem = X.colptr(col); - for(uword col=0; col < n_cols; ++col) - { - const eT* colmem = U.M.colptr(col); - - uword count = 0; - - for(uword row=0; row < n_rows; ++row) - { - count += (colmem[row] != eT(0)) ? uword(1) : uword(0); - } - - out_mem[col] = (n_rows == count) ? uword(1) : uword(0); - } - } - else - { - for(uword col=0; col < n_cols; ++col) + uword count = 0; + + for(uword row=0; row < n_rows; ++row) { - uword count = 0; - - for(uword row=0; row < n_rows; ++row) - { - if(P.at(row,col) != eT(0)) { ++count; } - } - - out_mem[col] = (n_rows == count) ? uword(1) : uword(0); + count += (colmem[row] != eT(0)) ? uword(1) : uword(0); } + + out_mem[col] = (n_rows == count) ? uword(1) : uword(0); } } else @@ -338,31 +317,72 @@ op_all::apply_helper(Mat& out, const Proxy& P, const uword dim) // internal dual use of 'out': keep the counts for each row - if(is_Mat::stored_type>::value) + for(uword col=0; col < n_cols; ++col) { - const unwrap::stored_type> U(P.Q); + const eT* colmem = X.colptr(col); - for(uword col=0; col < n_cols; ++col) + for(uword row=0; row < n_rows; ++row) { - const eT* colmem = U.M.colptr(col); - - for(uword row=0; row < n_rows; ++row) - { - out_mem[row] += (colmem[row] != eT(0)) ? uword(1) : uword(0); - } + out_mem[row] += (colmem[row] != eT(0)) ? uword(1) : uword(0); } } - else + + // see what the counts tell us + + for(uword row=0; row < n_rows; ++row) { - for(uword col=0; col < n_cols; ++col) + out_mem[row] = (n_cols == out_mem[row]) ? uword(1) : uword(0); + } + } + } + + + +template +inline +void +op_all::apply_proxy_noalias(Mat& out, const Proxy& P, const uword dim) + { + arma_debug_sigprint(); + + typedef typename Proxy::elem_type eT; + + const uword n_rows = P.get_n_rows(); + const uword n_cols = P.get_n_cols(); + + if(dim == 0) // traverse rows (ie. process each column) + { + out.zeros(1, n_cols); + + if(out.n_elem == 0) { return; } + + uword* out_mem = out.memptr(); + + for(uword col=0; col < n_cols; ++col) + { + uword count = 0; + + for(uword row=0; row < n_rows; ++row) { - for(uword row=0; row < n_rows; ++row) - { - if(P.at(row,col) != eT(0)) { ++out_mem[row]; } - } + if(P.at(row,col) != eT(0)) { ++count; } } + + out_mem[col] = (n_rows == count) ? uword(1) : uword(0); } + } + else + { + out.zeros(n_rows, 1); + + uword* out_mem = out.memptr(); + // internal dual use of 'out': keep the counts for each row + + for(uword col=0; col < n_cols; ++col) + for(uword row=0; row < n_rows; ++row) + { + if(P.at(row,col) != eT(0)) { ++out_mem[row]; } + } // see what the counts tell us @@ -370,7 +390,6 @@ op_all::apply_helper(Mat& out, const Proxy& P, const uword dim) { out_mem[row] = (n_cols == out_mem[row]) ? uword(1) : uword(0); } - } } @@ -385,19 +404,39 @@ op_all::apply(Mat& out, const mtOp& X) const uword dim = X.aux_uword_a; - const Proxy P(X.m); - - if(P.is_alias(out) == false) + if( (is_Mat::value) || (is_Mat::stored_type>::value) || (arma_config::openmp && Proxy::use_mp) ) { - op_all::apply_helper(out, P, dim); + const quasi_unwrap U(X.m); + + if(U.is_alias(out) == false) + { + op_all::apply_mat_noalias(out, U.M, dim); + } + else + { + Mat tmp; + + op_all::apply_mat_noalias(tmp, U.M, dim); + + out.steal_mem(tmp); + } } else { - Mat out2; - - op_all::apply_helper(out2, P, dim); + const Proxy P(X.m); - out.steal_mem(out2); + if(P.is_alias(out) == false) + { + op_all::apply_proxy_noalias(out, P, dim); + } + else + { + Mat tmp; + + op_all::apply_proxy_noalias(tmp, P, dim); + + out.steal_mem(tmp); + } } } diff --git a/inst/include/armadillo_bits/op_any_bones.hpp b/inst/include/armadillo_bits/op_any_bones.hpp index ffb197bd..6156dd48 100644 --- a/inst/include/armadillo_bits/op_any_bones.hpp +++ b/inst/include/armadillo_bits/op_any_bones.hpp @@ -68,9 +68,11 @@ class op_any static inline bool any_vec(T1& X); - template - static inline void apply_helper(Mat& out, const Proxy& P, const uword dim); + template + static inline void apply_mat_noalias(Mat& out, const Mat& X, const uword dim); + template + static inline void apply_proxy_noalias(Mat& out, const Proxy& P, const uword dim); template static inline void apply(Mat& out, const mtOp& X); diff --git a/inst/include/armadillo_bits/op_any_meat.hpp b/inst/include/armadillo_bits/op_any_meat.hpp index 3d5ff8aa..6e769b76 100644 --- a/inst/include/armadillo_bits/op_any_meat.hpp +++ b/inst/include/armadillo_bits/op_any_meat.hpp @@ -269,17 +269,15 @@ op_any::any_vec(T1& X) -template +template inline void -op_any::apply_helper(Mat& out, const Proxy& P, const uword dim) +op_any::apply_mat_noalias(Mat& out, const Mat& X, const uword dim) { arma_debug_sigprint(); - const uword n_rows = P.get_n_rows(); - const uword n_cols = P.get_n_cols(); - - typedef typename Proxy::elem_type eT; + const uword n_rows = X.n_rows; + const uword n_cols = X.n_cols; if(dim == 0) // traverse rows (ie. process each column) { @@ -287,59 +285,73 @@ op_any::apply_helper(Mat& out, const Proxy& P, const uword dim) uword* out_mem = out.memptr(); - if(is_Mat::stored_type>::value) + for(uword col=0; col < n_cols; ++col) { - const unwrap::stored_type> U(P.Q); + const eT* colmem = X.colptr(col); - for(uword col=0; col < n_cols; ++col) + for(uword row=0; row < n_rows; ++row) { - const eT* colmem = U.M.colptr(col); - - for(uword row=0; row < n_rows; ++row) - { - if(colmem[row] != eT(0)) { out_mem[col] = uword(1); break; } - } + if(colmem[row] != eT(0)) { out_mem[col] = uword(1); break; } } } - else + } + else + { + out.zeros(n_rows, 1); + + uword* out_mem = out.memptr(); + + for(uword col=0; col < n_cols; ++col) { - for(uword col=0; col < n_cols; ++col) + const eT* colmem = X.colptr(col); + + for(uword row=0; row < n_rows; ++row) { - for(uword row=0; row < n_rows; ++row) - { - if(P.at(row,col) != eT(0)) { out_mem[col] = uword(1); break; } - } + if(colmem[row] != eT(0)) { out_mem[row] = uword(1); } } } } - else + } + + + +template +inline +void +op_any::apply_proxy_noalias(Mat& out, const Proxy& P, const uword dim) + { + arma_debug_sigprint(); + + typedef typename Proxy::elem_type eT; + + const uword n_rows = P.get_n_rows(); + const uword n_cols = P.get_n_cols(); + + if(dim == 0) // traverse rows (ie. process each column) { - out.zeros(n_rows, 1); + out.zeros(1, n_cols); uword* out_mem = out.memptr(); - if(is_Mat::stored_type>::value) + for(uword col=0; col < n_cols; ++col) { - const unwrap::stored_type> U(P.Q); - - for(uword col=0; col < n_cols; ++col) + for(uword row=0; row < n_rows; ++row) { - const eT* colmem = U.M.colptr(col); - - for(uword row=0; row < n_rows; ++row) - { - if(colmem[row] != eT(0)) { out_mem[row] = uword(1); } - } + if(P.at(row,col) != eT(0)) { out_mem[col] = uword(1); break; } } } - else + } + else + { + out.zeros(n_rows, 1); + + uword* out_mem = out.memptr(); + + for(uword col=0; col < n_cols; ++col) { - for(uword col=0; col < n_cols; ++col) + for(uword row=0; row < n_rows; ++row) { - for(uword row=0; row < n_rows; ++row) - { - if(P.at(row,col) != eT(0)) { out_mem[row] = uword(1); } - } + if(P.at(row,col) != eT(0)) { out_mem[row] = uword(1); } } } } @@ -356,19 +368,39 @@ op_any::apply(Mat& out, const mtOp& X) const uword dim = X.aux_uword_a; - const Proxy P(X.m); - - if(P.is_alias(out) == false) + if( (is_Mat::value) || (is_Mat::stored_type>::value) || (arma_config::openmp && Proxy::use_mp) ) { - op_any::apply_helper(out, P, dim); + const quasi_unwrap U(X.m); + + if(U.is_alias(out) == false) + { + op_any::apply_mat_noalias(out, U.M, dim); + } + else + { + Mat tmp; + + op_any::apply_mat_noalias(tmp, U.M, dim); + + out.steal_mem(tmp); + } } else { - Mat out2; + const Proxy P(X.m); - op_any::apply_helper(out2, P, dim); - - out.steal_mem(out2); + if(P.is_alias(out) == false) + { + op_any::apply_proxy_noalias(out, P, dim); + } + else + { + Mat tmp; + + op_any::apply_proxy_noalias(tmp, P, dim); + + out.steal_mem(tmp); + } } } diff --git a/inst/include/armadillo_bits/op_htrans_meat.hpp b/inst/include/armadillo_bits/op_htrans_meat.hpp index 0968a50b..480d2900 100644 --- a/inst/include/armadillo_bits/op_htrans_meat.hpp +++ b/inst/include/armadillo_bits/op_htrans_meat.hpp @@ -306,49 +306,45 @@ op_htrans::apply_direct(Mat& out, const T1& X) typedef typename T1::elem_type eT; // allow detection of in-place transpose - if(is_Mat::value || (arma_config::openmp && Proxy::use_mp)) + if(is_Mat::value) { const unwrap U(X); op_htrans::apply_mat(out, U.M); } else + if((is_Mat::stored_type>::value) || (arma_config::openmp && Proxy::use_mp)) { - const Proxy P(X); + const quasi_unwrap U(X); - const bool is_alias = P.is_alias(out); + if(U.is_alias(out)) + { + Mat tmp; + + op_htrans::apply_mat_noalias(tmp, U.M); + + out.steal_mem(tmp); + } + else + { + op_htrans::apply_mat_noalias(out, U.M); + } + } + else + { + const Proxy P(X); - if(is_Mat::stored_type>::value) + if(P.is_alias(out)) { - const quasi_unwrap::stored_type> U(P.Q); + Mat tmp; - if(is_alias) - { - Mat tmp; - - op_htrans::apply_mat_noalias(tmp, U.M); - - out.steal_mem(tmp); - } - else - { - op_htrans::apply_mat_noalias(out, U.M); - } + op_htrans::apply_proxy(tmp, P); + + out.steal_mem(tmp); } else { - if(is_alias) - { - Mat tmp; - - op_htrans::apply_proxy(tmp, P); - - out.steal_mem(tmp); - } - else - { - op_htrans::apply_proxy(out, P); - } + op_htrans::apply_proxy(out, P); } } } diff --git a/inst/include/armadillo_bits/op_mean_meat.hpp b/inst/include/armadillo_bits/op_mean_meat.hpp index 3f06b263..1bb0c586 100644 --- a/inst/include/armadillo_bits/op_mean_meat.hpp +++ b/inst/include/armadillo_bits/op_mean_meat.hpp @@ -58,7 +58,7 @@ op_mean::apply_noalias(Mat& out, const Proxy& P, con { arma_debug_sigprint(); - if(is_Mat::stored_type>::value) + if((is_Mat::stored_type>::value) || (arma_config::openmp && Proxy::use_mp)) { op_mean::apply_noalias_unwrap(out, P, dim); } @@ -245,7 +245,7 @@ op_mean::apply_noalias(Cube& out, const ProxyCube& P { arma_debug_sigprint(); - if(is_Cube::stored_type>::value) + if((is_Cube::stored_type>::value) || (arma_config::openmp && ProxyCube::use_mp)) { op_mean::apply_noalias_unwrap(out, P, dim); } diff --git a/inst/include/armadillo_bits/op_misc_meat.hpp b/inst/include/armadillo_bits/op_misc_meat.hpp index 0f6a7b96..69c837ca 100644 --- a/inst/include/armadillo_bits/op_misc_meat.hpp +++ b/inst/include/armadillo_bits/op_misc_meat.hpp @@ -116,10 +116,15 @@ op_imag::apply( Mat& out, const mtOp P(X.m); + if(is_cx::no) { out.zeros(P.get_n_rows(), P.get_n_cols()); return; } + + // aliasing not possible at this point, as eT must be std::complex + const uword n_rows = P.get_n_rows(); const uword n_cols = P.get_n_cols(); @@ -136,7 +141,7 @@ op_imag::apply( Mat& out, const mtOp& out, const mtOp& out, const mtOpCube P(X.m); + if(is_cx::no) { out.zeros(P.get_n_rows(), P.get_n_cols(), P.get_n_slices()); return; } + + // aliasing not possible at this point, as eT must be std::complex + const uword n_rows = P.get_n_rows(); const uword n_cols = P.get_n_cols(); const uword n_slices = P.get_n_slices(); @@ -180,7 +190,7 @@ op_imag::apply( Cube& out, const mtOpCube& out, const mtOpCube inline void -op_reshape::apply(Mat& actual_out, const Op& in) +op_reshape::apply(Mat& out, const Op& in) { arma_debug_sigprint(); @@ -34,41 +34,54 @@ op_reshape::apply(Mat& actual_out, const Op::value || (arma_config::openmp && Proxy::use_mp)) + if(is_Mat::value) { const unwrap U(in.m); const Mat& A = U.M; - if(&actual_out == &A) + if(&out == &A) { - op_reshape::apply_mat_inplace(actual_out, new_n_rows, new_n_cols); + op_reshape::apply_mat_inplace(out, new_n_rows, new_n_cols); } else { - op_reshape::apply_mat_noalias(actual_out, A, new_n_rows, new_n_cols); + op_reshape::apply_mat_noalias(out, A, new_n_rows, new_n_cols); } } else + if((is_Mat::stored_type>::value) || (arma_config::openmp && Proxy::use_mp)) { - const Proxy P(in.m); - - const bool is_alias = P.is_alias(actual_out); - - Mat tmp; - Mat& out = (is_alias) ? tmp : actual_out; + const quasi_unwrap U(in.m); - if(is_Mat::stored_type>::value) + if(U.is_alias(out)) { - const quasi_unwrap::stored_type> U(P.Q); + Mat tmp; + op_reshape::apply_mat_noalias(tmp, U.M, new_n_rows, new_n_cols); + + out.steal_mem(tmp); + } + else + { op_reshape::apply_mat_noalias(out, U.M, new_n_rows, new_n_cols); } + } + else + { + const Proxy P(in.m); + + if(P.is_alias(out)) + { + Mat tmp; + + op_reshape::apply_proxy_noalias(tmp, P, new_n_rows, new_n_cols); + + out.steal_mem(tmp); + } else { op_reshape::apply_proxy_noalias(out, P, new_n_rows, new_n_cols); } - - if(is_alias) { actual_out.steal_mem(tmp); } } } @@ -95,7 +108,7 @@ op_reshape::apply_mat_inplace(Mat& A, const uword new_n_rows, const uword ne if(is_into_empty || is_into_colvec || is_into_rowvec || is_rowcol_swap) { A.set_size(new_n_rows, new_n_cols); return; } - Mat B; + Mat B(new_n_rows, new_n_cols, arma_nozeros_indicator()); op_reshape::apply_mat_noalias(B, A, new_n_rows, new_n_cols); @@ -226,7 +239,7 @@ op_reshape::apply_cube_inplace(Cube& A, const uword new_n_rows, const uword if(is_into_empty || is_into_colvec || is_into_rowvec || is_rowcol_swap) { A.set_size(new_n_rows, new_n_cols, new_n_slices); return; } - Cube B; + Cube B(new_n_rows, new_n_cols, new_n_slices, arma_nozeros_indicator()); op_reshape::apply_cube_noalias(B, A, new_n_rows, new_n_cols, new_n_slices); diff --git a/inst/include/armadillo_bits/op_resize_meat.hpp b/inst/include/armadillo_bits/op_resize_meat.hpp index 5db5d425..360d093b 100644 --- a/inst/include/armadillo_bits/op_resize_meat.hpp +++ b/inst/include/armadillo_bits/op_resize_meat.hpp @@ -34,16 +34,36 @@ op_resize::apply(Mat& out, const Op& in) const uword new_n_rows = in.aux_uword_a; const uword new_n_cols = in.aux_uword_b; - const unwrap tmp(in.m); - const Mat& A = tmp.M; - - if(&out == &A) + if(is_Mat::value) { - op_resize::apply_mat_inplace(out, new_n_rows, new_n_cols); + const unwrap U(in.m); + const Mat& A = U.M; + + if(&out == &A) + { + op_resize::apply_mat_inplace(out, new_n_rows, new_n_cols); + } + else + { + op_resize::apply_mat_noalias(out, A, new_n_rows, new_n_cols); + } } else { - op_resize::apply_mat_noalias(out, A, new_n_rows, new_n_cols); + const quasi_unwrap U(in.m); + + if(U.is_alias(out)) + { + Mat tmp; + + op_resize::apply_mat_noalias(tmp, U.M, new_n_rows, new_n_cols); + + out.steal_mem(tmp); + } + else + { + op_resize::apply_mat_noalias(out, U.M, new_n_rows, new_n_cols); + } } } @@ -63,7 +83,7 @@ op_resize::apply_mat_inplace(Mat& A, const uword new_n_rows, const uword new if(A.is_empty()) { A.zeros(new_n_rows, new_n_cols); return; } - Mat B; + Mat B(new_n_rows, new_n_cols, arma_nozeros_indicator()); op_resize::apply_mat_noalias(B, A, new_n_rows, new_n_cols); @@ -137,7 +157,7 @@ op_resize::apply_cube_inplace(Cube& A, const uword new_n_rows, const uword n if(A.is_empty()) { A.zeros(new_n_rows, new_n_cols, new_n_slices); return; } - Cube B; + Cube B(new_n_rows, new_n_cols, new_n_slices, arma_nozeros_indicator()); op_resize::apply_cube_noalias(B, A, new_n_rows, new_n_cols, new_n_slices); diff --git a/inst/include/armadillo_bits/op_sort_bones.hpp b/inst/include/armadillo_bits/op_sort_bones.hpp index 37449fe9..d21d9c31 100644 --- a/inst/include/armadillo_bits/op_sort_bones.hpp +++ b/inst/include/armadillo_bits/op_sort_bones.hpp @@ -27,19 +27,13 @@ class op_sort public: template - inline static void copy_row(eT* X, const Mat& A, const uword row); - - template - inline static void copy_row(Mat& A, const eT* X, const uword row); - - template - inline static void direct_sort(eT* X, const uword N, const uword sort_type = 0); + inline static void direct_sort(eT* X, const uword N, const uword sort_mode = 0); template inline static void direct_sort_ascending(eT* X, const uword N); template - inline static void apply_noalias(Mat& out, const Mat& X, const uword sort_type, const uword dim); + inline static void apply_noalias(Mat& out, const Mat& X, const uword sort_mode, const uword dim); template inline static void apply(Mat& out, const Op& in); diff --git a/inst/include/armadillo_bits/op_sort_index_bones.hpp b/inst/include/armadillo_bits/op_sort_index_bones.hpp index 7229ed5b..2eb9cd46 100644 --- a/inst/include/armadillo_bits/op_sort_index_bones.hpp +++ b/inst/include/armadillo_bits/op_sort_index_bones.hpp @@ -27,7 +27,10 @@ class op_sort_index public: template - static inline bool apply_noalias(Mat& out, const Proxy& P, const uword sort_type); + static inline bool apply_noalias_proxy(Mat& out, const Proxy& P, const uword sort_mode); + + template + static inline void apply_noalias_mat(Mat& out, const Mat& X, const uword sort_mode); template static inline void apply(Mat& out, const mtOp& in); @@ -41,7 +44,7 @@ class op_stable_sort_index public: template - static inline bool apply_noalias(Mat& out, const Proxy& P, const uword sort_type); + static inline bool apply_noalias(Mat& out, const Proxy& P, const uword sort_mode); template static inline void apply(Mat& out, const mtOp& in); @@ -134,4 +137,20 @@ struct arma_sort_index_helper_descend< std::complex > +template +struct arma_sort_index_helper_prepare + { + arma_inline eT operator() (const eT val) const { return val; } + }; + + + +template +struct arma_sort_index_helper_prepare< std::complex > + { + arma_inline T operator() (const std::complex& val) const { return std::abs(val); } + }; + + + //! @} diff --git a/inst/include/armadillo_bits/op_sort_index_meat.hpp b/inst/include/armadillo_bits/op_sort_index_meat.hpp index 1d02be7f..356dbd3b 100644 --- a/inst/include/armadillo_bits/op_sort_index_meat.hpp +++ b/inst/include/armadillo_bits/op_sort_index_meat.hpp @@ -24,17 +24,20 @@ template inline bool -arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_type) +arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_mode) { arma_debug_sigprint(); typedef typename T1::elem_type eT; + typedef typename T1::pod_type T; const uword n_elem = P.get_n_elem(); out.set_size(n_elem, 1); - std::vector< arma_sort_index_packet > packet_vec(n_elem); + const arma_sort_index_helper_prepare prepare; + + std::vector< arma_sort_index_packet > packet_vec(n_elem); if(Proxy::use_at == false) { @@ -44,7 +47,7 @@ arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_typ if(arma_isnan(val)) { out.soft_reset(); return false; } - packet_vec[i].val = val; + packet_vec[i].val = prepare(val); packet_vec[i].index = i; } } @@ -62,7 +65,7 @@ arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_typ if(arma_isnan(val)) { out.soft_reset(); return false; } - packet_vec[i].val = val; + packet_vec[i].val = prepare(val); packet_vec[i].index = i; ++i; @@ -70,11 +73,11 @@ arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_typ } - if(sort_type == 0) + if(sort_mode == 0) { // ascend - arma_sort_index_helper_ascend comparator; + arma_sort_index_helper_ascend comparator; if(sort_stable == false) { @@ -89,7 +92,7 @@ arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_typ { // descend - arma_sort_index_helper_descend comparator; + arma_sort_index_helper_descend comparator; if(sort_stable == false) { @@ -113,14 +116,36 @@ arma_sort_index_helper(Mat& out, const Proxy& P, const uword sort_typ +// + + + template inline bool -op_sort_index::apply_noalias(Mat& out, const Proxy& P, const uword sort_type) +op_sort_index::apply_noalias_proxy(Mat& out, const Proxy& P, const uword sort_mode) + { + arma_debug_sigprint(); + + return arma_sort_index_helper(out, P, sort_mode); + } + + + +template +inline +void +op_sort_index::apply_noalias_mat(Mat& out, const Mat& X, const uword sort_mode) { arma_debug_sigprint(); - return arma_sort_index_helper(out, P, sort_type); + if(X.n_elem == 0) { out.set_size(0,1); return; } + + const Proxy< Mat > P(X); + + const bool all_non_nan = op_sort_index::apply_noalias_proxy(out, P, sort_mode); + + arma_conform_check( (all_non_nan == false), "sort_index(): detected NaN" ); } @@ -136,21 +161,21 @@ op_sort_index::apply(Mat& out, const mtOp& in) if(P.get_n_elem() == 0) { out.set_size(0,1); return; } - const uword sort_type = in.aux_uword_a; + const uword sort_mode = in.aux_uword_a; bool all_non_nan = false; if(P.is_alias(out)) { - Mat out2; + Mat tmp; - all_non_nan = op_sort_index::apply_noalias(out2, P, sort_type); + all_non_nan = op_sort_index::apply_noalias_proxy(tmp, P, sort_mode); - out.steal_mem(out2); + out.steal_mem(tmp); } else { - all_non_nan = op_sort_index::apply_noalias(out, P, sort_type); + all_non_nan = op_sort_index::apply_noalias_proxy(out, P, sort_mode); } arma_conform_check( (all_non_nan == false), "sort_index(): detected NaN" ); @@ -158,14 +183,18 @@ op_sort_index::apply(Mat& out, const mtOp& in) +// + + + template inline bool -op_stable_sort_index::apply_noalias(Mat& out, const Proxy& P, const uword sort_type) +op_stable_sort_index::apply_noalias(Mat& out, const Proxy& P, const uword sort_mode) { arma_debug_sigprint(); - return arma_sort_index_helper(out, P, sort_type); + return arma_sort_index_helper(out, P, sort_mode); } @@ -181,21 +210,21 @@ op_stable_sort_index::apply(Mat& out, const mtOp out2; + Mat tmp; - all_non_nan = op_stable_sort_index::apply_noalias(out2, P, sort_type); + all_non_nan = op_stable_sort_index::apply_noalias(tmp, P, sort_mode); - out.steal_mem(out2); + out.steal_mem(tmp); } else { - all_non_nan = op_stable_sort_index::apply_noalias(out, P, sort_type); + all_non_nan = op_stable_sort_index::apply_noalias(out, P, sort_mode); } arma_conform_check( (all_non_nan == false), "stable_sort_index(): detected NaN" ); diff --git a/inst/include/armadillo_bits/op_sort_meat.hpp b/inst/include/armadillo_bits/op_sort_meat.hpp index 2e83055f..df3c40a5 100644 --- a/inst/include/armadillo_bits/op_sort_meat.hpp +++ b/inst/include/armadillo_bits/op_sort_meat.hpp @@ -24,11 +24,11 @@ template inline void -op_sort::direct_sort(eT* X, const uword n_elem, const uword sort_type) +op_sort::direct_sort(eT* X, const uword n_elem, const uword sort_mode) { arma_debug_sigprint(); - if(sort_type == 0) + if(sort_mode == 0) { arma_lt_comparator comparator; @@ -52,54 +52,8 @@ op_sort::direct_sort_ascending(eT* X, const uword n_elem) arma_debug_sigprint(); arma_lt_comparator comparator; - - std::sort(&X[0], &X[n_elem], comparator); - } - - - -template -inline -void -op_sort::copy_row(eT* X, const Mat& A, const uword row) - { - const uword N = A.n_cols; - - uword i,j; - - for(i=0, j=1; j -inline -void -op_sort::copy_row(Mat& A, const eT* X, const uword row) - { - const uword N = A.n_cols; - - uword i,j; - - for(i=0, j=1; j& A, const eT* X, const uword row) template inline void -op_sort::apply_noalias(Mat& out, const Mat& X, const uword sort_type, const uword dim) +op_sort::apply_noalias(Mat& out, const Mat& X, const uword sort_mode, const uword dim) { arma_debug_sigprint(); - if((X.n_rows * X.n_cols) <= 1) { out = X; return; } + if(X.n_elem <= 1) { out = X; return; } - if(dim == 0) // sort the contents of each column + if(is_cx::yes) { - arma_debug_print("op_sort::apply(): dim = 0"); + arma_debug_print("op_sort::apply(): complex version"); - out = X; - - const uword n_rows = out.n_rows; - const uword n_cols = out.n_cols; + if(dim == 0) // sort the contents of each column + { + arma_debug_print("op_sort::apply(): dim = 0"); + + const uword n_rows = X.n_rows; + const uword n_cols = X.n_cols; - for(uword col=0; col < n_cols; ++col) + out.set_size(n_rows, n_cols); + + uvec indices; + + for(uword col=0; col < n_cols; ++col) + { + const Col X_col( const_cast(X.colptr(col)), n_rows, false ); + + op_sort_index::apply_noalias_mat(indices, X_col, sort_mode); + + const uword* indices_mem = indices.memptr(); + const eT* X_col_mem = X_col.memptr(); + eT* out_col_mem = out.colptr(col); + + for(uword i=0; i < n_rows; ++i) { out_col_mem[i] = X_col_mem[ indices_mem[i] ]; } + } + } + else + if(dim == 1) // sort the contents of each row { - op_sort::direct_sort( out.colptr(col), n_rows, sort_type ); + arma_debug_print("op_sort::apply(): dim = 1"); + + Mat Y; + + op_strans::apply_mat_noalias(Y, X); + + const uword n_rows = Y.n_rows; + const uword n_cols = Y.n_cols; + + Mat tmp(n_rows, n_cols); + + uvec indices; + + for(uword col=0; col < n_cols; ++col) + { + const Col Y_col( const_cast(Y.colptr(col)), n_rows, false ); + + op_sort_index::apply_noalias_mat(indices, Y_col, sort_mode); + + const uword* indices_mem = indices.memptr(); + const eT* Y_col_mem = Y_col.memptr(); + eT* tmp_col_mem = tmp.colptr(col); + + for(uword i=0; i < n_rows; ++i) { tmp_col_mem[i] = Y_col_mem[ indices_mem[i] ]; } + } + + op_strans::apply_mat_noalias(out, tmp); } } else - if(dim == 1) // sort the contents of each row { - if(X.n_rows == 1) // a row vector + arma_debug_print("op_sort::apply(): plain version"); + + if(dim == 0) // sort the contents of each column { - arma_debug_print("op_sort::apply(): dim = 1, vector specific"); + arma_debug_print("op_sort::apply(): dim = 0"); out = X; - op_sort::direct_sort(out.memptr(), out.n_elem, sort_type); - } - else // not a row vector - { - arma_debug_print("op_sort::apply(): dim = 1, generic"); - - out.copy_size(X); const uword n_rows = out.n_rows; const uword n_cols = out.n_cols; - podarray tmp_array(n_cols); - - for(uword row=0; row < n_rows; ++row) + for(uword col=0; col < n_cols; ++col) + { + op_sort::direct_sort( out.colptr(col), n_rows, sort_mode ); + } + } + else + if(dim == 1) // sort the contents of each row + { + if(X.n_rows == 1) // a row vector { - op_sort::copy_row(tmp_array.memptr(), X, row); + arma_debug_print("op_sort::apply(): dim = 1, vector specific"); - op_sort::direct_sort( tmp_array.memptr(), n_cols, sort_type ); + out = X; - op_sort::copy_row(out, tmp_array.memptr(), row); + op_sort::direct_sort(out.memptr(), out.n_elem, sort_mode); + } + else // not a row vector + { + arma_debug_print("op_sort::apply(): dim = 1, generic"); + + Mat Y; + + op_strans::apply_mat_noalias(Y, X); + + const uword n_rows = Y.n_rows; + const uword n_cols = Y.n_cols; + + for(uword col=0; col < n_cols; ++col) + { + op_sort::direct_sort( Y.colptr(col), n_rows, sort_mode ); + } + + op_strans::apply_mat_noalias(out, Y); } } } @@ -174,10 +192,10 @@ op_sort::apply(Mat& out, const Op& in) const quasi_unwrap U(in.m); const Mat& X = U.M; - const uword sort_type = in.aux_uword_a; + const uword sort_mode = in.aux_uword_a; const uword dim = in.aux_uword_b; - arma_conform_check( (sort_type > 1), "sort(): parameter 'sort_type' must be 0 or 1" ); + arma_conform_check( (sort_mode > 1), "sort(): parameter 'sort_mode' must be 0 or 1" ); arma_conform_check( (dim > 1), "sort(): parameter 'dim' must be 0 or 1" ); arma_conform_check( (X.internal_has_nan()), "sort(): detected NaN" ); @@ -185,13 +203,13 @@ op_sort::apply(Mat& out, const Op& in) { Mat tmp; - op_sort::apply_noalias(tmp, X, sort_type, dim); + op_sort::apply_noalias(tmp, X, sort_mode, dim); out.steal_mem(tmp); } else { - op_sort::apply_noalias(out, X, sort_type, dim); + op_sort::apply_noalias(out, X, sort_mode, dim); } } @@ -209,31 +227,54 @@ op_sort_vec::apply(Mat& out, const Op& i const unwrap U(in.m); // not using quasi_unwrap, to ensure there is no aliasing with subviews const Mat& X = U.M; - const uword sort_type = in.aux_uword_a; + const uword sort_mode = in.aux_uword_a; - arma_conform_check( (sort_type > 1), "sort(): parameter 'sort_type' must be 0 or 1" ); + arma_conform_check( (sort_mode > 1), "sort(): parameter 'sort_mode' must be 0 or 1" ); arma_conform_check( (X.internal_has_nan()), "sort(): detected NaN" ); - out = X; // not checking for aliasing, to allow inplace sorting of vectors - - if(out.n_elem <= 1) { return; } + if(X.n_elem <= 1) { out = X; return; } - eT* out_mem = out.memptr(); - - eT* start_ptr = out_mem; - eT* endp1_ptr = &out_mem[out.n_elem]; - - if(sort_type == 0) + if(is_cx::yes) { - arma_lt_comparator comparator; + uvec indices; + + op_sort_index::apply_noalias_mat(indices, X, sort_mode); + + const uword N = indices.n_elem; + + arma_check( (N != X.n_elem), "internal error: op_sort_vec::apply(): N != X.n_elem" ); + + Mat tmp(X.n_rows, X.n_cols, arma_nozeros_indicator()); // in case there is aliasing + + const uword* indices_mem = indices.memptr(); + const eT* X_mem = X.memptr(); + eT* tmp_mem = tmp.memptr(); - std::sort(start_ptr, endp1_ptr, comparator); + for(uword i=0; i < N; ++i) { tmp_mem[i] = X_mem[ indices_mem[i] ]; } + + out.steal_mem(tmp); } else { - arma_gt_comparator comparator; + out = X; // not checking for aliasing, to allow inplace sorting of vectors - std::sort(start_ptr, endp1_ptr, comparator); + eT* out_mem = out.memptr(); + + eT* start_ptr = out_mem; + eT* endp1_ptr = &out_mem[out.n_elem]; + + if(sort_mode == 0) + { + arma_lt_comparator comparator; + + std::sort(start_ptr, endp1_ptr, comparator); + } + else + { + arma_gt_comparator comparator; + + std::sort(start_ptr, endp1_ptr, comparator); + } } } diff --git a/inst/include/armadillo_bits/op_strans_meat.hpp b/inst/include/armadillo_bits/op_strans_meat.hpp index 42d5f08d..7481ea6f 100644 --- a/inst/include/armadillo_bits/op_strans_meat.hpp +++ b/inst/include/armadillo_bits/op_strans_meat.hpp @@ -386,49 +386,45 @@ op_strans::apply_direct(Mat& out, const T1& X) typedef typename T1::elem_type eT; // allow detection of in-place transpose - if(is_Mat::value || (arma_config::openmp && Proxy::use_mp)) + if(is_Mat::value) { const unwrap U(X); op_strans::apply_mat(out, U.M); } else + if((is_Mat::stored_type>::value) || (is_subview_col::value) || (arma_config::openmp && Proxy::use_mp)) { - const Proxy P(X); + const quasi_unwrap U(X); - const bool is_alias = P.is_alias(out); + if(U.is_alias(out)) + { + Mat tmp; + + op_strans::apply_mat_noalias(tmp, U.M); + + out.steal_mem(tmp); + } + else + { + op_strans::apply_mat_noalias(out, U.M); + } + } + else + { + const Proxy P(X); - if(is_Mat::stored_type>::value) + if(P.is_alias(out)) { - const quasi_unwrap::stored_type> U(P.Q); + Mat tmp; - if(is_alias) - { - Mat tmp; - - op_strans::apply_mat_noalias(tmp, U.M); - - out.steal_mem(tmp); - } - else - { - op_strans::apply_mat_noalias(out, U.M); - } + op_strans::apply_proxy(tmp, P); + + out.steal_mem(tmp); } else { - if(is_alias) - { - Mat tmp; - - op_strans::apply_proxy(tmp, P); - - out.steal_mem(tmp); - } - else - { - op_strans::apply_proxy(out, P); - } + op_strans::apply_proxy(out, P); } } } diff --git a/inst/include/armadillo_bits/op_sum_bones.hpp b/inst/include/armadillo_bits/op_sum_bones.hpp index a72c7869..c580a204 100644 --- a/inst/include/armadillo_bits/op_sum_bones.hpp +++ b/inst/include/armadillo_bits/op_sum_bones.hpp @@ -28,31 +28,34 @@ class op_sum // dense matrices template - arma_hot inline static void apply(Mat& out, const Op& in); + inline static void apply(Mat& out, const Op< T1, op_sum >& in); template - arma_hot inline static void apply_noalias(Mat& out, const Proxy& P, const uword dim); - + inline static void apply(Mat& out, const Op< eOp, op_sum >& in); + template - arma_hot inline static void apply_noalias_unwrap(Mat& out, const Proxy& P, const uword dim); + inline static void apply(Mat& out, const Op< eOp, op_sum >& in); + + template + inline static void apply_mat_noalias(Mat& out, const Mat& X, const uword dim); + template + inline static void apply_mat_square_noalias(Mat& out, const Mat& X, const uword dim); + template - arma_hot inline static void apply_noalias_proxy(Mat& out, const Proxy& P, const uword dim); + inline static void apply_proxy_noalias(Mat& out, const Proxy& P, const uword dim); // cubes template - arma_hot inline static void apply(Cube& out, const OpCube& in); + inline static void apply(Cube& out, const OpCube& in); - template - arma_hot inline static void apply_noalias(Cube& out, const ProxyCube& P, const uword dim); - - template - arma_hot inline static void apply_noalias_unwrap(Cube& out, const ProxyCube& P, const uword dim); + template + inline static void apply_cube_noalias(Cube& out, const Cube& X, const uword dim); template - arma_hot inline static void apply_noalias_proxy(Cube& out, const ProxyCube& P, const uword dim); + inline static void apply_proxy_noalias(Cube& out, const ProxyCube& P, const uword dim); }; diff --git a/inst/include/armadillo_bits/op_sum_meat.hpp b/inst/include/armadillo_bits/op_sum_meat.hpp index 16becb36..d2f1eb8a 100644 --- a/inst/include/armadillo_bits/op_sum_meat.hpp +++ b/inst/include/armadillo_bits/op_sum_meat.hpp @@ -31,21 +31,42 @@ op_sum::apply(Mat& out, const Op& in) typedef typename T1::elem_type eT; const uword dim = in.aux_uword_a; - arma_conform_check( (dim > 1), "sum(): parameter 'dim' must be 0 or 1" ); - const Proxy P(in.m); + arma_conform_check( (dim > 1), "sum(): parameter 'dim' must be 0 or 1" ); - if(P.is_alias(out) == false) + if((is_Mat::value) || (is_Mat::stored_type>::value) || (arma_config::openmp && Proxy::use_mp)) { - op_sum::apply_noalias(out, P, dim); + const quasi_unwrap U(in.m); + + if(U.is_alias(out)) + { + Mat tmp; + + op_sum::apply_mat_noalias(tmp, U.M, dim); + + out.steal_mem(tmp); + } + else + { + op_sum::apply_mat_noalias(out, U.M, dim); + } } else { - Mat tmp; - - op_sum::apply_noalias(tmp, P, dim); + const Proxy P(in.m); - out.steal_mem(tmp); + if(P.is_alias(out)) + { + Mat tmp; + + op_sum::apply_proxy_noalias(tmp, P, dim); + + out.steal_mem(tmp); + } + else + { + op_sum::apply_proxy_noalias(out, P, dim); + } } } @@ -54,17 +75,71 @@ op_sum::apply(Mat& out, const Op& in) template inline void -op_sum::apply_noalias(Mat& out, const Proxy& P, const uword dim) +op_sum::apply(Mat& out, const Op< eOp, op_sum >& in) { arma_debug_sigprint(); - if(is_Mat::stored_type>::value || (arma_config::openmp && Proxy::use_mp)) + typedef typename T1::elem_type eT; + + typedef eOp inner_expr_type; + + typedef typename inner_expr_type::proxy_type::stored_type inner_expr_P_stored_type; + + const uword dim = in.aux_uword_a; + + arma_conform_check( (dim > 1), "sum(): parameter 'dim' must be 0 or 1" ); + + if(is_Mat::value) + { + const quasi_unwrap U(in.m.P.Q); + + if(U.is_alias(out)) + { + Mat tmp; + + op_sum::apply_mat_square_noalias(tmp, U.M, dim); + + out.steal_mem(tmp); + } + else + { + op_sum::apply_mat_square_noalias(out, U.M, dim); + } + } + else + if(arma_config::openmp && Proxy::use_mp) { - op_sum::apply_noalias_unwrap(out, P, dim); + const quasi_unwrap U(in.m); // force evaluation of compound inner expression + + if(U.is_alias(out)) + { + Mat tmp; + + op_sum::apply_mat_noalias(tmp, U.M, dim); + + out.steal_mem(tmp); + } + else + { + op_sum::apply_mat_noalias(out, U.M, dim); + } } else { - op_sum::apply_noalias_proxy(out, P, dim); + const Proxy P(in.m); + + if(P.is_alias(out)) + { + Mat tmp; + + op_sum::apply_proxy_noalias(tmp, P, dim); + + out.steal_mem(tmp); + } + else + { + op_sum::apply_proxy_noalias(out, P, dim); + } } } @@ -73,17 +148,82 @@ op_sum::apply_noalias(Mat& out, const Proxy& P, cons template inline void -op_sum::apply_noalias_unwrap(Mat& out, const Proxy& P, const uword dim) +op_sum::apply(Mat& out, const Op< eOp, op_sum >& in) { arma_debug_sigprint(); typedef typename T1::elem_type eT; - typedef typename Proxy::stored_type P_stored_type; + if(in.m.aux == eT(2)) + { + typedef Op< eOp, op_sum > modified_whole_expr_type; + + op_sum::apply(out, reinterpret_cast(in) ); + + return; + } + + if((in.m.aux == eT(0.5)) && is_non_integral::value) + { + typedef Op< eOp, op_sum > modified_whole_expr_type; + + op_sum::apply(out, reinterpret_cast(in) ); + + return; + } + + typedef eOp inner_expr_type; + + typedef typename inner_expr_type::proxy_type::stored_type inner_expr_P_stored_type; + + const uword dim = in.aux_uword_a; - const unwrap tmp(P.Q); + arma_conform_check( (dim > 1), "sum(): parameter 'dim' must be 0 or 1" ); - const typename unwrap::stored_type& X = tmp.M; + if( (is_Mat::value) || (arma_config::openmp && Proxy::use_mp) ) + { + const quasi_unwrap U(in.m); // force evaluation of eop_pow + + if(U.is_alias(out)) + { + Mat tmp; + + op_sum::apply_mat_noalias(tmp, U.M, dim); + + out.steal_mem(tmp); + } + else + { + op_sum::apply_mat_noalias(out, U.M, dim); + } + } + else + { + const Proxy P(in.m); + + if(P.is_alias(out)) + { + Mat tmp; + + op_sum::apply_proxy_noalias(tmp, P, dim); + + out.steal_mem(tmp); + } + else + { + op_sum::apply_proxy_noalias(out, P, dim); + } + } + } + + + +template +inline +void +op_sum::apply_mat_noalias(Mat& out, const Mat& X, const uword dim) + { + arma_debug_sigprint(); const uword X_n_rows = X.n_rows; const uword X_n_cols = X.n_cols; @@ -124,10 +264,56 @@ op_sum::apply_noalias_unwrap(Mat& out, const Proxy& +template +inline +void +op_sum::apply_mat_square_noalias(Mat& out, const Mat& X, const uword dim) + { + arma_debug_sigprint(); + + const uword X_n_rows = X.n_rows; + const uword X_n_cols = X.n_cols; + + const uword out_n_rows = (dim == 0) ? uword(1) : X_n_rows; + const uword out_n_cols = (dim == 0) ? X_n_cols : uword(1); + + out.set_size(out_n_rows, out_n_cols); + + if(X.n_elem == 0) { out.zeros(); return; } + + const eT* X_colptr = X.memptr(); + eT* out_mem = out.memptr(); + + if(dim == 0) + { + for(uword col=0; col < X_n_cols; ++col) + { + out_mem[col] = op_dot::direct_dot(X_n_rows, X_colptr, X_colptr); + + X_colptr += X_n_rows; + } + } + else + { + for(uword row=0; row < X_n_rows; ++row) { const eT tmp = X_colptr[row]; out_mem[row] = tmp*tmp; } + + X_colptr += X_n_rows; + + for(uword col=1; col < X_n_cols; ++col) + { + for(uword row=0; row < X_n_rows; ++row) { const eT tmp = X_colptr[row]; out_mem[row] += tmp*tmp; } + + X_colptr += X_n_rows; + } + } + } + + + template inline void -op_sum::apply_noalias_proxy(Mat& out, const Proxy& P, const uword dim) +op_sum::apply_proxy_noalias(Mat& out, const Proxy& P, const uword dim) { arma_debug_sigprint(); @@ -244,60 +430,54 @@ op_sum::apply(Cube& out, const OpCube& in) typedef typename T1::elem_type eT; const uword dim = in.aux_uword_a; - arma_conform_check( (dim > 2), "sum(): parameter 'dim' must be 0 or 1 or 2" ); - const ProxyCube P(in.m); + arma_conform_check( (dim > 2), "sum(): parameter 'dim' must be 0 or 1 or 2" ); - if(P.is_alias(out) == false) + if((is_Cube::value) || (is_Cube::stored_type>::value) || (arma_config::openmp && ProxyCube::use_mp)) { - op_sum::apply_noalias(out, P, dim); - } - else - { - Cube tmp; - - op_sum::apply_noalias(tmp, P, dim); + const unwrap_cube U(in.m); - out.steal_mem(tmp); - } - } - - - -template -inline -void -op_sum::apply_noalias(Cube& out, const ProxyCube& P, const uword dim) - { - arma_debug_sigprint(); - - if(is_Cube::stored_type>::value || (arma_config::openmp && ProxyCube::use_mp)) - { - op_sum::apply_noalias_unwrap(out, P, dim); + if(U.is_alias(out)) + { + Cube tmp; + + op_sum::apply_cube_noalias(tmp, U.M, dim); + + out.steal_mem(tmp); + } + else + { + op_sum::apply_cube_noalias(out, U.M, dim); + } } else { - op_sum::apply_noalias_proxy(out, P, dim); + const ProxyCube P(in.m); + + if(P.is_alias(out)) + { + Cube tmp; + + op_sum::apply_proxy_noalias(tmp, P, dim); + + out.steal_mem(tmp); + } + else + { + op_sum::apply_proxy_noalias(out, P, dim); + } } } -template +template inline void -op_sum::apply_noalias_unwrap(Cube& out, const ProxyCube& P, const uword dim) +op_sum::apply_cube_noalias(Cube& out, const Cube& X, const uword dim) { arma_debug_sigprint(); - typedef typename T1::elem_type eT; - - typedef typename ProxyCube::stored_type P_stored_type; - - const unwrap_cube tmp(P.Q); - - const Cube& X = tmp.M; - const uword X_n_rows = X.n_rows; const uword X_n_cols = X.n_cols; const uword X_n_slices = X.n_slices; @@ -350,7 +530,7 @@ op_sum::apply_noalias_unwrap(Cube& out, const ProxyCube< template inline void -op_sum::apply_noalias_proxy(Cube& out, const ProxyCube& P, const uword dim) +op_sum::apply_proxy_noalias(Cube& out, const ProxyCube& P, const uword dim) { arma_debug_sigprint(); diff --git a/inst/include/armadillo_bits/op_trimat_bones.hpp b/inst/include/armadillo_bits/op_trimat_bones.hpp index f500cbd4..de0137f5 100644 --- a/inst/include/armadillo_bits/op_trimat_bones.hpp +++ b/inst/include/armadillo_bits/op_trimat_bones.hpp @@ -37,10 +37,10 @@ class op_trimat inline static void apply(Mat& out, const Op& in); template - inline static void apply_unwrap(Mat& out, const Mat& A, const bool upper); + inline static void apply_mat_noalias(Mat& out, const Mat& A, const bool upper); template - inline static void apply_proxy(Mat& out, const Proxy& P, const bool upper); + inline static void apply_proxy_noalias(Mat& out, const Proxy& P, const bool upper); }; diff --git a/inst/include/armadillo_bits/op_trimat_meat.hpp b/inst/include/armadillo_bits/op_trimat_meat.hpp index 76b13d7c..404db924 100644 --- a/inst/include/armadillo_bits/op_trimat_meat.hpp +++ b/inst/include/armadillo_bits/op_trimat_meat.hpp @@ -68,49 +68,54 @@ op_trimat::apply(Mat& out, const Op& in) const bool upper = (in.aux_uword_a == 0); // allow detection of in-place operation - if(is_Mat::value || (arma_config::openmp && Proxy::use_mp)) + if(is_Mat::value) { const unwrap U(in.m); - op_trimat::apply_unwrap(out, U.M, upper); + if(&out == &(U.M)) + { + arma_conform_check( (U.M.is_square() == false), "trimatu()/trimatl(): given matrix must be square sized" ); + + op_trimat::fill_zeros(out, upper); + } + else + { + op_trimat::apply_mat_noalias(out, U.M, upper); + } } else + if((is_Mat::stored_type>::value) || (arma_config::openmp && Proxy::use_mp)) { - const Proxy P(in.m); + const quasi_unwrap U(in.m); - const bool is_alias = P.is_alias(out); + if(U.is_alias(out)) + { + Mat tmp; + + op_trimat::apply_mat_noalias(tmp, U.M, upper); + + out.steal_mem(tmp); + } + else + { + op_trimat::apply_mat_noalias(out, U.M, upper); + } + } + else + { + const Proxy P(in.m); - if(is_Mat::stored_type>::value) + if(P.is_alias(out)) { - const quasi_unwrap::stored_type> U(P.Q); + Mat tmp; - if(is_alias) - { - Mat tmp; - - op_trimat::apply_unwrap(tmp, U.M, upper); - - out.steal_mem(tmp); - } - else - { - op_trimat::apply_unwrap(out, U.M, upper); - } + op_trimat::apply_proxy_noalias(tmp, P, upper); + + out.steal_mem(tmp); } else { - if(is_alias) - { - Mat tmp; - - op_trimat::apply_proxy(tmp, P, upper); - - out.steal_mem(tmp); - } - else - { - op_trimat::apply_proxy(out, P, upper); - } + op_trimat::apply_proxy_noalias(out, P, upper); } } } @@ -120,39 +125,36 @@ op_trimat::apply(Mat& out, const Op& in) template inline void -op_trimat::apply_unwrap(Mat& out, const Mat& A, const bool upper) +op_trimat::apply_mat_noalias(Mat& out, const Mat& A, const bool upper) { arma_debug_sigprint(); arma_conform_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square sized" ); - if(&out != &A) + out.copy_size(A); + + const uword N = A.n_rows; + + if(upper) { - out.copy_size(A); - - const uword N = A.n_rows; - - if(upper) + // upper triangular: copy the diagonal and the elements above the diagonal + for(uword i=0; i& out, const Mat& A, const bool upper) template inline void -op_trimat::apply_proxy(Mat& out, const Proxy& P, const bool upper) +op_trimat::apply_proxy_noalias(Mat& out, const Proxy& P, const bool upper) { arma_debug_sigprint(); diff --git a/inst/include/armadillo_bits/op_vectorise_meat.hpp b/inst/include/armadillo_bits/op_vectorise_meat.hpp index 243c0819..2435e60e 100644 --- a/inst/include/armadillo_bits/op_vectorise_meat.hpp +++ b/inst/include/armadillo_bits/op_vectorise_meat.hpp @@ -44,7 +44,7 @@ op_vectorise_col::apply_direct(Mat& out, const T1& expr) typedef typename T1::elem_type eT; // allow detection of in-place operation - if(is_Mat::value || (arma_config::openmp && Proxy::use_mp)) + if(is_Mat::value) { const unwrap U(expr); @@ -80,42 +80,38 @@ op_vectorise_col::apply_direct(Mat& out, const T1& expr) } } else + if((is_Mat::stored_type>::value) || (arma_config::openmp && Proxy::use_mp)) { - const Proxy P(expr); + const quasi_unwrap U(expr); - const bool is_alias = P.is_alias(out); + if(U.is_alias(out)) + { + Mat tmp(U.M.memptr(), U.M.n_elem, 1); + + out.steal_mem(tmp); + } + else + { + out.set_size(U.M.n_elem, 1); + + arrayops::copy(out.memptr(), U.M.memptr(), U.M.n_elem); + } + } + else + { + const Proxy P(expr); - if(is_Mat::stored_type>::value) + if(P.is_alias(out)) { - const quasi_unwrap::stored_type> U(P.Q); + Mat tmp; - if(is_alias) - { - Mat tmp(U.M.memptr(), U.M.n_elem, 1); - - out.steal_mem(tmp); - } - else - { - out.set_size(U.M.n_elem, 1); - - arrayops::copy(out.memptr(), U.M.memptr(), U.M.n_elem); - } + op_vectorise_col::apply_proxy(tmp, P); + + out.steal_mem(tmp); } else { - if(is_alias) - { - Mat tmp; - - op_vectorise_col::apply_proxy(tmp, P); - - out.steal_mem(tmp); - } - else - { - op_vectorise_col::apply_proxy(out, P); - } + op_vectorise_col::apply_proxy(out, P); } } } @@ -343,7 +339,7 @@ op_vectorise_cube_col::apply(Mat& out, const CubeToMatOp } else { - if(is_Cube::value || (arma_config::openmp && ProxyCube::use_mp)) + if((is_Cube::value) || (is_Cube::stored_type>::value) || (arma_config::openmp && ProxyCube::use_mp)) { op_vectorise_cube_col::apply_unwrap(out, in.m); } @@ -409,13 +405,6 @@ op_vectorise_cube_col::apply_proxy(Mat& out, const T1& e const ProxyCube P(expr); - if(is_Cube::stored_type>::value) - { - op_vectorise_cube_col::apply_unwrap(out, P.Q); - - return; - } - const uword N = P.get_n_elem(); out.set_size(N, 1); diff --git a/inst/include/armadillo_bits/spdiagview_bones.hpp b/inst/include/armadillo_bits/spdiagview_bones.hpp index 238e8a38..c936dc0d 100644 --- a/inst/include/armadillo_bits/spdiagview_bones.hpp +++ b/inst/include/armadillo_bits/spdiagview_bones.hpp @@ -29,7 +29,7 @@ class spdiagview : public SpBase< eT, spdiagview > typedef eT elem_type; typedef typename get_pod_type::result pod_type; - arma_aligned const SpMat& m; + const SpMat& m; static constexpr bool is_row = false; static constexpr bool is_col = true; diff --git a/inst/include/armadillo_bits/spop_misc_meat.hpp b/inst/include/armadillo_bits/spop_misc_meat.hpp index 3560ffce..14fc9330 100644 --- a/inst/include/armadillo_bits/spop_misc_meat.hpp +++ b/inst/include/armadillo_bits/spop_misc_meat.hpp @@ -288,8 +288,11 @@ namespace priv { struct functor_imag { + template + arma_inline eT operator()(const eT ) const { return eT(0); } + template - arma_inline T operator()(const std::complex& val) const { return val.imag(); } + arma_inline T operator()(const std::complex& val) const { return val.imag(); } }; } @@ -302,7 +305,16 @@ spop_imag::apply(SpMat& out, const mtSpOp::no) + { + const SpProxy P(in.m); + + out.zeros(P.get_n_rows(), P.get_n_cols()); + } + else + { + out.init_xform_mt(in.m, priv::functor_imag()); + } } diff --git a/inst/include/armadillo_bits/traits.hpp b/inst/include/armadillo_bits/traits.hpp index e3bf703e..e1e2d5ed 100644 --- a/inst/include/armadillo_bits/traits.hpp +++ b/inst/include/armadillo_bits/traits.hpp @@ -33,8 +33,8 @@ struct get_pod_type< std::complex > template struct is_Mat_fixed_only { - typedef char yes[1]; - typedef char no[2]; + using yes = char[1]; + using no = char[2]; template static yes& check(typename X::Mat_fixed_type*); template static no& check(...); @@ -47,8 +47,8 @@ struct is_Mat_fixed_only template struct is_Row_fixed_only { - typedef char yes[1]; - typedef char no[2]; + using yes = char[1]; + using no = char[2]; template static yes& check(typename X::Row_fixed_type*); template static no& check(...); @@ -61,8 +61,8 @@ struct is_Row_fixed_only template struct is_Col_fixed_only { - typedef char yes[1]; - typedef char no[2]; + using yes = char[1]; + using no = char[2]; template static yes& check(typename X::Col_fixed_type*); template static no& check(...); @@ -1283,8 +1283,8 @@ struct has_op_inv_any< Glue, glue_times> > template struct has_nested_op_traits { - typedef char yes[1]; - typedef char no[2]; + using yes = char[1]; + using no = char[2]; template static yes& check(typename X::template traits*); template static no& check(...); @@ -1295,8 +1295,8 @@ struct has_nested_op_traits template struct has_nested_glue_traits { - typedef char yes[1]; - typedef char no[2]; + using yes = char[1]; + using no = char[2]; template static yes& check(typename X::template traits*); template static no& check(...); diff --git a/inst/include/armadillo_bits/unwrap.hpp b/inst/include/armadillo_bits/unwrap.hpp index b6be60bc..77ead018 100644 --- a/inst/include/armadillo_bits/unwrap.hpp +++ b/inst/include/armadillo_bits/unwrap.hpp @@ -21,7 +21,7 @@ // TODO: document the conditions and restrictions for the use of each unwrap variant: -// TODO: unwrap, unwrap_check, quasi_unwrap, partial_unwrap, partial_unwrap_check +// TODO: unwrap, unwrap_check, quasi_unwrap, partial_unwrap template @@ -2425,1038 +2425,6 @@ struct partial_unwrap< eOp, eop_neg> > -// - - - -template -struct partial_unwrap_check_default - { - typedef typename T1::elem_type eT; - typedef Mat stored_type; - - inline - partial_unwrap_check_default(const T1& A, const Mat&) - : M(A) - { - arma_debug_sigprint(); - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = false; - - const Mat M; - }; - - -template -struct partial_unwrap_check_fixed - { - typedef typename T1::elem_type eT; - typedef T1 stored_type; - - inline explicit - partial_unwrap_check_fixed(const T1& A, const Mat& B) - : M_local( (&A == &B) ? new T1(A) : nullptr ) - , M ( (&A == &B) ? (*M_local) : A ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check_fixed() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = false; - - const T1* M_local; - const T1& M; - }; - - - -template -struct partial_unwrap_check_redirect {}; - -template -struct partial_unwrap_check_redirect { typedef partial_unwrap_check_default result; }; - -template -struct partial_unwrap_check_redirect { typedef partial_unwrap_check_fixed result; }; - -template -struct partial_unwrap_check : public partial_unwrap_check_redirect::value>::result - { - typedef typename T1::elem_type eT; - - inline partial_unwrap_check(const T1& A, const Mat& B) - : partial_unwrap_check_redirect::value>::result(A, B) - { - } - }; - - - -template -struct partial_unwrap_check< Mat > - { - typedef Mat stored_type; - - inline - partial_unwrap_check(const Mat& A, const Mat& B) - : M_local ( (&A == &B) ? new Mat(A) : nullptr ) - , M ( (&A == &B) ? (*M_local) : A ) - { - arma_debug_sigprint(); - } - - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = false; - - // the order below is important - const Mat* M_local; - const Mat& M; - }; - - - -template -struct partial_unwrap_check< Row > - { - typedef Row stored_type; - - inline - partial_unwrap_check(const Row& A, const Mat& B) - : M_local ( (&A == &B) ? new Row(A) : nullptr ) - , M ( (&A == &B) ? (*M_local) : A ) - { - arma_debug_sigprint(); - } - - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = false; - - // the order below is important - const Row* M_local; - const Row& M; - }; - - - -template -struct partial_unwrap_check< Col > - { - typedef Col stored_type; - - inline - partial_unwrap_check(const Col& A, const Mat& B) - : M_local ( (&A == &B) ? new Col(A) : nullptr ) - , M ( (&A == &B) ? (*M_local) : A ) - { - arma_debug_sigprint(); - } - - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = false; - - // the order below is important - const Col* M_local; - const Col& M; - }; - - - -// NOTE: we can get away with this shortcut as the partial_unwrap_check class is only used by the glue_times class, -// NOTE: which relies on partial_unwrap_check to check for aliasing -template -struct partial_unwrap_check< subview_col > - { - typedef Col stored_type; - - inline - partial_unwrap_check(const subview_col& A, const Mat& B) - : M ( const_cast( A.colmem ), A.n_rows, (&(A.m) == &B), false ) - { - arma_debug_sigprint(); - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = false; - - const Col M; - }; - - - -template -struct partial_unwrap_check_htrans_default - { - typedef typename T1::elem_type eT; - typedef Mat stored_type; - - inline - partial_unwrap_check_htrans_default(const Op& A, const Mat&) - : M(A.m) - { - arma_debug_sigprint(); - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = true; - static constexpr bool do_times = false; - - const Mat M; - }; - - -template -struct partial_unwrap_check_htrans_fixed - { - typedef typename T1::elem_type eT; - typedef T1 stored_type; - - inline explicit - partial_unwrap_check_htrans_fixed(const Op& A, const Mat& B) - : M_local( (&(A.m) == &B) ? new T1(A.m) : nullptr ) - , M ( (&(A.m) == &B) ? (*M_local) : A.m ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check_htrans_fixed() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = true; - static constexpr bool do_times = false; - - const T1* M_local; - const T1& M; - }; - - - -template -struct partial_unwrap_check_htrans_redirect {}; - -template -struct partial_unwrap_check_htrans_redirect { typedef partial_unwrap_check_htrans_default result; }; - -template -struct partial_unwrap_check_htrans_redirect { typedef partial_unwrap_check_htrans_fixed result; }; - - -template -struct partial_unwrap_check< Op > : public partial_unwrap_check_htrans_redirect::value>::result - { - typedef typename T1::elem_type eT; - - inline partial_unwrap_check(const Op& A, const Mat& B) - : partial_unwrap_check_htrans_redirect::value>::result(A, B) - { - } - }; - - - -template -struct partial_unwrap_check< Op< Mat, op_htrans> > - { - typedef Mat stored_type; - - inline - partial_unwrap_check(const Op< Mat, op_htrans>& A, const Mat& B) - : M_local ( (&A.m == &B) ? new Mat(A.m) : nullptr ) - , M ( (&A.m == &B) ? (*M_local) : A.m ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = true; - static constexpr bool do_times = false; - - // the order below is important - const Mat* M_local; - const Mat& M; - }; - - - -template -struct partial_unwrap_check< Op< Row, op_htrans> > - { - typedef Row stored_type; - - inline - partial_unwrap_check(const Op< Row, op_htrans>& A, const Mat& B) - : M_local ( (&A.m == &B) ? new Row(A.m) : nullptr ) - , M ( (&A.m == &B) ? (*M_local) : A.m ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = true; - static constexpr bool do_times = false; - - // the order below is important - const Row* M_local; - const Row& M; - }; - - - -template -struct partial_unwrap_check< Op< Col, op_htrans> > - { - typedef Col stored_type; - - inline - partial_unwrap_check(const Op< Col, op_htrans>& A, const Mat& B) - : M_local ( (&A.m == &B) ? new Col(A.m) : nullptr ) - , M ( (&A.m == &B) ? (*M_local) : A.m ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = true; - static constexpr bool do_times = false; - - // the order below is important - const Col* M_local; - const Col& M; - }; - - - -// NOTE: we can get away with this shortcut as the partial_unwrap_check class is only used by the glue_times class, -// NOTE: which relies on partial_unwrap_check to check for aliasing -template -struct partial_unwrap_check< Op< subview_col, op_htrans> > - { - typedef Col stored_type; - - inline - partial_unwrap_check(const Op< subview_col, op_htrans>& A, const Mat& B) - : M ( const_cast( A.m.colmem ), A.m.n_rows, (&(A.m.m) == &B), false ) - { - arma_debug_sigprint(); - } - - constexpr eT get_val() const { return eT(1); } - - static constexpr bool do_trans = true; - static constexpr bool do_times = false; - - const Col M; - }; - - - -template -struct partial_unwrap_check_htrans2_default - { - typedef typename T1::elem_type eT; - typedef Mat stored_type; - - inline - partial_unwrap_check_htrans2_default(const Op& A, const Mat&) - : val(A.aux) - , M (A.m) - { - arma_debug_sigprint(); - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = true; - static constexpr bool do_times = true; - - const eT val; - const Mat M; - }; - - - -template -struct partial_unwrap_check_htrans2_fixed - { - typedef typename T1::elem_type eT; - typedef T1 stored_type; - - inline explicit - partial_unwrap_check_htrans2_fixed(const Op& A, const Mat& B) - : val (A.aux) - , M_local( (&(A.m) == &B) ? new T1(A.m) : nullptr ) - , M ( (&(A.m) == &B) ? (*M_local) : A.m ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check_htrans2_fixed() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = true; - static constexpr bool do_times = true; - - const eT val; - const T1* M_local; - const T1& M; - }; - - - -template -struct partial_unwrap_check_htrans2_redirect {}; - -template -struct partial_unwrap_check_htrans2_redirect { typedef partial_unwrap_check_htrans2_default result; }; - -template -struct partial_unwrap_check_htrans2_redirect { typedef partial_unwrap_check_htrans2_fixed result; }; - - -template -struct partial_unwrap_check< Op > : public partial_unwrap_check_htrans2_redirect::value>::result - { - typedef typename T1::elem_type eT; - - inline partial_unwrap_check(const Op& A, const Mat& B) - : partial_unwrap_check_htrans2_redirect::value>::result(A, B) - { - } - }; - - - -template -struct partial_unwrap_check< Op< Mat, op_htrans2> > - { - typedef Mat stored_type; - - inline - partial_unwrap_check(const Op< Mat, op_htrans2>& A, const Mat& B) - : val (A.aux) - , M_local ( (&A.m == &B) ? new Mat(A.m) : nullptr ) - , M ( (&A.m == &B) ? (*M_local) : A.m ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = true; - static constexpr bool do_times = true; - - // the order below is important - const eT val; - const Mat* M_local; - const Mat& M; - }; - - - -template -struct partial_unwrap_check< Op< Row, op_htrans2> > - { - typedef Row stored_type; - - inline - partial_unwrap_check(const Op< Row, op_htrans2>& A, const Mat& B) - : val (A.aux) - , M_local ( (&A.m == &B) ? new Row(A.m) : nullptr ) - , M ( (&A.m == &B) ? (*M_local) : A.m ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = true; - static constexpr bool do_times = true; - - // the order below is important - const eT val; - const Row* M_local; - const Row& M; - }; - - - -template -struct partial_unwrap_check< Op< Col, op_htrans2> > - { - typedef Col stored_type; - - inline - partial_unwrap_check(const Op< Col, op_htrans2>& A, const Mat& B) - : val (A.aux) - , M_local ( (&A.m == &B) ? new Col(A.m) : nullptr ) - , M ( (&A.m == &B) ? (*M_local) : A.m ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = true; - static constexpr bool do_times = true; - - // the order below is important - const eT val; - const Col* M_local; - const Col& M; - }; - - - -// NOTE: we can get away with this shortcut as the partial_unwrap_check class is only used by the glue_times class, -// NOTE: which relies on partial_unwrap_check to check for aliasing -template -struct partial_unwrap_check< Op< subview_col, op_htrans2> > - { - typedef Col stored_type; - - inline - partial_unwrap_check(const Op< subview_col, op_htrans2>& A, const Mat& B) - : val( A.aux ) - , M ( const_cast( A.m.colmem ), A.m.n_rows, (&(A.m.m) == &B), false ) - { - arma_debug_sigprint(); - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = true; - static constexpr bool do_times = true; - - const eT val; - const Col M; - }; - - - -template -struct partial_unwrap_check_scalar_times_default - { - typedef typename T1::elem_type eT; - typedef Mat stored_type; - - inline - partial_unwrap_check_scalar_times_default(const eOp& A, const Mat&) - : val(A.aux) - , M (A.P.Q) - { - arma_debug_sigprint(); - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const eT val; - const Mat M; - }; - - - -template -struct partial_unwrap_check_scalar_times_fixed - { - typedef typename T1::elem_type eT; - typedef T1 stored_type; - - inline explicit - partial_unwrap_check_scalar_times_fixed(const eOp& A, const Mat& B) - : val ( A.aux ) - , M_local( (&(A.P.Q) == &B) ? new T1(A.P.Q) : nullptr ) - , M ( (&(A.P.Q) == &B) ? (*M_local) : A.P.Q ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check_scalar_times_fixed() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const eT val; - const T1* M_local; - const T1& M; - }; - - - -template -struct partial_unwrap_check_scalar_times_redirect {}; - -template -struct partial_unwrap_check_scalar_times_redirect { typedef partial_unwrap_check_scalar_times_default result; }; - -template -struct partial_unwrap_check_scalar_times_redirect { typedef partial_unwrap_check_scalar_times_fixed result; }; - - -template -struct partial_unwrap_check< eOp > : public partial_unwrap_check_scalar_times_redirect::value>::result - { - typedef typename T1::elem_type eT; - - inline partial_unwrap_check(const eOp& A, const Mat& B) - : partial_unwrap_check_scalar_times_redirect::value>::result(A, B) - { - } - }; - - - -template -struct partial_unwrap_check< eOp, eop_scalar_times> > - { - typedef Mat stored_type; - - inline - partial_unwrap_check(const eOp,eop_scalar_times>& A, const Mat& B) - : val (A.aux) - , M_local( (&(A.P.Q) == &B) ? new Mat(A.P.Q) : nullptr ) - , M ( (&(A.P.Q) == &B) ? *M_local : A.P.Q ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const eT val; - const Mat* M_local; - const Mat& M; - }; - - - -template -struct partial_unwrap_check< eOp, eop_scalar_times> > - { - typedef Row stored_type; - - inline - partial_unwrap_check(const eOp,eop_scalar_times>& A, const Mat& B) - : val(A.aux) - , M_local( (&(A.P.Q) == &B) ? new Row(A.P.Q) : nullptr ) - , M ( (&(A.P.Q) == &B) ? *M_local : A.P.Q ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const eT val; - const Row* M_local; - const Row& M; - }; - - - -template -struct partial_unwrap_check< eOp, eop_scalar_times> > - { - typedef Col stored_type; - - inline - partial_unwrap_check(const eOp,eop_scalar_times>& A, const Mat& B) - : val ( A.aux ) - , M_local( (&(A.P.Q) == &B) ? new Col(A.P.Q) : nullptr ) - , M ( (&(A.P.Q) == &B) ? *M_local : A.P.Q ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const eT val; - const Col* M_local; - const Col& M; - }; - - - -// NOTE: we can get away with this shortcut as the partial_unwrap_check class is only used by the glue_times class, -// NOTE: which relies on partial_unwrap_check to check for aliasing -template -struct partial_unwrap_check< eOp, eop_scalar_times> > - { - typedef Col stored_type; - - inline - partial_unwrap_check(const eOp,eop_scalar_times>& A, const Mat& B) - : val( A.aux ) - , M ( const_cast( A.P.Q.colmem ), A.P.Q.n_rows, (&(A.P.Q.m) == &B), false ) - { - arma_debug_sigprint(); - } - - arma_inline eT get_val() const { return val; } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const eT val; - const Col M; - }; - - - -template -struct partial_unwrap_check_neg_default - { - typedef typename T1::elem_type eT; - typedef Mat stored_type; - - inline - partial_unwrap_check_neg_default(const eOp& A, const Mat&) - : M(A.P.Q) - { - arma_debug_sigprint(); - } - - constexpr eT get_val() const { return eT(-1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const Mat M; - }; - - - -template -struct partial_unwrap_check_neg_fixed - { - typedef typename T1::elem_type eT; - typedef T1 stored_type; - - inline explicit - partial_unwrap_check_neg_fixed(const eOp& A, const Mat& B) - : M_local( (&(A.P.Q) == &B) ? new T1(A.P.Q) : nullptr ) - , M ( (&(A.P.Q) == &B) ? (*M_local) : A.P.Q ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check_neg_fixed() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(-1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const T1* M_local; - const T1& M; - }; - - - -template -struct partial_unwrap_check_neg_redirect {}; - -template -struct partial_unwrap_check_neg_redirect { typedef partial_unwrap_check_neg_default result; }; - -template -struct partial_unwrap_check_neg_redirect { typedef partial_unwrap_check_neg_fixed result; }; - - -template -struct partial_unwrap_check< eOp > : public partial_unwrap_check_neg_redirect::value>::result - { - typedef typename T1::elem_type eT; - - inline partial_unwrap_check(const eOp& A, const Mat& B) - : partial_unwrap_check_neg_redirect::value>::result(A, B) - { - } - }; - - - -template -struct partial_unwrap_check< eOp, eop_neg> > - { - typedef Mat stored_type; - - inline - partial_unwrap_check(const eOp,eop_neg>& A, const Mat& B) - : M_local( (&(A.P.Q) == &B) ? new Mat(A.P.Q) : nullptr ) - , M ( (&(A.P.Q) == &B) ? *M_local : A.P.Q ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(-1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const Mat* M_local; - const Mat& M; - }; - - - -template -struct partial_unwrap_check< eOp, eop_neg> > - { - typedef Row stored_type; - - inline - partial_unwrap_check(const eOp,eop_neg>& A, const Mat& B) - : M_local( (&(A.P.Q) == &B) ? new Row(A.P.Q) : nullptr ) - , M ( (&(A.P.Q) == &B) ? *M_local : A.P.Q ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(-1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const Row* M_local; - const Row& M; - }; - - - -template -struct partial_unwrap_check< eOp, eop_neg> > - { - typedef Col stored_type; - - inline - partial_unwrap_check(const eOp,eop_neg>& A, const Mat& B) - : M_local( (&(A.P.Q) == &B) ? new Col(A.P.Q) : nullptr ) - , M ( (&(A.P.Q) == &B) ? *M_local : A.P.Q ) - { - arma_debug_sigprint(); - } - - inline - ~partial_unwrap_check() - { - arma_debug_sigprint(); - - if(M_local) { delete M_local; } - } - - constexpr eT get_val() const { return eT(-1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const Col* M_local; - const Col& M; - }; - - - -// NOTE: we can get away with this shortcut as the partial_unwrap_check class is only used by the glue_times class, -// NOTE: which relies on partial_unwrap_check to check for aliasing -template -struct partial_unwrap_check< eOp, eop_neg> > - { - typedef Col stored_type; - - inline - partial_unwrap_check(const eOp,eop_neg>& A, const Mat& B) - : M ( const_cast( A.P.Q.colmem ), A.P.Q.n_rows, (&(A.P.Q.m) == &B), false ) - { - arma_debug_sigprint(); - } - - constexpr eT get_val() const { return eT(-1); } - - static constexpr bool do_trans = false; - static constexpr bool do_times = true; - - const Col M; - }; - - - // // // From 451618dfae93faf3d321d4f313b056259f621353 Mon Sep 17 00:00:00 2001 From: Dirk Eddelbuettel Date: Tue, 11 Feb 2025 13:29:10 -0600 Subject: [PATCH 2/5] Test release 14.3.90-1 for reverse dependency checks --- ChangeLog | 4 ++++ DESCRIPTION | 4 ++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/ChangeLog b/ChangeLog index 229bdf57..1a625360 100644 --- a/ChangeLog +++ b/ChangeLog @@ -4,6 +4,10 @@ support additional DESCRIPTION fields * man/RcppArmadillo.package.skeleton.Rd: Document +2025-02-10 Dirk Eddelbuettel + + * inst/include/armadillo_bits/: Armadillo 14.3.90 as a RC for 14.4.* + 2025-02-05 Dirk Eddelbuettel * DESCRIPTION (Version, Date): RcppArmadillo 14.2.3-1 diff --git a/DESCRIPTION b/DESCRIPTION index 7a394c7d..a6e41c60 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: RcppArmadillo Type: Package Title: 'Rcpp' Integration for the 'Armadillo' Templated Linear Algebra Library -Version: 14.2.3-1 -Date: 2025-02-05 +Version: 14.3.90-1 +Date: 2025-02-10 Authors@R: c(person("Dirk", "Eddelbuettel", role = c("aut", "cre"), email = "edd@debian.org", comment = c(ORCID = "0000-0001-6419-907X")), person("Romain", "Francois", role = "aut", From 1dde93356ca38ec4ac001551b44a89d45821f1ef Mon Sep 17 00:00:00 2001 From: Dirk Eddelbuettel Date: Thu, 13 Feb 2025 08:31:00 -0600 Subject: [PATCH 3/5] Release 14.3.91-1 synced with 14.3.91 --- ChangeLog | 8 ++++ DESCRIPTION | 4 +- inst/include/armadillo_bits/Mat_meat.hpp | 16 ++++---- inst/include/armadillo_bits/arma_version.hpp | 4 +- inst/include/armadillo_bits/fn_accu.hpp | 6 +-- inst/include/armadillo_bits/fn_conv_to.hpp | 23 +++++++++++ inst/include/armadillo_bits/fn_elem.hpp | 42 ++++++++++++++++++++ 7 files changed, 88 insertions(+), 15 deletions(-) diff --git a/ChangeLog b/ChangeLog index 1a625360..eda83ff6 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,13 @@ +2025-02-13 Dirk Eddelbuettel + + * inst/include/armadillo_bits/: Sync again with updated sources for + Armadillo 14.3.91 as a RC for 14.4.* + 2025-02-11 Dirk Eddelbuettel + * inst/include/armadillo_bits/: Sync again with updated sources for + Armadillo 14.3.90 as a RC for 14.4.* + * R/RcppArmadillo.package.skeleton.R: Generalise helper function to support additional DESCRIPTION fields * man/RcppArmadillo.package.skeleton.Rd: Document diff --git a/DESCRIPTION b/DESCRIPTION index a6e41c60..62d4852b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: RcppArmadillo Type: Package Title: 'Rcpp' Integration for the 'Armadillo' Templated Linear Algebra Library -Version: 14.3.90-1 -Date: 2025-02-10 +Version: 14.3.91-1 +Date: 2025-02-13 Authors@R: c(person("Dirk", "Eddelbuettel", role = c("aut", "cre"), email = "edd@debian.org", comment = c(ORCID = "0000-0001-6419-907X")), person("Romain", "Francois", role = "aut", diff --git a/inst/include/armadillo_bits/Mat_meat.hpp b/inst/include/armadillo_bits/Mat_meat.hpp index e18b1d74..596b39b0 100644 --- a/inst/include/armadillo_bits/Mat_meat.hpp +++ b/inst/include/armadillo_bits/Mat_meat.hpp @@ -474,9 +474,9 @@ Mat::Mat(const char* text) init( std::string(text) ); } - - - + + + //! create the matrix from a textual description template inline @@ -489,8 +489,8 @@ Mat::operator=(const char* text) return *this; } - - + + //! create the matrix from a textual description template @@ -508,9 +508,9 @@ Mat::Mat(const std::string& text) init(text); } - - - + + + //! create the matrix from a textual description template inline diff --git a/inst/include/armadillo_bits/arma_version.hpp b/inst/include/armadillo_bits/arma_version.hpp index c3f06874..2e82ba7d 100644 --- a/inst/include/armadillo_bits/arma_version.hpp +++ b/inst/include/armadillo_bits/arma_version.hpp @@ -23,8 +23,8 @@ #define ARMA_VERSION_MAJOR 14 #define ARMA_VERSION_MINOR 3 -#define ARMA_VERSION_PATCH 90 -#define ARMA_VERSION_NAME "unstable" +#define ARMA_VERSION_PATCH 91 +#define ARMA_VERSION_NAME "14.4-RC1" diff --git a/inst/include/armadillo_bits/fn_accu.hpp b/inst/include/armadillo_bits/fn_accu.hpp index 57dbb36c..8990bd44 100644 --- a/inst/include/armadillo_bits/fn_accu.hpp +++ b/inst/include/armadillo_bits/fn_accu.hpp @@ -288,7 +288,7 @@ accu(const eOp& expr) if((is_Mat::value) || (is_subview_col::value)) { - const quasi_unwrap U(expr.P.Q); + const quasi_unwrap U(expr.P.Q); const eT* X_mem = U.M.memptr(); @@ -987,8 +987,8 @@ accu(const eGlueCube& expr) typedef eGlueCube expr_type; - typedef typename expr_type::proxy_type::stored_type P1_stored_type; - typedef typename expr_type::proxy_type::stored_type P2_stored_type; + typedef typename expr_type::proxy1_type::stored_type P1_stored_type; + typedef typename expr_type::proxy2_type::stored_type P2_stored_type; if(is_Cube::value && is_Cube::value) { diff --git a/inst/include/armadillo_bits/fn_conv_to.hpp b/inst/include/armadillo_bits/fn_conv_to.hpp index 3639fdd3..15985121 100644 --- a/inst/include/armadillo_bits/fn_conv_to.hpp +++ b/inst/include/armadillo_bits/fn_conv_to.hpp @@ -28,6 +28,9 @@ class conv_to { public: + template + arma_frown("use as_scalar() instead") inline static out_eT from(const in_eT& in, const typename arma_scalar_only::result* junk = nullptr); + template arma_frown("use as_scalar() instead") inline static out_eT from(const Base& in, const typename arma_not_cx::result* junk = nullptr); @@ -43,6 +46,26 @@ class conv_to +template +template +arma_warn_unused +inline +out_eT +conv_to::from(const in_eT& in, const typename arma_scalar_only::result* junk) + { + arma_debug_sigprint(); + arma_ignore(junk); + + arma_type_check(( is_supported_elem_type::value == false )); + + // NOTE: this is meant only as a workaround for old user code; + // NOTE: it doesn't handle conversions from complex to real + + return out_eT(in); + } + + + template template arma_warn_unused diff --git a/inst/include/armadillo_bits/fn_elem.hpp b/inst/include/armadillo_bits/fn_elem.hpp index ad4b0d11..934c921c 100644 --- a/inst/include/armadillo_bits/fn_elem.hpp +++ b/inst/include/armadillo_bits/fn_elem.hpp @@ -603,6 +603,20 @@ sqrt(const T1& A) +// workaround for old user code +template +arma_frown("use std::sqrt(arma::as_scalar(X)) instead") +inline +typename arma_scalar_only::result +sqrt(const eT& A) + { + arma_debug_sigprint(); + + return eop_aux::sqrt(A); + } + + + template arma_warn_unused arma_inline @@ -767,6 +781,20 @@ pow(const T1& A, const typename T1::elem_type exponent) +// workaround for old user code +template +arma_frown("use std::pow(arma::as_scalar(X),Y) instead") +inline +typename arma_scalar_only::result +pow(const eT& A, const eT exponent) + { + arma_debug_sigprint(); + + return eop_aux::pow(A, exponent); + } + + + template arma_warn_unused arma_inline @@ -797,6 +825,20 @@ pow(const T1& A, const typename T1::elem_type::value_type exponent) +// workaround for old user code +template +arma_frown("use std::pow(arma::as_scalar(X),Y) instead") +inline +typename enable_if2< is_cx::yes, eT >::result +pow(const eT& A, const typename eT::value_type exponent) + { + arma_debug_sigprint(); + + return eop_aux::pow(A, eT(exponent)); + } + + + template arma_warn_unused arma_inline From 0d7537eca39588886afd6e8de3650d6b6cfd2358 Mon Sep 17 00:00:00 2001 From: Dirk Eddelbuettel Date: Sat, 15 Feb 2025 10:22:56 -0600 Subject: [PATCH 4/5] RcppArmadillo 14.3.92-1 with RC-2 of Armadillo 14.4 --- ChangeLog | 5 ++ DESCRIPTION | 4 +- inst/include/armadillo_bits/arma_forward.hpp | 1 + inst/include/armadillo_bits/arma_version.hpp | 4 +- inst/include/armadillo_bits/fn_accu.hpp | 18 +++- inst/include/armadillo_bits/fn_dot.hpp | 2 +- inst/include/armadillo_bits/fn_elem.hpp | 58 +++++-------- .../armadillo_bits/op_sp_sum_bones.hpp | 3 + .../include/armadillo_bits/op_sp_sum_meat.hpp | 82 ++++++++++++++++--- 9 files changed, 122 insertions(+), 55 deletions(-) diff --git a/ChangeLog b/ChangeLog index eda83ff6..19ac41a2 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,8 @@ +2025-02-15 Dirk Eddelbuettel + + * inst/include/armadillo_bits/: Sync again with updated sources for + Armadillo 14.3.92 as an RC for 14.4.* + 2025-02-13 Dirk Eddelbuettel * inst/include/armadillo_bits/: Sync again with updated sources for diff --git a/DESCRIPTION b/DESCRIPTION index 62d4852b..7b9f9693 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: RcppArmadillo Type: Package Title: 'Rcpp' Integration for the 'Armadillo' Templated Linear Algebra Library -Version: 14.3.91-1 -Date: 2025-02-13 +Version: 14.3.92-1 +Date: 2025-02-15 Authors@R: c(person("Dirk", "Eddelbuettel", role = c("aut", "cre"), email = "edd@debian.org", comment = c(ORCID = "0000-0001-6419-907X")), person("Romain", "Francois", role = "aut", diff --git a/inst/include/armadillo_bits/arma_forward.hpp b/inst/include/armadillo_bits/arma_forward.hpp index b56d4e4d..0eb857ee 100644 --- a/inst/include/armadillo_bits/arma_forward.hpp +++ b/inst/include/armadillo_bits/arma_forward.hpp @@ -130,6 +130,7 @@ class spop_strans; class spop_htrans; class spop_vectorise_row; class spop_vectorise_col; +class spop_square; class spop_rel_lt_pre; class spop_rel_lt_post; diff --git a/inst/include/armadillo_bits/arma_version.hpp b/inst/include/armadillo_bits/arma_version.hpp index 2e82ba7d..8d306868 100644 --- a/inst/include/armadillo_bits/arma_version.hpp +++ b/inst/include/armadillo_bits/arma_version.hpp @@ -23,8 +23,8 @@ #define ARMA_VERSION_MAJOR 14 #define ARMA_VERSION_MINOR 3 -#define ARMA_VERSION_PATCH 91 -#define ARMA_VERSION_NAME "14.4-RC1" +#define ARMA_VERSION_PATCH 92 +#define ARMA_VERSION_NAME "14.4-RC2" diff --git a/inst/include/armadillo_bits/fn_accu.hpp b/inst/include/armadillo_bits/fn_accu.hpp index 8990bd44..b7287dcf 100644 --- a/inst/include/armadillo_bits/fn_accu.hpp +++ b/inst/include/armadillo_bits/fn_accu.hpp @@ -1119,6 +1119,22 @@ accu(const SpGlue& expr) const SpProxy px(expr.A); const SpProxy py(expr.B); + arma_conform_assert_same_size(px.get_n_rows(), px.get_n_cols(), py.get_n_rows(), py.get_n_cols(), "element-wise multiplication"); + + typedef typename SpProxy::stored_type px_Q_type; + typedef typename SpProxy::stored_type py_Q_type; + + if(is_SpMat::value && is_SpMat::value) + { + const unwrap_spmat UX(px.Q); + const unwrap_spmat UY(py.Q); + + const SpMat& X = UX.M; + const SpMat& Y = UY.M; + + if(&X == &Y) { return op_dot::direct_dot(X.n_nonzero, X.values, X.values); } + } + typename SpProxy::const_iterator_type x_it = px.begin(); typename SpProxy::const_iterator_type x_it_end = px.end(); @@ -1195,7 +1211,7 @@ accu(const SpOp& expr) eT val = eT(0); - for(uword i=0; i < N; ++i) { const eT tmp = (*it); ++it; val += tmp*tmp; } + for(uword i=0; i < N; ++i) { const eT tmp = (*it); val += (tmp*tmp); ++it; } return val; } diff --git a/inst/include/armadillo_bits/fn_dot.hpp b/inst/include/armadillo_bits/fn_dot.hpp index a9c26ca5..0a9e2c1e 100644 --- a/inst/include/armadillo_bits/fn_dot.hpp +++ b/inst/include/armadillo_bits/fn_dot.hpp @@ -256,7 +256,7 @@ dot if( &A == &B ) { // We can do it directly! - return op_dot::direct_dot_arma(A.n_nonzero, A.values, A.values); + return op_dot::direct_dot(A.n_nonzero, A.values, A.values); } else { diff --git a/inst/include/armadillo_bits/fn_elem.hpp b/inst/include/armadillo_bits/fn_elem.hpp index 934c921c..51c9e8a1 100644 --- a/inst/include/armadillo_bits/fn_elem.hpp +++ b/inst/include/armadillo_bits/fn_elem.hpp @@ -603,20 +603,6 @@ sqrt(const T1& A) -// workaround for old user code -template -arma_frown("use std::sqrt(arma::as_scalar(X)) instead") -inline -typename arma_scalar_only::result -sqrt(const eT& A) - { - arma_debug_sigprint(); - - return eop_aux::sqrt(A); - } - - - template arma_warn_unused arma_inline @@ -781,17 +767,17 @@ pow(const T1& A, const typename T1::elem_type exponent) -// workaround for old user code -template -arma_frown("use std::pow(arma::as_scalar(X),Y) instead") -inline -typename arma_scalar_only::result -pow(const eT& A, const eT exponent) - { - arma_debug_sigprint(); - - return eop_aux::pow(A, exponent); - } +// // workaround for old user code +// template +// arma_frown("use std::pow(arma::as_scalar(X),Y) instead") +// inline +// typename arma_scalar_only::result +// pow(const eT& A, const eT exponent) +// { +// arma_debug_sigprint(); +// +// return eop_aux::pow(A, exponent); +// } @@ -825,17 +811,17 @@ pow(const T1& A, const typename T1::elem_type::value_type exponent) -// workaround for old user code -template -arma_frown("use std::pow(arma::as_scalar(X),Y) instead") -inline -typename enable_if2< is_cx::yes, eT >::result -pow(const eT& A, const typename eT::value_type exponent) - { - arma_debug_sigprint(); - - return eop_aux::pow(A, eT(exponent)); - } +// // workaround for old user code +// template +// arma_frown("use std::pow(arma::as_scalar(X),Y) instead") +// inline +// typename enable_if2< is_cx::yes, eT >::result +// pow(const eT& A, const typename eT::value_type exponent) +// { +// arma_debug_sigprint(); +// +// return eop_aux::pow(A, eT(exponent)); +// } diff --git a/inst/include/armadillo_bits/op_sp_sum_bones.hpp b/inst/include/armadillo_bits/op_sp_sum_bones.hpp index 8c9c67a0..63db4566 100644 --- a/inst/include/armadillo_bits/op_sp_sum_bones.hpp +++ b/inst/include/armadillo_bits/op_sp_sum_bones.hpp @@ -27,6 +27,9 @@ class op_sp_sum template inline static void apply(Mat& out, const mtSpReduceOp& in); + + template + inline static void apply(Mat& out, const mtSpReduceOp, op_sp_sum>& in); }; diff --git a/inst/include/armadillo_bits/op_sp_sum_meat.hpp b/inst/include/armadillo_bits/op_sp_sum_meat.hpp index 8e9a22da..539242d9 100644 --- a/inst/include/armadillo_bits/op_sp_sum_meat.hpp +++ b/inst/include/armadillo_bits/op_sp_sum_meat.hpp @@ -54,21 +54,19 @@ op_sp_sum::apply(Mat& out, const mtSpReduceOp& out, const mtSpReduceOp +inline +void +op_sp_sum::apply(Mat& out, const mtSpReduceOp, op_sp_sum>& in) + { + arma_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const uword dim = in.aux_uword_a; + + arma_conform_check( (dim > 1), "sum(): parameter 'dim' must be 0 or 1" ); + + const SpProxy p(in.m.m); + + const uword p_n_rows = p.get_n_rows(); + const uword p_n_cols = p.get_n_cols(); + + if(dim == 0) { out.zeros(1, p_n_cols); } + if(dim == 1) { out.zeros(p_n_rows, 1); } + + if(p.get_n_nonzero() == 0) { return; } + + eT* out_mem = out.memptr(); + + if(dim == 0) // find the sum of squares in each column + { + if(SpProxy::use_iterator) { - out_mem[it.row()] += (*it); - ++it; + typename SpProxy::const_iterator_type it = p.begin(); + + const uword N = p.get_n_nonzero(); + + for(uword i=0; i < N; ++i) { const eT val = (*it); out_mem[it.col()] += (val*val); ++it; } } + else + { + const eT* values = p.get_values(); + const uword* colptrs = p.get_col_ptrs(); + + for(uword col = 0; col < p_n_cols; ++col) + { + const eT* coldata = &(values[ colptrs[col] ]); + const uword N = colptrs[col + 1] - colptrs[col]; + + out_mem[col] = op_dot::direct_dot(N, coldata, coldata); + } + } + } + else + if(dim == 1) // find the sum of squares in each row + { + typename SpProxy::const_iterator_type it = p.begin(); + + const uword N = p.get_n_nonzero(); + + for(uword i=0; i < N; ++i) { const eT val = (*it); out_mem[it.row()] += (val*val); ++it; } } } From 5db4dcbef1af04cfc03866b597954e716b5aaeb4 Mon Sep 17 00:00:00 2001 From: Dirk Eddelbuettel Date: Mon, 17 Feb 2025 07:07:50 -0600 Subject: [PATCH 5/5] Sync with Armadillo 14.4.0 --- inst/include/armadillo_bits/arma_version.hpp | 6 ++--- inst/include/armadillo_bits/fn_elem.hpp | 28 -------------------- 2 files changed, 3 insertions(+), 31 deletions(-) diff --git a/inst/include/armadillo_bits/arma_version.hpp b/inst/include/armadillo_bits/arma_version.hpp index 8d306868..c9c50d49 100644 --- a/inst/include/armadillo_bits/arma_version.hpp +++ b/inst/include/armadillo_bits/arma_version.hpp @@ -22,9 +22,9 @@ #define ARMA_VERSION_MAJOR 14 -#define ARMA_VERSION_MINOR 3 -#define ARMA_VERSION_PATCH 92 -#define ARMA_VERSION_NAME "14.4-RC2" +#define ARMA_VERSION_MINOR 4 +#define ARMA_VERSION_PATCH 0 +#define ARMA_VERSION_NAME "Filtered Espresso" diff --git a/inst/include/armadillo_bits/fn_elem.hpp b/inst/include/armadillo_bits/fn_elem.hpp index 51c9e8a1..ad4b0d11 100644 --- a/inst/include/armadillo_bits/fn_elem.hpp +++ b/inst/include/armadillo_bits/fn_elem.hpp @@ -767,20 +767,6 @@ pow(const T1& A, const typename T1::elem_type exponent) -// // workaround for old user code -// template -// arma_frown("use std::pow(arma::as_scalar(X),Y) instead") -// inline -// typename arma_scalar_only::result -// pow(const eT& A, const eT exponent) -// { -// arma_debug_sigprint(); -// -// return eop_aux::pow(A, exponent); -// } - - - template arma_warn_unused arma_inline @@ -811,20 +797,6 @@ pow(const T1& A, const typename T1::elem_type::value_type exponent) -// // workaround for old user code -// template -// arma_frown("use std::pow(arma::as_scalar(X),Y) instead") -// inline -// typename enable_if2< is_cx::yes, eT >::result -// pow(const eT& A, const typename eT::value_type exponent) -// { -// arma_debug_sigprint(); -// -// return eop_aux::pow(A, eT(exponent)); -// } - - - template arma_warn_unused arma_inline