diff --git a/ChangeLog b/ChangeLog index 9f855357..325890fc 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,10 +1,14 @@ +2026-03-16 Dirk Eddelbuettel + + * inst/include/current/: Sync with Armadillo 15.2.4 + 2026-01-07 Dirk Eddelbuettel * .github/workflows/ci.yaml: Switch to actions/checkout@v6 2026-01-04 Dirk Eddelbuettel - * NAMESPACE: Also import stats:na.omit and utils::packageVersion for + * NAMESPACE: Also import stats::na.omit and utils::packageVersion for use in init.R 2026-01-03 Dirk Eddelbuettel diff --git a/DESCRIPTION b/DESCRIPTION index 97f6e613..126768e2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: RcppArmadillo Type: Package Title: 'Rcpp' Integration for the 'Armadillo' Templated Linear Algebra Library -Version: 15.2.3-1.1 -Date: 2026-01-03 +Version: 15.2.4-0 +Date: 2026-03-16 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/R/init.R b/R/init.R index 58d952c7..3f9be66f 100644 --- a/R/init.R +++ b/R/init.R @@ -32,7 +32,8 @@ .onAttach <- function(libname, pkgname) { if (interactive()) { - packageStartupMessage("RcppArmadillo ", packageVersion("RcppArmadillo"), + packageStartupMessage("RcppArmadillo ", + packageDescription("RcppArmadillo", fields="Version"), " using ", .pkgenv[["omp_threads"]], " cores. See ", "'help(\"RcppArmadillo-package\")' for details.") } diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index ca6fbbfc..28dc5848 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -3,9 +3,14 @@ \newcommand{\ghpr}{\href{https://github.com/RcppCore/RcppArmadillo/pull/#1}{##1}} \newcommand{\ghit}{\href{https://github.com/RcppCore/RcppArmadillo/issues/#1}{##1}} -\section{Changes in RcppArmadillo version 15.2.3-1.1 (2026-01-03)}{ +\section{Changes in RcppArmadillo version 15.2.4-0 (2026-03-16)}{ \itemize{ - \item Refined OpenMP setup (Dirk in \ghpr{500} + \item Upgraded to Armadillo release 15.2.4 (Medium Roast Deluxe) + \itemize{ + \item Workarounds for bugs in GCC and Clang sanitisers (ASAN false positives) + \item Faster handling of blank sparse matrices + } + \item Refined OpenMP setup (Dirk in \ghpr{500}) } } diff --git a/inst/include/current/armadillo_bits/Col_meat.hpp b/inst/include/current/armadillo_bits/Col_meat.hpp index ba41bc41..91d93eb3 100644 --- a/inst/include/current/armadillo_bits/Col_meat.hpp +++ b/inst/include/current/armadillo_bits/Col_meat.hpp @@ -1344,6 +1344,10 @@ Col::fixed::fixed(const fill::fill_class&) if(is_same_type::yes) { Mat::eye(); } if(is_same_type::yes) { Mat::randu(); } if(is_same_type::yes) { Mat::randn(); } + + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::nan () ); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::pos_inf() ); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::neg_inf() ); } } diff --git a/inst/include/current/armadillo_bits/Cube_bones.hpp b/inst/include/current/armadillo_bits/Cube_bones.hpp index cb662f90..3c97d4ae 100644 --- a/inst/include/current/armadillo_bits/Cube_bones.hpp +++ b/inst/include/current/armadillo_bits/Cube_bones.hpp @@ -360,6 +360,9 @@ class Cube : public BaseCube< eT, Cube > inline Cube& fill(const eT val); + template + inline Cube& fill(const fill::fill_class& f); + inline Cube& zeros(); inline Cube& zeros(const uword new_n_rows, const uword new_n_cols, const uword new_n_slices); inline Cube& zeros(const SizeCube& s); diff --git a/inst/include/current/armadillo_bits/Cube_meat.hpp b/inst/include/current/armadillo_bits/Cube_meat.hpp index 62e2c562..2571abc1 100644 --- a/inst/include/current/armadillo_bits/Cube_meat.hpp +++ b/inst/include/current/armadillo_bits/Cube_meat.hpp @@ -52,7 +52,7 @@ Cube::Cube() , n_elem(0) , n_alloc(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); } @@ -172,7 +172,7 @@ Cube::Cube(const SizeCube& s, const arma_initmode_indicator&) template template inline -Cube::Cube(const uword in_n_rows, const uword in_n_cols, const uword in_n_slices, const fill::fill_class&) +Cube::Cube(const uword in_n_rows, const uword in_n_cols, const uword in_n_slices, const fill::fill_class& f) : n_rows(in_n_rows) , n_cols(in_n_cols) , n_elem_slice(in_n_rows*in_n_cols) @@ -186,12 +186,7 @@ Cube::Cube(const uword in_n_rows, const uword in_n_cols, const uword in_n_sl init_cold(); - if(is_same_type::yes) { (*this).zeros(); } - if(is_same_type::yes) { (*this).ones(); } - if(is_same_type::yes) { (*this).randu(); } - if(is_same_type::yes) { (*this).randn(); } - - arma_static_check( (is_same_type::yes), "Cube::Cube(): unsupported fill type" ); + (*this).fill(f); } @@ -199,7 +194,7 @@ Cube::Cube(const uword in_n_rows, const uword in_n_cols, const uword in_n_sl template template inline -Cube::Cube(const SizeCube& s, const fill::fill_class&) +Cube::Cube(const SizeCube& s, const fill::fill_class& f) : n_rows(s.n_rows) , n_cols(s.n_cols) , n_elem_slice(s.n_rows*s.n_cols) @@ -213,12 +208,7 @@ Cube::Cube(const SizeCube& s, const fill::fill_class&) init_cold(); - if(is_same_type::yes) { (*this).zeros(); } - if(is_same_type::yes) { (*this).ones(); } - if(is_same_type::yes) { (*this).randu(); } - if(is_same_type::yes) { (*this).randn(); } - - arma_static_check( (is_same_type::yes), "Cube::Cube(): unsupported fill type" ); + (*this).fill(f); } @@ -934,7 +924,7 @@ Cube::Cube , n_elem(0) , n_alloc(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -1064,7 +1054,7 @@ Cube::Cube(const subview_cube_slices& X) , n_elem(0) , n_alloc(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2751,7 +2741,7 @@ Cube::Cube(const OpCube& X) , n_elem(0) , n_alloc(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -3048,7 +3038,7 @@ Cube::Cube(const mtOpCube& X) , n_elem(0) , n_alloc(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -3144,7 +3134,7 @@ Cube::Cube(const GlueCube& X) , n_elem(0) , n_alloc(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -3401,7 +3391,7 @@ Cube::Cube(const mtGlueCube& X) , n_elem(0) , n_alloc(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -3887,7 +3877,7 @@ arma_inline eT* Cube::slice_memptr(const uword uslice) { - return const_cast( &mem[ uslice*n_elem_slice ] ); + return access::rwp( mem + (uslice*n_elem_slice) ); } @@ -3898,7 +3888,7 @@ arma_inline const eT* Cube::slice_memptr(const uword uslice) const { - return &mem[ uslice*n_elem_slice ]; + return mem + (uslice*n_elem_slice); } @@ -3909,7 +3899,7 @@ arma_inline eT* Cube::slice_colptr(const uword uslice, const uword col) { - return const_cast( &mem[ uslice*n_elem_slice + col*n_rows] ); + return access::rwp( mem + (uslice*n_elem_slice + col*n_rows) ); } @@ -3920,7 +3910,7 @@ arma_inline const eT* Cube::slice_colptr(const uword uslice, const uword col) const { - return &mem[ uslice*n_elem_slice + col*n_rows ]; + return mem + (uslice*n_elem_slice + col*n_rows); } @@ -4227,6 +4217,30 @@ Cube::fill(const eT val) +template +template +inline +Cube& +Cube::fill(const fill::fill_class&) + { + arma_debug_sigprint(); + + arma_static_check( (is_same_type::yes), "Cube::fill(): unsupported fill type" ); + + if(is_same_type::yes) { (*this).zeros(); } + if(is_same_type::yes) { (*this).ones(); } + if(is_same_type::yes) { (*this).randu(); } + if(is_same_type::yes) { (*this).randn(); } + + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::nan () ); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::pos_inf() ); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::neg_inf() ); } + + return *this; + } + + + template inline Cube& @@ -5458,18 +5472,13 @@ template template template inline -Cube::fixed::fixed(const fill::fill_class&) +Cube::fixed::fixed(const fill::fill_class& f) { arma_debug_sigprint_this(this); mem_setup(); - if(is_same_type::yes) { Cube::zeros(); } - if(is_same_type::yes) { Cube::ones(); } - if(is_same_type::yes) { Cube::randu(); } - if(is_same_type::yes) { Cube::randn(); } - - arma_static_check( (is_same_type::yes), "Cube::fixed::fixed(): unsupported fill type" ); + (*this).fill(f); } diff --git a/inst/include/current/armadillo_bits/Mat_meat.hpp b/inst/include/current/armadillo_bits/Mat_meat.hpp index 4e4fa3b6..cd3ee13e 100644 --- a/inst/include/current/armadillo_bits/Mat_meat.hpp +++ b/inst/include/current/armadillo_bits/Mat_meat.hpp @@ -255,7 +255,7 @@ Mat::Mat(const arma_vec_indicator&, const uhword in_vec_state) , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); } @@ -456,7 +456,7 @@ Mat::Mat(const char* text) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -490,7 +490,7 @@ Mat::Mat(const std::string& text) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -678,7 +678,7 @@ Mat::Mat(const std::initializer_list& list) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -710,7 +710,7 @@ Mat::Mat(const std::initializer_list< std::initializer_list >& list) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -877,7 +877,7 @@ Mat::Mat(const Mat& in_mat, const arma_vec_indicator&, const uhword in_v , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint(arma_str::format("this: %x; in_mat: %x") % this % &in_mat); @@ -1512,7 +1512,7 @@ Mat::Mat(const BaseCube& X, const arma_vec_indicator&, const uhword i , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -1531,7 +1531,7 @@ Mat::Mat(const BaseCube& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -1978,7 +1978,7 @@ Mat::Mat(const Base::pod_type,T1>& A, const Base::Mat(const Base::pod_type,T1>& A, const Base::Mat(const subview& X, const arma_vec_indicator&, const uhword in_ve , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2266,7 +2266,7 @@ Mat::Mat(const subview_cube& x, const arma_vec_indicator&, const uhword , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2285,7 +2285,7 @@ Mat::Mat(const subview_cube& x) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2395,7 +2395,7 @@ Mat::Mat(const diagview& X, const arma_vec_indicator&, const uhword in_v , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2540,7 +2540,7 @@ Mat::Mat(const subview_elem1& X, const arma_vec_indicator&, const uhw , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2559,7 +2559,7 @@ Mat::Mat(const subview_elem1& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2668,7 +2668,7 @@ Mat::Mat(const subview_elem2& X, const arma_vec_indicator&, const , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2687,7 +2687,7 @@ Mat::Mat(const subview_elem2& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2796,7 +2796,7 @@ Mat::Mat(const SpBase& m, const arma_vec_indicator&, const uhword in , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2815,7 +2815,7 @@ Mat::Mat(const SpBase& m) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -2876,6 +2876,8 @@ Mat::operator+=(const SpBase& m) arma_conform_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols(), "addition"); + if(p.get_n_nonzero() == 0) { return *this; } + typename SpProxy::const_iterator_type it = p.begin(); typename SpProxy::const_iterator_type it_end = p.end(); @@ -2898,6 +2900,8 @@ Mat::operator-=(const SpBase& m) arma_conform_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols(), "subtraction"); + if(p.get_n_nonzero() == 0) { return *this; } + typename SpProxy::const_iterator_type it = p.begin(); typename SpProxy::const_iterator_type it_end = p.end(); @@ -2916,9 +2920,9 @@ Mat::operator*=(const SpBase& m) { arma_debug_sigprint(); - Mat z = (*this) * m.get_ref(); + Mat tmp = (*this) * m.get_ref(); - steal_mem(z); + steal_mem(tmp); return *this; } @@ -2999,7 +3003,7 @@ Mat::Mat(const SpSubview& X, const arma_vec_indicator&, const uhword in_ , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -3017,7 +3021,7 @@ Mat::Mat(const SpSubview& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -3200,7 +3204,7 @@ Mat::Mat(const spdiagview& X, const arma_vec_indicator&, const uhword in , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -5128,7 +5132,7 @@ Mat::Mat(const Gen& X, const arma_vec_indicator&, const uhword , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -5278,7 +5282,7 @@ Mat::Mat(const Op& X, const arma_vec_indicator&, const uhword i , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -5300,7 +5304,7 @@ Mat::Mat(const Op& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -5429,7 +5433,7 @@ Mat::Mat(const eOp& X, const arma_vec_indicator&, const uhword , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -5657,7 +5661,7 @@ Mat::Mat(const mtOp& X, const arma_vec_indicator&, const uh , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -5676,7 +5680,7 @@ Mat::Mat(const mtOp& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -5785,7 +5789,7 @@ Mat::Mat(const CubeToMatOp& X, const arma_vec_indicator&, const , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -5806,7 +5810,7 @@ Mat::Mat(const CubeToMatOp& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -5929,7 +5933,7 @@ Mat::Mat(const SpToDOp& X, const arma_vec_indicator&, const uhw , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -5950,7 +5954,7 @@ Mat::Mat(const SpToDOp& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -6079,7 +6083,7 @@ Mat::Mat(const mtSpReduceOp& X, const arma_vec_indicator&, , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -6098,7 +6102,7 @@ Mat::Mat(const mtSpReduceOp& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -6207,7 +6211,7 @@ Mat::Mat(const Glue& X, const arma_vec_indicator&, const , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -6230,7 +6234,7 @@ Mat::Mat(const Glue& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -6396,7 +6400,7 @@ Mat::Mat(const eGlue& X, const arma_vec_indicator&, cons , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -6603,7 +6607,7 @@ Mat::Mat(const mtGlue& X, const arma_vec_indicator&, , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -6622,7 +6626,7 @@ Mat::Mat(const mtGlue& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -6733,7 +6737,7 @@ Mat::Mat(const SpToDGlue& X, const arma_vec_indicator&, c , n_alloc(0) , vec_state(in_vec_state) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -6755,7 +6759,7 @@ Mat::Mat(const SpToDGlue& X) , n_alloc(0) , vec_state(0) , mem_state(0) - , mem() + , mem(nullptr) { arma_debug_sigprint_this(this); @@ -7475,7 +7479,7 @@ arma_inline eT* Mat::colptr(const uword in_col) { - return & access::rw(mem[in_col*n_rows]); + return access::rwp( mem + (in_col*n_rows) ); } @@ -7486,7 +7490,7 @@ arma_inline const eT* Mat::colptr(const uword in_col) const { - return & mem[in_col*n_rows]; + return mem + (in_col*n_rows); } @@ -7926,6 +7930,10 @@ Mat::fill(const fill::fill_class&) if(is_same_type::yes) { (*this).randu(); } if(is_same_type::yes) { (*this).randn(); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::nan () ); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::pos_inf() ); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::neg_inf() ); } + return *this; } @@ -10099,6 +10107,10 @@ Mat::fixed::fixed(const fill::fill_class::yes) { Mat::eye(); } if(is_same_type::yes) { Mat::randu(); } if(is_same_type::yes) { Mat::randn(); } + + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::nan () ); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::pos_inf() ); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::neg_inf() ); } } @@ -10533,7 +10545,7 @@ Mat::fixed::colptr(const uword in_col) { eT* mem_actual = (use_extra) ? mem_local_extra : mem_local; - return & access::rw(mem_actual[in_col*fixed_n_rows]); + return access::rwp( mem_actual + (in_col*fixed_n_rows) ); } @@ -10546,7 +10558,7 @@ Mat::fixed::colptr(const uword in_col) const { const eT* mem_actual = (use_extra) ? mem_local_extra : mem_local; - return & mem_actual[in_col*fixed_n_rows]; + return mem_actual + (in_col*fixed_n_rows); } diff --git a/inst/include/current/armadillo_bits/Row_meat.hpp b/inst/include/current/armadillo_bits/Row_meat.hpp index 987a0f48..012423a1 100644 --- a/inst/include/current/armadillo_bits/Row_meat.hpp +++ b/inst/include/current/armadillo_bits/Row_meat.hpp @@ -1351,6 +1351,10 @@ Row::fixed::fixed(const fill::fill_class&) if(is_same_type::yes) { Mat::eye(); } if(is_same_type::yes) { Mat::randu(); } if(is_same_type::yes) { Mat::randn(); } + + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::nan () ); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::pos_inf() ); } + if(is_same_type::yes) { (*this).fill( priv::Datum_helper::neg_inf() ); } } diff --git a/inst/include/current/armadillo_bits/SpMat_meat.hpp b/inst/include/current/armadillo_bits/SpMat_meat.hpp index 266f6d95..c17fed4d 100644 --- a/inst/include/current/armadillo_bits/SpMat_meat.hpp +++ b/inst/include/current/armadillo_bits/SpMat_meat.hpp @@ -688,15 +688,22 @@ SpMat::operator=(const SpMat& x) template inline SpMat& -SpMat::operator+=(const SpMat& x) +SpMat::operator+=(const SpMat& X) { arma_debug_sigprint(); sync_csc(); - SpMat out = (*this) + x; - - steal_mem(out); + if(X.n_nonzero == 0) + { + arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "addition"); + } + else + { + SpMat tmp = (*this) + X; + + steal_mem(tmp); + } return *this; } @@ -706,15 +713,22 @@ SpMat::operator+=(const SpMat& x) template inline SpMat& -SpMat::operator-=(const SpMat& x) +SpMat::operator-=(const SpMat& X) { arma_debug_sigprint(); sync_csc(); - SpMat out = (*this) - x; - - steal_mem(out); + if(X.n_nonzero == 0) + { + arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "subtraction"); + } + else + { + SpMat tmp = (*this) - X; + + steal_mem(tmp); + } return *this; } @@ -724,15 +738,15 @@ SpMat::operator-=(const SpMat& x) template inline SpMat& -SpMat::operator*=(const SpMat& y) +SpMat::operator*=(const SpMat& X) { arma_debug_sigprint(); sync_csc(); - SpMat z = (*this) * y; + SpMat tmp = (*this) * X; - steal_mem(z); + steal_mem(tmp); return *this; } @@ -743,15 +757,24 @@ SpMat::operator*=(const SpMat& y) template inline SpMat& -SpMat::operator%=(const SpMat& y) +SpMat::operator%=(const SpMat& X) { arma_debug_sigprint(); sync_csc(); - SpMat z = (*this) % y; - - steal_mem(z); + if(X.n_nonzero == 0) + { + arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "element-wise multiplication"); + + (*this).zeros(); + } + else + { + SpMat tmp = (*this) % X; + + steal_mem(tmp); + } return *this; } @@ -1338,9 +1361,16 @@ SpMat::operator+=(const SpSubview& X) sync_csc(); - SpMat tmp = (*this) + X; - - steal_mem(tmp); + if(X.n_nonzero == 0) + { + arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "addition"); + } + else + { + SpMat tmp = (*this) + X; + + steal_mem(tmp); + } return *this; } @@ -1356,9 +1386,16 @@ SpMat::operator-=(const SpSubview& X) sync_csc(); - SpMat tmp = (*this) - X; - - steal_mem(tmp); + if(X.n_nonzero == 0) + { + arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "subtraction"); + } + else + { + SpMat tmp = (*this) - X; + + steal_mem(tmp); + } return *this; } @@ -1368,15 +1405,15 @@ SpMat::operator-=(const SpSubview& X) template inline SpMat& -SpMat::operator*=(const SpSubview& y) +SpMat::operator*=(const SpSubview& X) { arma_debug_sigprint(); sync_csc(); - SpMat z = (*this) * y; + SpMat tmp = (*this) * X; - steal_mem(z); + steal_mem(tmp); return *this; } @@ -1386,15 +1423,24 @@ SpMat::operator*=(const SpSubview& y) template inline SpMat& -SpMat::operator%=(const SpSubview& x) +SpMat::operator%=(const SpSubview& X) { arma_debug_sigprint(); sync_csc(); - SpMat tmp = (*this) % x; - - steal_mem(tmp); + if(X.n_nonzero == 0) + { + arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "element-wise multiplication"); + + (*this).zeros(); + } + else + { + SpMat tmp = (*this) % X; + + steal_mem(tmp); + } return *this; } @@ -1404,16 +1450,18 @@ SpMat::operator%=(const SpSubview& x) template inline SpMat& -SpMat::operator/=(const SpSubview& x) +SpMat::operator/=(const SpSubview& X) { arma_debug_sigprint(); - arma_conform_assert_same_size(n_rows, n_cols, x.n_rows, x.n_cols, "element-wise division"); + // NOTE: use of this function is not advised; it is implemented only for completeness + + arma_conform_assert_same_size(n_rows, n_cols, X.n_rows, X.n_cols, "element-wise division"); - // There is no pretty way to do this. - for(uword elem = 0; elem < n_elem; elem++) + for(uword c = 0; c < n_cols; ++c) + for(uword r = 0; r < n_rows; ++r) { - at(elem) /= x(elem); + at(r, c) /= X.at(r, c); } return *this; @@ -1507,9 +1555,9 @@ SpMat::operator*=(const SpSubview_col_list& X) sync_csc(); - SpMat z = (*this) * X; + SpMat tmp = (*this) * X; - steal_mem(z); + steal_mem(tmp); return *this; } diff --git a/inst/include/current/armadillo_bits/SpSubview_meat.hpp b/inst/include/current/armadillo_bits/SpSubview_meat.hpp index b1c19cbb..22b06c69 100644 --- a/inst/include/current/armadillo_bits/SpSubview_meat.hpp +++ b/inst/include/current/armadillo_bits/SpSubview_meat.hpp @@ -44,6 +44,8 @@ SpSubview::SpSubview(const SpMat& in_m, const uword in_row1, const uword m.sync_csc(); + if( (n_elem == 0) || (m.n_nonzero == 0) ) { return; } // (*this).n_nonzero already set to zero + // count the number of non-zeros in the subview uword count = 0; @@ -127,6 +129,8 @@ SpSubview::operator+=(const eT val) tmp.fill(val); + if(n_nonzero == 0) { return (*this).operator=(tmp); } + return (*this).operator=( (*this) + tmp ); } @@ -319,6 +323,15 @@ SpSubview::operator+=(const Base& x) { arma_debug_sigprint(); + if(n_nonzero == 0) + { + const quasi_unwrap U(x.get_ref()); + + arma_conform_assert_same_size(n_rows, n_cols, U.M.n_rows, U.M.n_cols, "addition"); + + return (*this).operator=(U.M); + } + return (*this).operator=( (*this) + x.get_ref() ); } @@ -369,6 +382,8 @@ SpSubview::operator%=(const Base& x) arma_conform_assert_same_size(sv.n_rows, sv.n_cols, B.n_rows, B.n_cols, "element-wise multiplication"); + if(n_nonzero == 0) { return *this; } + SpMat& sv_m = access::rw(sv.m); sv_m.sync_csc(); @@ -555,6 +570,15 @@ SpSubview::operator+=(const SpBase& x) { arma_debug_sigprint(); + if(n_nonzero == 0) + { + const unwrap_spmat U(x.get_ref()); + + arma_conform_assert_same_size(n_rows, n_cols, U.M.n_rows, U.M.n_cols, "addition"); + + return (*this).operator_equ_common(U.M); + } + // TODO: implement dedicated machinery return (*this).operator=( (*this) + x.get_ref() ); } @@ -596,6 +620,15 @@ SpSubview::operator%=(const SpBase& x) { arma_debug_sigprint(); + if(n_nonzero == 0) + { + const SpProxy P(x.get_ref()); + + arma_conform_assert_same_size(n_rows, n_cols, P.get_n_rows(), P.get_n_cols(), "element-wise multiplication"); + + return *this; + } + // TODO: implement dedicated machinery return (*this).operator=( (*this) % x.get_ref() ); } @@ -648,6 +681,8 @@ SpSubview::for_each(functor F) m.sync_csc(); m.invalidate_cache(); + if(n_nonzero == 0) { return; } + const uword lstart_row = aux_row1; const uword lend_row = aux_row1 + n_rows; @@ -704,6 +739,8 @@ SpSubview::for_each(functor F) const m.sync_csc(); + if(n_nonzero == 0) { return; } + const uword lstart_row = aux_row1; const uword lend_row = aux_row1 + n_rows; @@ -743,6 +780,8 @@ SpSubview::transform(functor F) m.sync_csc(); m.invalidate_cache(); + if(n_nonzero == 0) { return; } + const uword lstart_row = aux_row1; const uword lend_row = aux_row1 + n_rows; @@ -813,6 +852,8 @@ SpSubview::replace(const eT old_val, const eT new_val) m.sync_csc(); m.invalidate_cache(); + if(n_nonzero == 0) { return; } + const uword lstart_row = aux_row1; const uword lend_row = aux_row1 + n_rows; diff --git a/inst/include/current/armadillo_bits/arma_version.hpp b/inst/include/current/armadillo_bits/arma_version.hpp index e92a4e9e..c8e72351 100644 --- a/inst/include/current/armadillo_bits/arma_version.hpp +++ b/inst/include/current/armadillo_bits/arma_version.hpp @@ -23,7 +23,7 @@ #define ARMA_VERSION_MAJOR 15 #define ARMA_VERSION_MINOR 2 -#define ARMA_VERSION_PATCH 3 +#define ARMA_VERSION_PATCH 4 #define ARMA_VERSION_NAME "Medium Roast Deluxe" diff --git a/inst/include/current/armadillo_bits/config.hpp b/inst/include/current/armadillo_bits/config.hpp index 07829d6f..f7e936dc 100644 --- a/inst/include/current/armadillo_bits/config.hpp +++ b/inst/include/current/armadillo_bits/config.hpp @@ -330,18 +330,6 @@ #undef ARMA_OPTIMISE_POWEXPR #endif -#if defined(ARMA_DONT_CHECK_CONFORMANCE) - #if defined(ARMA_CHECK_CONFORMANCE) && (ARMA_WARN_LEVEL >= 2) - #pragma message ("WARNING: conformance checks disabled") - #endif - - #undef ARMA_CHECK_CONFORMANCE -#endif - -#if defined(ARMA_DONT_CHECK_NONFINITE) - #undef ARMA_CHECK_NONFINITE -#endif - #if defined(ARMA_NO_DEBUG) #undef ARMA_DEBUG #undef ARMA_EXTRA_DEBUG @@ -361,6 +349,18 @@ #define ARMA_WARN_LEVEL 3 #endif +#if defined(ARMA_DONT_CHECK_CONFORMANCE) + #if defined(ARMA_CHECK_CONFORMANCE) && (ARMA_WARN_LEVEL >= 2) + #pragma message ("WARNING: conformance checks disabled") + #endif + + #undef ARMA_CHECK_CONFORMANCE +#endif + +#if defined(ARMA_DONT_CHECK_NONFINITE) + #undef ARMA_CHECK_NONFINITE +#endif + #if defined(ARMA_DONT_PRINT_EXCEPTIONS) #undef ARMA_PRINT_EXCEPTIONS #endif diff --git a/inst/include/current/armadillo_bits/constants.hpp b/inst/include/current/armadillo_bits/constants.hpp index a92f06c8..6ad85fb3 100644 --- a/inst/include/current/armadillo_bits/constants.hpp +++ b/inst/include/current/armadillo_bits/constants.hpp @@ -62,11 +62,14 @@ namespace priv } + // + + template static constexpr typename arma_real_only::result - inf(typename arma_real_only::result* junk = nullptr) + pos_inf(typename arma_real_only::result* junk = nullptr) { arma_ignore(junk); @@ -78,13 +81,13 @@ namespace priv static constexpr typename arma_cx_only::result - inf(typename arma_cx_only::result* junk = nullptr) + pos_inf(typename arma_cx_only::result* junk = nullptr) { arma_ignore(junk); typedef typename get_pod_type::result T; - return eT( Datum_helper::inf(), Datum_helper::inf() ); + return eT( Datum_helper::pos_inf(), Datum_helper::pos_inf() ); } @@ -92,12 +95,54 @@ namespace priv static constexpr typename arma_integral_only::result - inf(typename arma_integral_only::result* junk = nullptr) + pos_inf(typename arma_integral_only::result* junk = nullptr) { arma_ignore(junk); return std::numeric_limits::max(); } + + + // + + + template + static + constexpr + typename arma_real_only::result + neg_inf(typename arma_real_only::result* junk = nullptr) + { + arma_ignore(junk); + + return (std::numeric_limits::has_infinity) ? eT(-std::numeric_limits::infinity()) : eT(std::numeric_limits::lowest()); + } + + + template + static + constexpr + typename arma_cx_only::result + neg_inf(typename arma_cx_only::result* junk = nullptr) + { + arma_ignore(junk); + + typedef typename get_pod_type::result T; + + return eT( Datum_helper::neg_inf(), Datum_helper::neg_inf() ); + } + + + template + static + constexpr + typename arma_integral_only::result + neg_inf(typename arma_integral_only::result* junk = nullptr) + { + arma_ignore(junk); + + return std::numeric_limits::lowest(); + } + }; } @@ -125,7 +170,9 @@ struct Datum static const eT log_min; //!< log of the minimum representable value static const eT log_max; //!< log of the maximum representable value static const eT nan; //!< "not a number" - static const eT inf; //!< infinity + static const eT inf; //!< positive infinity + static const eT pos_inf; //!< positive infinity + static const eT neg_inf; //!< negative infinity // @@ -176,7 +223,9 @@ template const eT Datum::eps = std::numeric_limits: template const eT Datum::log_min = std::log(std::numeric_limits::min()); template const eT Datum::log_max = std::log(std::numeric_limits::max()); template const eT Datum::nan = priv::Datum_helper::nan(); -template const eT Datum::inf = priv::Datum_helper::inf(); +template const eT Datum::inf = priv::Datum_helper::pos_inf(); +template const eT Datum::pos_inf = priv::Datum_helper::pos_inf(); +template const eT Datum::neg_inf = priv::Datum_helper::neg_inf(); template const eT Datum::m_u = eT(1.66053906892e-27); template const eT Datum::N_A = eT(6.02214076e23); diff --git a/inst/include/current/armadillo_bits/fill.hpp b/inst/include/current/armadillo_bits/fill.hpp index a61e6ec9..4ddcbac1 100644 --- a/inst/include/current/armadillo_bits/fill.hpp +++ b/inst/include/current/armadillo_bits/fill.hpp @@ -22,22 +22,29 @@ namespace fill { - struct fill_none {}; - struct fill_zeros {}; - struct fill_ones {}; - struct fill_eye {}; - struct fill_randu {}; - struct fill_randn {}; + struct fill_none {}; + struct fill_zeros {}; + struct fill_ones {}; + struct fill_eye {}; + struct fill_randu {}; + struct fill_randn {}; + struct fill_nan {}; + struct fill_pos_inf {}; + struct fill_neg_inf {}; template struct fill_class { inline constexpr fill_class() {} }; - static constexpr fill_class none; - static constexpr fill_class zeros; - static constexpr fill_class ones; - static constexpr fill_class eye; - static constexpr fill_class randu; - static constexpr fill_class randn; + static constexpr fill_class none; + static constexpr fill_class zeros; + static constexpr fill_class ones; + static constexpr fill_class eye; + static constexpr fill_class randu; + static constexpr fill_class randn; + static constexpr fill_class nan; + static constexpr fill_class inf; + static constexpr fill_class pos_inf; + static constexpr fill_class neg_inf; // diff --git a/inst/include/current/armadillo_bits/op_accu_meat.hpp b/inst/include/current/armadillo_bits/op_accu_meat.hpp index 3cd71064..fe7dc918 100644 --- a/inst/include/current/armadillo_bits/op_accu_meat.hpp +++ b/inst/include/current/armadillo_bits/op_accu_meat.hpp @@ -80,6 +80,8 @@ op_accu_mat::apply_proxy_at(const Proxy& P) const uword n_rows = P.get_n_rows(); const uword n_cols = P.get_n_cols(); + if(n_rows == 0) { return eT(0); } + eT val = eT(0); if(n_rows != 1) @@ -146,6 +148,8 @@ op_accu_mat::apply_omit_helper(const Proxy& P, functor is_omitted) const uword n_rows = P.get_n_rows(); const uword n_cols = P.get_n_cols(); + if(n_rows == 0) { return eT_zero; } + for(uword c=0; c < n_cols; ++c) for(uword r=0; r < n_rows; ++r) { @@ -602,6 +606,8 @@ op_accu_mat::apply(const subview& X) const uword X_n_rows = X.n_rows; const uword X_n_cols = X.n_cols; + if( (X_n_rows == 0) || (X_n_cols == 0) ) { return eT(0); } + if(X_n_rows == 1) { const uword X_m_n_rows = X.m.n_rows; @@ -750,6 +756,8 @@ op_accu_cube::apply_proxy_at(const ProxyCube& P) const uword n_cols = P.get_n_cols(); const uword n_slices = P.get_n_slices(); + if( (n_rows == 0) || (n_cols == 0) ) { return eT(0); } + eT val1 = eT(0); eT val2 = eT(0); diff --git a/inst/include/current/armadillo_bits/op_htrans_meat.hpp b/inst/include/current/armadillo_bits/op_htrans_meat.hpp index 99b86ace..7c79bd1d 100644 --- a/inst/include/current/armadillo_bits/op_htrans_meat.hpp +++ b/inst/include/current/armadillo_bits/op_htrans_meat.hpp @@ -65,6 +65,7 @@ op_htrans::apply_mat_noalias(Mat& out, const Mat& A, const typename arma op_htrans::apply_mat_noalias_large(out, A); } else + if( (A_n_rows != 0) && (A_n_cols != 0) ) { eT* outptr = out.memptr(); diff --git a/inst/include/current/armadillo_bits/op_strans_meat.hpp b/inst/include/current/armadillo_bits/op_strans_meat.hpp index 4e3458ea..216a9bc5 100644 --- a/inst/include/current/armadillo_bits/op_strans_meat.hpp +++ b/inst/include/current/armadillo_bits/op_strans_meat.hpp @@ -200,6 +200,7 @@ op_strans::apply_mat_noalias(Mat& out, const TA& A) op_strans::apply_mat_noalias_large(out, A); } else + if( (A_n_rows != 0) && (A_n_cols != 0) ) { eT* outptr = out.memptr(); diff --git a/inst/include/current/armadillo_bits/op_vectorise_meat.hpp b/inst/include/current/armadillo_bits/op_vectorise_meat.hpp index 724ce50c..f185ac4c 100644 --- a/inst/include/current/armadillo_bits/op_vectorise_meat.hpp +++ b/inst/include/current/armadillo_bits/op_vectorise_meat.hpp @@ -176,6 +176,8 @@ op_vectorise_col::apply_subview(Mat& out, const subview& sv) out.set_size(sv.n_elem, 1); + if(sv.n_elem == 0) { return; } + eT* out_mem = out.memptr(); for(uword col=0; col < sv_n_cols; ++col) @@ -201,6 +203,8 @@ op_vectorise_col::apply_proxy(Mat& out, const Proxy& out.set_size(N, 1); + if(N == 0) { return; } + eT* outmem = out.memptr(); if(Proxy::use_at == false) @@ -305,6 +309,8 @@ op_vectorise_row::apply_proxy(Mat& out, const Proxy& out.set_size(1, n_elem); + if(n_elem == 0) { return; } + eT* outmem = out.memptr(); if(n_cols == 1) @@ -471,6 +477,8 @@ op_vectorise_cube_col::apply_proxy(Mat& out, const T1& e out.set_size(N, 1); + if(N == 0) { return; } + eT* outmem = out.memptr(); if(ProxyCube::use_at == false) diff --git a/inst/include/current/armadillo_bits/subview_bones.hpp b/inst/include/current/armadillo_bits/subview_bones.hpp index f3837606..4f727737 100644 --- a/inst/include/current/armadillo_bits/subview_bones.hpp +++ b/inst/include/current/armadillo_bits/subview_bones.hpp @@ -142,6 +142,9 @@ class subview : public Base< eT, subview > arma_inline eT* colptr(const uword in_col); arma_inline const eT* colptr(const uword in_col) const; + arma_inline eT* startptr(); + arma_inline const eT* startptr() const; + template inline bool check_overlap(const subview& x) const; diff --git a/inst/include/current/armadillo_bits/subview_cube_each_meat.hpp b/inst/include/current/armadillo_bits/subview_cube_each_meat.hpp index 739885b9..2df51b5b 100644 --- a/inst/include/current/armadillo_bits/subview_cube_each_meat.hpp +++ b/inst/include/current/armadillo_bits/subview_cube_each_meat.hpp @@ -109,6 +109,8 @@ subview_cube_each1::operator= (const Base& in) const uword p_n_slices = p.n_slices; const uword p_n_elem_slice = p.n_elem_slice; + if(p_n_elem_slice == 0) { return; } + const eT* A_mem = A.memptr(); for(uword i=0; i < p_n_slices; ++i) { arrayops::copy( p.slice_memptr(i), A_mem, p_n_elem_slice ); } @@ -134,6 +136,8 @@ subview_cube_each1::operator+= (const Base& in) const uword p_n_slices = p.n_slices; const uword p_n_elem_slice = p.n_elem_slice; + if(p_n_elem_slice == 0) { return; } + const eT* A_mem = A.memptr(); for(uword i=0; i < p_n_slices; ++i) { arrayops::inplace_plus( p.slice_memptr(i), A_mem, p_n_elem_slice ); } @@ -159,6 +163,8 @@ subview_cube_each1::operator-= (const Base& in) const uword p_n_slices = p.n_slices; const uword p_n_elem_slice = p.n_elem_slice; + if(p_n_elem_slice == 0) { return; } + const eT* A_mem = A.memptr(); for(uword i=0; i < p_n_slices; ++i) { arrayops::inplace_minus( p.slice_memptr(i), A_mem, p_n_elem_slice ); } @@ -184,6 +190,8 @@ subview_cube_each1::operator%= (const Base& in) const uword p_n_slices = p.n_slices; const uword p_n_elem_slice = p.n_elem_slice; + if(p_n_elem_slice == 0) { return; } + const eT* A_mem = A.memptr(); for(uword i=0; i < p_n_slices; ++i) { arrayops::inplace_mul( p.slice_memptr(i), A_mem, p_n_elem_slice ); } @@ -209,6 +217,8 @@ subview_cube_each1::operator/= (const Base& in) const uword p_n_slices = p.n_slices; const uword p_n_elem_slice = p.n_elem_slice; + if(p_n_elem_slice == 0) { return; } + const eT* A_mem = A.memptr(); for(uword i=0; i < p_n_slices; ++i) { arrayops::inplace_div( p.slice_memptr(i), A_mem, p_n_elem_slice ); } @@ -300,6 +310,8 @@ subview_cube_each2::operator= (const Base& in) arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if(p_n_elem_slice == 0) { continue; } + arrayops::copy(p.slice_memptr(slice), A_mem, p_n_elem_slice); } } @@ -339,6 +351,8 @@ subview_cube_each2::operator+= (const Base& in) arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if(p_n_elem_slice == 0) { continue; } + arrayops::inplace_plus(p.slice_memptr(slice), A_mem, p_n_elem_slice); } } @@ -378,6 +392,8 @@ subview_cube_each2::operator-= (const Base& in) arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if(p_n_elem_slice == 0) { continue; } + arrayops::inplace_minus(p.slice_memptr(slice), A_mem, p_n_elem_slice); } } @@ -417,6 +433,8 @@ subview_cube_each2::operator%= (const Base& in) arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if(p_n_elem_slice == 0) { continue; } + arrayops::inplace_mul(p.slice_memptr(slice), A_mem, p_n_elem_slice); } } @@ -456,6 +474,8 @@ subview_cube_each2::operator/= (const Base& in) arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if(p_n_elem_slice == 0) { continue; } + arrayops::inplace_div(p.slice_memptr(slice), A_mem, p_n_elem_slice); } } @@ -492,12 +512,15 @@ subview_cube_each1_aux::operator_plus X.check_size(A); - for(uword i=0; i < p_n_slices; ++i) + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); - const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); - - out_slice = p_slice + A; + for(uword i=0; i < p_n_slices; ++i) + { + Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); + const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); + + out_slice = p_slice + A; + } } return out; @@ -529,12 +552,15 @@ subview_cube_each1_aux::operator_minus X.check_size(A); - for(uword i=0; i < p_n_slices; ++i) + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); - const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); - - out_slice = p_slice - A; + for(uword i=0; i < p_n_slices; ++i) + { + Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); + const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); + + out_slice = p_slice - A; + } } return out; @@ -566,12 +592,15 @@ subview_cube_each1_aux::operator_minus Y.check_size(A); - for(uword i=0; i < p_n_slices; ++i) + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); - const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); - - out_slice = A - p_slice; + for(uword i=0; i < p_n_slices; ++i) + { + Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); + const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); + + out_slice = A - p_slice; + } } return out; @@ -603,12 +632,15 @@ subview_cube_each1_aux::operator_schur X.check_size(A); - for(uword i=0; i < p_n_slices; ++i) + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); - const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); - - out_slice = p_slice % A; + for(uword i=0; i < p_n_slices; ++i) + { + Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); + const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); + + out_slice = p_slice % A; + } } return out; @@ -640,12 +672,15 @@ subview_cube_each1_aux::operator_div X.check_size(A); - for(uword i=0; i < p_n_slices; ++i) + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); - const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); - - out_slice = p_slice / A; + for(uword i=0; i < p_n_slices; ++i) + { + Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); + const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); + + out_slice = p_slice / A; + } } return out; @@ -677,12 +712,15 @@ subview_cube_each1_aux::operator_div Y.check_size(A); - for(uword i=0; i < p_n_slices; ++i) + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); - const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); - - out_slice = A / p_slice; + for(uword i=0; i < p_n_slices; ++i) + { + Mat out_slice( out.slice_memptr(i), p_n_rows, p_n_cols, false, true); + const Mat p_slice(const_cast(p.slice_memptr(i)), p_n_rows, p_n_cols, false, true); + + out_slice = A / p_slice; + } } return out; @@ -706,14 +744,26 @@ subview_cube_each1_aux::operator_times const unwrap tmp(Y.get_ref()); const Mat& M = tmp.M; + if(arma_config::check_conform) + { + if(C.n_cols != M.n_rows) { arma_stop_logic_error("each_slice(): incompatible sizes for matrix multiplication"); } + } + Cube out(C.n_rows, M.n_cols, C.n_slices, arma_nozeros_indicator()); - for(uword i=0; i < C.n_slices; ++i) + if( (C.n_elem == 0) || (M.n_elem == 0) ) { - Mat out_slice( out.slice_memptr(i), C.n_rows, M.n_cols, false, true); - const Mat C_slice(const_cast(C.slice_memptr(i)), C.n_rows, C.n_cols, false, true); - - out_slice = C_slice * M; + out.zeros(); + } + else + { + for(uword i=0; i < C.n_slices; ++i) + { + Mat out_slice( out.slice_memptr(i), C.n_rows, M.n_cols, false, true); + const Mat C_slice(const_cast(C.slice_memptr(i)), C.n_rows, C.n_cols, false, true); + + out_slice = C_slice * M; + } } return out; @@ -737,14 +787,26 @@ subview_cube_each1_aux::operator_times const Cube& C = Y.P; + if(arma_config::check_conform) + { + if(M.n_cols != C.n_rows) { arma_stop_logic_error("each_slice(): incompatible sizes for matrix multiplication"); } + } + Cube out(M.n_rows, C.n_cols, C.n_slices, arma_nozeros_indicator()); - for(uword i=0; i < C.n_slices; ++i) + if( (M.n_elem == 0) || (C.n_elem == 0) ) { - Mat out_slice( out.slice_memptr(i), M.n_rows, C.n_cols, false, true); - const Mat C_slice(const_cast(C.slice_memptr(i)), C.n_rows, C.n_cols, false, true); - - out_slice = M * C_slice; + out.zeros(); + } + else + { + for(uword i=0; i < C.n_slices; ++i) + { + Mat out_slice( out.slice_memptr(i), M.n_rows, C.n_cols, false, true); + const Mat C_slice(const_cast(C.slice_memptr(i)), C.n_rows, C.n_cols, false, true); + + out_slice = M * C_slice; + } } return out; @@ -795,6 +857,8 @@ subview_cube_each2_aux::operator_plus arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if(p_n_elem_slice == 0) { continue; } + arrayops::inplace_plus(out.slice_memptr(slice), A_mem, p_n_elem_slice); } @@ -840,6 +904,8 @@ subview_cube_each2_aux::operator_minus arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if(p_n_elem_slice == 0) { continue; } + arrayops::inplace_minus(out.slice_memptr(slice), A_mem, p_n_elem_slice); } @@ -884,6 +950,8 @@ subview_cube_each2_aux::operator_minus arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if( (p_n_rows == 0) || (p_n_cols == 0) ) { continue; } + Mat out_slice( out.slice_memptr(slice), p_n_rows, p_n_cols, false, true); const Mat p_slice(const_cast(p.slice_memptr(slice)), p_n_rows, p_n_cols, false, true); @@ -932,6 +1000,8 @@ subview_cube_each2_aux::operator_schur arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if(p_n_elem_slice == 0) { continue; } + arrayops::inplace_mul(out.slice_memptr(slice), A_mem, p_n_elem_slice); } @@ -977,6 +1047,8 @@ subview_cube_each2_aux::operator_div arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if(p_n_elem_slice == 0) { continue; } + arrayops::inplace_div(out.slice_memptr(slice), A_mem, p_n_elem_slice); } @@ -1021,6 +1093,8 @@ subview_cube_each2_aux::operator_div arma_conform_check_bounds( (slice >= p_n_slices), "each_slice(): index out of bounds" ); + if( (p_n_rows == 0) || (p_n_cols == 0) ) { continue; } + Mat out_slice( out.slice_memptr(slice), p_n_rows, p_n_cols, false, true); const Mat p_slice(const_cast(p.slice_memptr(slice)), p_n_rows, p_n_cols, false, true); diff --git a/inst/include/current/armadillo_bits/subview_cube_meat.hpp b/inst/include/current/armadillo_bits/subview_cube_meat.hpp index 651eafc0..76b4ae96 100644 --- a/inst/include/current/armadillo_bits/subview_cube_meat.hpp +++ b/inst/include/current/armadillo_bits/subview_cube_meat.hpp @@ -117,6 +117,8 @@ subview_cube::inplace_op(const eT val) const uword t_n_cols = t.n_cols; const uword t_n_slices = t.n_slices; + if( (t_n_rows == 0) || (t_n_cols == 0) ) { return; } + for(uword s=0; s < t_n_slices; ++s) for(uword c=0; c < t_n_cols; ++c) { @@ -150,6 +152,8 @@ subview_cube::inplace_op(const BaseCube& in, const char* identifier) arma_conform_assert_same_size(t, P, identifier); + if( (t_n_rows == 0) || (t_n_cols == 0) || (t_n_slices == 0) ) { return; } + const bool use_mp = arma_config::openmp && ProxyCube::use_mp && mp_gate::eval(t.n_elem); const bool has_overlap = P.has_overlap(t); @@ -258,6 +262,8 @@ subview_cube::inplace_op(const subview_cube& x, const char* identifier) const uword t_n_cols = t.n_cols; const uword t_n_slices = t.n_slices; + if( (t_n_rows == 0) || (t_n_cols == 0) ) { return; } + for(uword s=0; s < t_n_slices; ++s) for(uword c=0; c < t_n_cols; ++c) { @@ -1174,6 +1180,8 @@ subview_cube::replace(const eT old_val, const eT new_val) const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; + if( (local_n_rows == 0) || (local_n_cols == 0) ) { return; } + for(uword slice = 0; slice < local_n_slices; ++slice) { for(uword col = 0; col < local_n_cols; ++col) @@ -1196,6 +1204,8 @@ subview_cube::clean(const typename get_pod_type::result threshold) const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; + if( (local_n_rows == 0) || (local_n_cols == 0) ) { return; } + for(uword slice = 0; slice < local_n_slices; ++slice) { for(uword col = 0; col < local_n_cols; ++col) @@ -1228,6 +1238,8 @@ subview_cube::clamp(const eT min_val, const eT max_val) const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; + if( (local_n_rows == 0) || (local_n_cols == 0) ) { return; } + for(uword slice = 0; slice < local_n_slices; ++slice) { for(uword col = 0; col < local_n_cols; ++col) @@ -1250,6 +1262,8 @@ subview_cube::fill(const eT val) const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; + if( (local_n_rows == 0) || (local_n_cols == 0) ) { return; } + for(uword slice = 0; slice < local_n_slices; ++slice) { for(uword col = 0; col < local_n_cols; ++col) @@ -1272,6 +1286,8 @@ subview_cube::zeros() const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; + if( (local_n_rows == 0) || (local_n_cols == 0) ) { return; } + for(uword slice = 0; slice < local_n_slices; ++slice) { for(uword col = 0; col < local_n_cols; ++col) @@ -1306,6 +1322,8 @@ subview_cube::randu() const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; + if( (local_n_rows == 0) || (local_n_cols == 0) ) { return; } + for(uword slice = 0; slice < local_n_slices; ++slice) { for(uword col = 0; col < local_n_cols; ++col) @@ -1328,6 +1346,8 @@ subview_cube::randn() const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; + if( (local_n_rows == 0) || (local_n_cols == 0) ) { return; } + for(uword slice = 0; slice < local_n_slices; ++slice) { for(uword col = 0; col < local_n_cols; ++col) @@ -1352,11 +1372,14 @@ subview_cube::is_finite() const const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; - for(uword slice = 0; slice < local_n_slices; ++slice) + if( (local_n_rows != 0) && (local_n_cols != 0) ) { - for(uword col = 0; col < local_n_cols; ++col) + for(uword slice = 0; slice < local_n_slices; ++slice) { - if(arrayops::is_finite(slice_colptr(slice,col), local_n_rows) == false) { return false; } + for(uword col = 0; col < local_n_cols; ++col) + { + if(arrayops::is_finite(slice_colptr(slice,col), local_n_rows) == false) { return false; } + } } } @@ -1376,11 +1399,14 @@ subview_cube::is_zero(const typename get_pod_type::result tol) const const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; - for(uword slice = 0; slice < local_n_slices; ++slice) + if( (local_n_rows != 0) && (local_n_cols != 0) ) { - for(uword col = 0; col < local_n_cols; ++col) + for(uword slice = 0; slice < local_n_slices; ++slice) { - if(arrayops::is_zero(slice_colptr(slice,col), local_n_rows, tol) == false) { return false; } + for(uword col = 0; col < local_n_cols; ++col) + { + if(arrayops::is_zero(slice_colptr(slice,col), local_n_rows, tol) == false) { return false; } + } } } @@ -1402,11 +1428,14 @@ subview_cube::has_inf() const const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; - for(uword slice = 0; slice < local_n_slices; ++slice) + if( (local_n_rows != 0) && (local_n_cols != 0) ) { - for(uword col = 0; col < local_n_cols; ++col) + for(uword slice = 0; slice < local_n_slices; ++slice) { - if(arrayops::has_inf(slice_colptr(slice,col), local_n_rows)) { return true; } + for(uword col = 0; col < local_n_cols; ++col) + { + if(arrayops::has_inf(slice_colptr(slice,col), local_n_rows)) { return true; } + } } } @@ -1428,11 +1457,14 @@ subview_cube::has_nan() const const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; - for(uword slice = 0; slice < local_n_slices; ++slice) + if( (local_n_rows != 0) && (local_n_cols != 0) ) { - for(uword col = 0; col < local_n_cols; ++col) + for(uword slice = 0; slice < local_n_slices; ++slice) { - if(arrayops::has_nan(slice_colptr(slice,col), local_n_rows)) { return true; } + for(uword col = 0; col < local_n_cols; ++col) + { + if(arrayops::has_nan(slice_colptr(slice,col), local_n_rows)) { return true; } + } } } @@ -1454,11 +1486,14 @@ subview_cube::has_nonfinite() const const uword local_n_cols = n_cols; const uword local_n_slices = n_slices; - for(uword slice = 0; slice < local_n_slices; ++slice) + if( (local_n_rows != 0) && (local_n_cols != 0) ) { - for(uword col = 0; col < local_n_cols; ++col) + for(uword slice = 0; slice < local_n_slices; ++slice) { - if(arrayops::is_finite(slice_colptr(slice,col), local_n_rows) == false) { return true; } + for(uword col = 0; col < local_n_cols; ++col) + { + if(arrayops::is_finite(slice_colptr(slice,col), local_n_rows) == false) { return true; } + } } } @@ -1614,7 +1649,7 @@ arma_inline eT* subview_cube::slice_colptr(const uword in_slice, const uword in_col) { - return & access::rw((const_cast< Cube& >(m)).mem[ (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 ]); + return access::rwp( m.mem + ( (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 ) ); } @@ -1624,7 +1659,7 @@ arma_inline const eT* subview_cube::slice_colptr(const uword in_slice, const uword in_col) const { - return & m.mem[ (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 ]; + return m.mem + ( (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 ); } @@ -1704,7 +1739,7 @@ void subview_cube::extract(Cube& out, const subview_cube& in) { arma_debug_sigprint(); - + // NOTE: we're assuming that the cube has already been set to the correct size and there is no aliasing; // size setting and alias checking is done by either the Cube constructor or operator=() @@ -1714,6 +1749,8 @@ subview_cube::extract(Cube& out, const subview_cube& in) arma_debug_print(arma_str::format("out.n_rows: %u; out.n_cols: %u; out.n_slices: %u; in.m.n_rows: %u; in.m.n_cols: %u; in.m.n_slices: %u") % out.n_rows % out.n_cols % out.n_slices % in.m.n_rows % in.m.n_cols % in.m.n_slices); + if( (n_rows == 0) || (n_cols == 0) ) { return; } + if( (in.aux_row1 == 0) && (n_rows == in.m.n_rows) ) { for(uword s=0; s < n_slices; ++s) @@ -1747,6 +1784,8 @@ subview_cube::plus_inplace(Cube& out, const subview_cube& in) const uword n_cols = out.n_cols; const uword n_slices = out.n_slices; + if( (n_rows == 0) || (n_cols == 0) ) { return; } + for(uword slice = 0; slice::minus_inplace(Cube& out, const subview_cube& in) const uword n_cols = out.n_cols; const uword n_slices = out.n_slices; + if( (n_rows == 0) || (n_cols == 0) ) { return; } + for(uword slice = 0; slice::schur_inplace(Cube& out, const subview_cube& in) const uword n_cols = out.n_cols; const uword n_slices = out.n_slices; + if( (n_rows == 0) || (n_cols == 0) ) { return; } + for(uword slice = 0; slice::div_inplace(Cube& out, const subview_cube& in) const uword n_cols = out.n_cols; const uword n_slices = out.n_slices; + if( (n_rows == 0) || (n_cols == 0) ) { return; } + for(uword slice = 0; slice::operator= (const Base& in) const uword p_n_rows = p.n_rows; const uword p_n_cols = p.n_cols; + if(p_n_rows == 0) { return; } + if(mode == 0) // each column { for(uword i=0; i < p_n_cols; ++i) @@ -194,6 +196,8 @@ subview_each1::operator+= (const Base& in) const uword p_n_rows = p.n_rows; const uword p_n_cols = p.n_cols; + if(p_n_rows == 0) { return; } + if(mode == 0) // each column { for(uword i=0; i < p_n_cols; ++i) @@ -231,6 +235,8 @@ subview_each1::operator-= (const Base& in) const uword p_n_rows = p.n_rows; const uword p_n_cols = p.n_cols; + if(p_n_rows == 0) { return; } + if(mode == 0) // each column { for(uword i=0; i < p_n_cols; ++i) @@ -268,6 +274,8 @@ subview_each1::operator%= (const Base& in) const uword p_n_rows = p.n_rows; const uword p_n_cols = p.n_cols; + if(p_n_rows == 0) { return; } + if(mode == 0) // each column { for(uword i=0; i < p_n_cols; ++i) @@ -305,6 +313,8 @@ subview_each1::operator/= (const Base& in) const uword p_n_rows = p.n_rows; const uword p_n_cols = p.n_cols; + if(p_n_rows == 0) { return; } + if(mode == 0) // each column { for(uword i=0; i < p_n_cols; ++i) @@ -400,6 +410,8 @@ subview_each2::operator= (const Base& in) arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + arrayops::copy( p.colptr(col), A_mem, p_n_rows ); } } @@ -456,6 +468,8 @@ subview_each2::operator+= (const Base& in) arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + arrayops::inplace_plus( p.colptr(col), A_mem, p_n_rows ); } } @@ -509,6 +523,8 @@ subview_each2::operator-= (const Base& in) arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + arrayops::inplace_minus( p.colptr(col), A_mem, p_n_rows ); } } @@ -562,6 +578,8 @@ subview_each2::operator%= (const Base& in) arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + arrayops::inplace_mul( p.colptr(col), A_mem, p_n_rows ); } } @@ -615,6 +633,8 @@ subview_each2::operator/= (const Base& in) arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + arrayops::inplace_div( p.colptr(col), A_mem, p_n_rows ); } } @@ -666,32 +686,35 @@ subview_each1_aux::operator_plus const eT* A_mem = A.memptr(); - if(mode == 0) // each column + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - for(uword i=0; i < p_n_cols; ++i) + if(mode == 0) // each column { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = p_mem[row] + A_mem[row]; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = p_mem[row] + A_mem[row]; + } } } - } - - if(mode == 1) // each row - { - for(uword i=0; i < p_n_cols; ++i) + + if(mode == 1) // each row { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - const eT A_val = A_mem[i]; - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = p_mem[row] + A_val; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + const eT A_val = A_mem[i]; + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = p_mem[row] + A_val; + } } } } @@ -728,32 +751,35 @@ subview_each1_aux::operator_minus const eT* A_mem = A.memptr(); - if(mode == 0) // each column + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - for(uword i=0; i < p_n_cols; ++i) + if(mode == 0) // each column { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = p_mem[row] - A_mem[row]; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = p_mem[row] - A_mem[row]; + } } } - } - - if(mode == 1) // each row - { - for(uword i=0; i < p_n_cols; ++i) + + if(mode == 1) // each row { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - const eT A_val = A_mem[i]; - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = p_mem[row] - A_val; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + const eT A_val = A_mem[i]; + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = p_mem[row] - A_val; + } } } } @@ -790,32 +816,35 @@ subview_each1_aux::operator_minus const eT* A_mem = A.memptr(); - if(mode == 0) // each column + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - for(uword i=0; i < p_n_cols; ++i) + if(mode == 0) // each column { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = A_mem[row] - p_mem[row]; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = A_mem[row] - p_mem[row]; + } } } - } - - if(mode == 1) // each row - { - for(uword i=0; i < p_n_cols; ++i) + + if(mode == 1) // each row { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - const eT A_val = A_mem[i]; - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = A_val - p_mem[row]; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + const eT A_val = A_mem[i]; + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = A_val - p_mem[row]; + } } } } @@ -852,32 +881,35 @@ subview_each1_aux::operator_schur const eT* A_mem = A.memptr(); - if(mode == 0) // each column + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - for(uword i=0; i < p_n_cols; ++i) + if(mode == 0) // each column { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = p_mem[row] * A_mem[row]; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = p_mem[row] * A_mem[row]; + } } } - } - - if(mode == 1) // each row - { - for(uword i=0; i < p_n_cols; ++i) + + if(mode == 1) // each row { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - const eT A_val = A_mem[i]; - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = p_mem[row] * A_val; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + const eT A_val = A_mem[i]; + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = p_mem[row] * A_val; + } } } } @@ -914,32 +946,35 @@ subview_each1_aux::operator_div const eT* A_mem = A.memptr(); - if(mode == 0) // each column + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - for(uword i=0; i < p_n_cols; ++i) + if(mode == 0) // each column { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = p_mem[row] / A_mem[row]; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = p_mem[row] / A_mem[row]; + } } } - } - - if(mode == 1) // each row - { - for(uword i=0; i < p_n_cols; ++i) + + if(mode == 1) // each row { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - const eT A_val = A_mem[i]; - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = p_mem[row] / A_val; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + const eT A_val = A_mem[i]; + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = p_mem[row] / A_val; + } } } } @@ -976,32 +1011,35 @@ subview_each1_aux::operator_div const eT* A_mem = A.memptr(); - if(mode == 0) // each column + if( (p_n_rows != 0) && (p_n_cols != 0) ) { - for(uword i=0; i < p_n_cols; ++i) + if(mode == 0) // each column { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = A_mem[row] / p_mem[row]; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = A_mem[row] / p_mem[row]; + } } } - } - - if(mode == 1) // each row - { - for(uword i=0; i < p_n_cols; ++i) + + if(mode == 1) // each row { - const eT* p_mem = p.colptr(i); - eT* out_mem = out.colptr(i); - - const eT A_val = A_mem[i]; - - for(uword row=0; row < p_n_rows; ++row) + for(uword i=0; i < p_n_cols; ++i) { - out_mem[row] = A_val / p_mem[row]; + const eT* p_mem = p.colptr(i); + eT* out_mem = out.colptr(i); + + const eT A_val = A_mem[i]; + + for(uword row=0; row < p_n_rows; ++row) + { + out_mem[row] = A_val / p_mem[row]; + } } } } @@ -1058,6 +1096,8 @@ subview_each2_aux::operator_plus arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + arrayops::inplace_plus( out.colptr(col), A_mem, p_n_rows ); } } @@ -1120,6 +1160,8 @@ subview_each2_aux::operator_minus arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + arrayops::inplace_minus( out.colptr(col), A_mem, p_n_rows ); } } @@ -1182,6 +1224,8 @@ subview_each2_aux::operator_minus arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + const eT* p_mem = p.colptr(col); eT* out_mem = out.colptr(col); @@ -1250,6 +1294,8 @@ subview_each2_aux::operator_schur arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + arrayops::inplace_mul( out.colptr(col), A_mem, p_n_rows ); } } @@ -1312,6 +1358,8 @@ subview_each2_aux::operator_div arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + arrayops::inplace_div( out.colptr(col), A_mem, p_n_rows ); } } @@ -1374,6 +1422,8 @@ subview_each2_aux::operator_div arma_conform_check_bounds( (col >= p_n_cols), "each_col(): index out of bounds" ); + if(p_n_rows == 0) { continue; } + const eT* p_mem = p.colptr(col); eT* out_mem = out.colptr(col); diff --git a/inst/include/current/armadillo_bits/subview_meat.hpp b/inst/include/current/armadillo_bits/subview_meat.hpp index 1146b7bd..02bfdc15 100644 --- a/inst/include/current/armadillo_bits/subview_meat.hpp +++ b/inst/include/current/armadillo_bits/subview_meat.hpp @@ -95,6 +95,8 @@ subview::inplace_op(const eT val) const uword s_n_rows = s.n_rows; const uword s_n_cols = s.n_cols; + if( (s_n_rows == 0) || (s_n_cols == 0) ) { return; } + if(s_n_rows == 1) { Mat& A = const_cast< Mat& >(s.m); @@ -151,6 +153,8 @@ subview::inplace_op(const Base& in, const char* identifier) arma_conform_assert_same_size(s, P, identifier); + if( (s_n_rows == 0) || (s_n_cols == 0) ) { return; } + const bool use_mp = arma_config::openmp && Proxy::use_mp && mp_gate::eval(s.n_elem); const bool has_overlap = P.has_overlap(s); @@ -348,8 +352,10 @@ subview::inplace_op(const subview& x, const char* identifier) arma_conform_assert_same_size(s, x, identifier); - const uword s_n_cols = s.n_cols; const uword s_n_rows = s.n_rows; + const uword s_n_cols = s.n_cols; + + if( (s_n_rows == 0) || (s_n_cols == 0) ) { return; } if(s_n_rows == 1) { @@ -606,7 +612,8 @@ subview::operator=(const SpBase& x) // Clear the subview. zeros(); - // Iterate through the sparse subview and set the nonzero values appropriately. + if(p.get_n_nonzero() == 0) { return; } + typename SpProxy::const_iterator_type cit = p.begin(); typename SpProxy::const_iterator_type cit_end = p.end(); @@ -631,7 +638,8 @@ subview::operator+=(const SpBase& x) arma_conform_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols(), "addition"); - // Iterate through the sparse subview and add its values. + if(p.get_n_nonzero() == 0) { return; } + typename SpProxy::const_iterator_type cit = p.begin(); typename SpProxy::const_iterator_type cit_end = p.end(); @@ -656,7 +664,8 @@ subview::operator-=(const SpBase& x) arma_conform_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols(), "subtraction"); - // Iterate through the sparse subview and subtract its values. + if(p.get_n_nonzero() == 0) { return; } + typename SpProxy::const_iterator_type cit = p.begin(); typename SpProxy::const_iterator_type cit_end = p.end(); @@ -725,13 +734,12 @@ subview::operator/=(const SpBase& x) { arma_debug_sigprint(); + // NOTE: use of this function is not advised; it is implemented only for completeness + const SpProxy p(x.get_ref()); arma_conform_assert_same_size(n_rows, n_cols, p.get_n_rows(), p.get_n_cols(), "element-wise division"); - // This is probably going to fill your subview with a bunch of NaNs, - // so I'm not going to bother to implement it fast. - // You can have slow NaNs. They're fine too. for(uword c = 0; c < n_cols; ++c) for(uword r = 0; r < n_rows; ++r) { @@ -968,8 +976,10 @@ subview::replace(const eT old_val, const eT new_val) subview& s = *this; - const uword s_n_cols = s.n_cols; const uword s_n_rows = s.n_rows; + const uword s_n_cols = s.n_cols; + + if( (s_n_rows == 0) || (s_n_cols == 0) ) { return; } if(s_n_rows == 1) { @@ -1018,8 +1028,10 @@ subview::clean(const typename get_pod_type::result threshold) subview& s = *this; - const uword s_n_cols = s.n_cols; const uword s_n_rows = s.n_rows; + const uword s_n_cols = s.n_cols; + + if(s_n_rows == 0) { return; } for(uword ucol=0; ucol < s_n_cols; ++ucol) { @@ -1048,8 +1060,10 @@ subview::clamp(const eT min_val, const eT max_val) subview& s = *this; - const uword s_n_cols = s.n_cols; const uword s_n_rows = s.n_rows; + const uword s_n_cols = s.n_cols; + + if(s_n_rows == 0) { return; } for(uword ucol=0; ucol < s_n_cols; ++ucol) { @@ -1068,8 +1082,10 @@ subview::fill(const eT val) subview& s = *this; - const uword s_n_cols = s.n_cols; const uword s_n_rows = s.n_rows; + const uword s_n_cols = s.n_cols; + + if( (s_n_rows == 0) || (s_n_cols == 0) ) { return; } if(s_n_rows == 1) { @@ -1109,7 +1125,42 @@ subview::zeros() { arma_debug_sigprint(); - (*this).fill(eT(0)); + subview& s = *this; + + const uword s_n_rows = s.n_rows; + const uword s_n_cols = s.n_cols; + + if( (s_n_rows == 0) || (s_n_cols == 0) ) { return; } + + if(s_n_rows == 1) + { + Mat& A = const_cast< Mat& >(s.m); + + const uword A_n_rows = A.n_rows; + + eT* Aptr = &(A.at(s.aux_row1,s.aux_col1)); + + constexpr eT eT_zero = eT(0); + + for(uword ii=0; ii < s_n_cols; ++ii) + { + (*Aptr) = eT_zero; Aptr += A_n_rows; + } + } + else + { + if( (s.aux_row1 == 0) && (s_n_rows == s.m.n_rows) ) + { + arrayops::fill_zeros( s.colptr(0), s.n_elem ); + } + else + { + for(uword ucol=0; ucol < s_n_cols; ++ucol) + { + arrayops::fill_zeros( s.colptr(ucol), s_n_rows ); + } + } + } } @@ -1157,6 +1208,8 @@ subview::randu() const uword s_n_rows = s.n_rows; const uword s_n_cols = s.n_cols; + if( (s_n_rows == 0) || (s_n_cols == 0) ) { return; } + if(s_n_rows == 1) { podarray tmp(s_n_cols); @@ -1197,6 +1250,8 @@ subview::randn() const uword s_n_rows = s.n_rows; const uword s_n_cols = s.n_cols; + if( (s_n_rows == 0) || (s_n_cols == 0) ) { return; } + if(s_n_rows == 1) { podarray tmp(s_n_cols); @@ -1410,7 +1465,7 @@ arma_inline eT* subview::colptr(const uword in_col) { - return & access::rw((const_cast< Mat& >(m)).mem[ (in_col + aux_col1)*m.n_rows + aux_row1 ]); + return access::rwp( m.mem + ((in_col + aux_col1)*m.n_rows + aux_row1) ); } @@ -1420,7 +1475,27 @@ arma_inline const eT* subview::colptr(const uword in_col) const { - return & m.mem[ (in_col + aux_col1)*m.n_rows + aux_row1 ]; + return m.mem + ((in_col + aux_col1)*m.n_rows + aux_row1); + } + + + +template +arma_inline +eT* +subview::startptr() + { + return access::rwp( m.mem + (aux_col1*m.n_rows + aux_row1) ); + } + + + +template +arma_inline +const eT* +subview::startptr() const + { + return m.mem + (aux_col1*m.n_rows + aux_row1); } @@ -1483,9 +1558,12 @@ subview::is_finite() const const uword local_n_rows = n_rows; const uword local_n_cols = n_cols; - for(uword ii=0; ii::is_zero(const typename get_pod_type::result tol) const const uword local_n_rows = n_rows; const uword local_n_cols = n_cols; - for(uword ii=0; ii::has_inf() const const uword local_n_rows = n_rows; const uword local_n_cols = n_cols; - for(uword ii=0; ii::has_nan() const const uword local_n_rows = n_rows; const uword local_n_cols = n_cols; - for(uword ii=0; ii::has_nonfinite() const const uword local_n_rows = n_rows; const uword local_n_cols = n_cols; - for(uword ii=0; ii::extract(Mat& out, const subview& in) // NOTE: we're assuming that the matrix has already been set to the correct size and there is no aliasing; // size setting and alias checking is done by either the Mat constructor or operator=() - const uword n_rows = in.n_rows; // number of rows in the subview - const uword n_cols = in.n_cols; // number of columns in the subview - arma_debug_print(arma_str::format("out.n_rows: %u; out.n_cols: %u; in.m.n_rows: %u; in.m.n_cols: %u") % out.n_rows % out.n_cols % in.m.n_rows % in.m.n_cols ); + const uword n_rows = in.n_rows; + const uword n_cols = in.n_cols; + + if( (n_rows == 0) || (n_cols == 0) ) { return; } + + if(n_cols == 1) + { + arma_debug_print("subview::extract(): copying col"); + + // in.colptr(0) is the first column of the subview, taking into account any row offset + arrayops::copy( out.memptr(), in.colptr(0), n_rows ); + + return; + } - if(in.is_vec()) + if(n_rows == 1) { - if(n_cols == 1) // a column vector + arma_debug_print("subview::extract(): copying row"); + + eT* out_mem = out.memptr(); + + const uword X_n_rows = in.m.n_rows; + + const eT* Xptr = &(in.m.at(in.aux_row1,in.aux_col1)); + + uword j; + + for(j=1; j < n_cols; j+=2) { - arma_debug_print("subview::extract(): copying col"); + const eT tmp1 = (*Xptr); Xptr += X_n_rows; + const eT tmp2 = (*Xptr); Xptr += X_n_rows; - // in.colptr(0) is the first column of the subview, taking into account any row offset - arrayops::copy( out.memptr(), in.colptr(0), n_rows ); + (*out_mem) = tmp1; out_mem++; + (*out_mem) = tmp2; out_mem++; } - else - if(n_rows == 1) // a row vector + + if((j-1) < n_cols) { - arma_debug_print("subview::extract(): copying row)"); - - eT* out_mem = out.memptr(); - - const uword X_n_rows = in.m.n_rows; - - const eT* Xptr = &(in.m.at(in.aux_row1,in.aux_col1)); - - uword j; - - for(j=1; j < n_cols; j+=2) - { - const eT tmp1 = (*Xptr); Xptr += X_n_rows; - const eT tmp2 = (*Xptr); Xptr += X_n_rows; - - (*out_mem) = tmp1; out_mem++; - (*out_mem) = tmp2; out_mem++; - } - - if((j-1) < n_cols) - { - (*out_mem) = (*Xptr); - } + (*out_mem) = (*Xptr); } + + return; } - else // general submatrix + + if( (in.aux_row1 == 0) && (n_rows == in.m.n_rows) ) { - arma_debug_print("subview::extract(): general submatrix"); + arma_debug_print("subview::extract(): contiguous submatrix"); - if( (in.aux_row1 == 0) && (n_rows == in.m.n_rows) ) - { - arrayops::copy( out.memptr(), in.colptr(0), in.n_elem ); - } - else - { - for(uword col=0; col < n_cols; ++col) - { - arrayops::copy( out.colptr(col), in.colptr(col), n_rows ); - } - } + arrayops::copy( out.memptr(), in.colptr(0), in.n_elem ); + + return; + } + + arma_debug_print("subview::extract(): general submatrix"); + + for(uword col=0; col < n_cols; ++col) + { + arrayops::copy( out.colptr(col), in.colptr(col), n_rows ); } } @@ -1666,6 +1758,8 @@ subview::plus_inplace(Mat& out, const subview& in) const uword n_rows = in.n_rows; const uword n_cols = in.n_cols; + if( (n_rows == 0) || (n_cols == 0) ) { return; } + if(n_rows == 1) { eT* out_mem = out.memptr(); @@ -1714,6 +1808,8 @@ subview::minus_inplace(Mat& out, const subview& in) const uword n_rows = in.n_rows; const uword n_cols = in.n_cols; + if( (n_rows == 0) || (n_cols == 0) ) { return; } + if(n_rows == 1) { eT* out_mem = out.memptr(); @@ -1762,6 +1858,8 @@ subview::schur_inplace(Mat& out, const subview& in) const uword n_rows = in.n_rows; const uword n_cols = in.n_cols; + if( (n_rows == 0) || (n_cols == 0) ) { return; } + if(n_rows == 1) { eT* out_mem = out.memptr(); @@ -1810,6 +1908,8 @@ subview::div_inplace(Mat& out, const subview& in) const uword n_rows = in.n_rows; const uword n_cols = in.n_cols; + if( (n_rows == 0) || (n_cols == 0) ) { return; } + if(n_rows == 1) { eT* out_mem = out.memptr(); @@ -2377,9 +2477,12 @@ subview::each_col(const std::function< void(Col&) >& F) { arma_debug_sigprint(); - for(uword ii=0; ii < n_cols; ++ii) + const uword local_n_rows = n_rows; + const uword local_n_cols = n_cols; + + for(uword ii=0; ii < local_n_cols; ++ii) { - Col tmp(colptr(ii), n_rows, false, true); + Col tmp(colptr(ii), local_n_rows, false, true); F(tmp); } } @@ -2393,9 +2496,12 @@ subview::each_col(const std::function< void(const Col&) >& F) const { arma_debug_sigprint(); - for(uword ii=0; ii < n_cols; ++ii) + const uword local_n_rows = n_rows; + const uword local_n_cols = n_cols; + + for(uword ii=0; ii < local_n_cols; ++ii) { - const Col tmp(colptr(ii), n_rows, false, true); + const Col tmp(colptr(ii), local_n_rows, false, true); F(tmp); } } @@ -2410,20 +2516,23 @@ subview::each_row(const std::function< void(Row&) >& F) { arma_debug_sigprint(); - podarray array1(n_cols); - podarray array2(n_cols); + const uword local_n_rows = n_rows; + const uword local_n_cols = n_cols; - Row tmp1( array1.memptr(), n_cols, false, true ); - Row tmp2( array2.memptr(), n_cols, false, true ); + podarray array1(local_n_cols); + podarray array2(local_n_cols); + + Row tmp1( array1.memptr(), local_n_cols, false, true ); + Row tmp2( array2.memptr(), local_n_cols, false, true ); eT* tmp1_mem = tmp1.memptr(); eT* tmp2_mem = tmp2.memptr(); uword ii, jj; - for(ii=0, jj=1; jj < n_rows; ii+=2, jj+=2) + for(ii=0, jj=1; jj < local_n_rows; ii+=2, jj+=2) { - for(uword col_id = 0; col_id < n_cols; ++col_id) + for(uword col_id = 0; col_id < local_n_cols; ++col_id) { const eT* col_mem = colptr(col_id); @@ -2434,7 +2543,7 @@ subview::each_row(const std::function< void(Row&) >& F) F(tmp1); F(tmp2); - for(uword col_id = 0; col_id < n_cols; ++col_id) + for(uword col_id = 0; col_id < local_n_cols; ++col_id) { eT* col_mem = colptr(col_id); @@ -2443,7 +2552,7 @@ subview::each_row(const std::function< void(Row&) >& F) } } - if(ii < n_rows) + if(ii < local_n_rows) { tmp1 = (*this).row(ii); @@ -2462,20 +2571,23 @@ subview::each_row(const std::function< void(const Row&) >& F) const { arma_debug_sigprint(); - podarray array1(n_cols); - podarray array2(n_cols); + const uword local_n_rows = n_rows; + const uword local_n_cols = n_cols; + + podarray array1(local_n_cols); + podarray array2(local_n_cols); - Row tmp1( array1.memptr(), n_cols, false, true ); - Row tmp2( array2.memptr(), n_cols, false, true ); + Row tmp1( array1.memptr(), local_n_cols, false, true ); + Row tmp2( array2.memptr(), local_n_cols, false, true ); eT* tmp1_mem = tmp1.memptr(); eT* tmp2_mem = tmp2.memptr(); uword ii, jj; - for(ii=0, jj=1; jj < n_rows; ii+=2, jj+=2) + for(ii=0, jj=1; jj < local_n_rows; ii+=2, jj+=2) { - for(uword col_id = 0; col_id < n_cols; ++col_id) + for(uword col_id = 0; col_id < local_n_cols; ++col_id) { const eT* col_mem = colptr(col_id); @@ -2487,7 +2599,7 @@ subview::each_row(const std::function< void(const Row&) >& F) const F(tmp2); } - if(ii < n_rows) + if(ii < local_n_rows) { tmp1 = (*this).row(ii); @@ -3264,7 +3376,7 @@ template inline subview_col::subview_col(const Mat& in_m, const uword in_col) : subview(in_m, 0, in_col, in_m.n_rows, 1) - , colmem(subview::colptr(0)) + , colmem(subview::startptr()) { arma_debug_sigprint(); } @@ -3275,7 +3387,7 @@ template inline subview_col::subview_col(const Mat& in_m, const uword in_col, const uword in_row1, const uword in_n_rows) : subview(in_m, in_row1, in_col, in_n_rows, 1) - , colmem(subview::colptr(0)) + , colmem(subview::startptr()) { arma_debug_sigprint(); } @@ -4268,7 +4380,7 @@ template inline subview_row::subview_row(const Mat& in_m, const uword in_row) : subview(in_m, in_row, 0, 1, in_m.n_cols) - , rowmem(subview::colptr(0)) + , rowmem(subview::startptr()) { arma_debug_sigprint(); } @@ -4279,7 +4391,7 @@ template inline subview_row::subview_row(const Mat& in_m, const uword in_row, const uword in_col1, const uword in_n_cols) : subview(in_m, in_row, in_col1, 1, in_n_cols) - , rowmem(subview::colptr(0)) + , rowmem(subview::startptr()) { arma_debug_sigprint(); }