From 3f920db2132df7d6079ce061c337f4e14a88c4a2 Mon Sep 17 00:00:00 2001 From: Dirk Eddelbuettel Date: Tue, 6 Feb 2024 13:45:08 -0600 Subject: [PATCH] Armadillo 12.8.0 --- inst/include/armadillo | 2 + inst/include/armadillo_bits/arma_config.hpp | 7 + inst/include/armadillo_bits/arma_version.hpp | 6 +- .../include/armadillo_bits/compiler_check.hpp | 10 + inst/include/armadillo_bits/debug.hpp | 4 +- inst/include/armadillo_bits/diskio_meat.hpp | 116 +++++++-- inst/include/armadillo_bits/eop_aux.hpp | 4 + .../include/armadillo_bits/eop_core_bones.hpp | 1 + inst/include/armadillo_bits/eop_core_meat.hpp | 3 + inst/include/armadillo_bits/fn_conv_to.hpp | 240 +++++++++++++++++- inst/include/armadillo_bits/fn_elem.hpp | 42 +++ inst/include/armadillo_bits/fn_shift.hpp | 35 +++ inst/include/armadillo_bits/op_pinv_meat.hpp | 15 +- inst/include/armadillo_bits/op_rank_meat.hpp | 19 +- inst/include/armadillo_bits/op_shift_meat.hpp | 2 + .../include/armadillo_bits/sp_auxlib_meat.hpp | 55 +--- .../armadillo_bits/spop_misc_bones.hpp | 11 + .../include/armadillo_bits/spop_misc_meat.hpp | 27 +- .../armadillo_bits/spop_shift_bones.hpp | 35 +++ .../armadillo_bits/spop_shift_meat.hpp | 93 +++++++ inst/include/armadillo_bits/traits.hpp | 39 +++ 21 files changed, 669 insertions(+), 97 deletions(-) create mode 100644 inst/include/armadillo_bits/spop_shift_bones.hpp create mode 100644 inst/include/armadillo_bits/spop_shift_meat.hpp diff --git a/inst/include/armadillo b/inst/include/armadillo index 59694d1c..ed4ed3f5 100644 --- a/inst/include/armadillo +++ b/inst/include/armadillo @@ -369,6 +369,7 @@ namespace arma #include "armadillo_bits/spop_vectorise_bones.hpp" #include "armadillo_bits/spop_norm_bones.hpp" #include "armadillo_bits/spop_vecnorm_bones.hpp" + #include "armadillo_bits/spop_shift_bones.hpp" #include "armadillo_bits/spglue_plus_bones.hpp" #include "armadillo_bits/spglue_minus_bones.hpp" @@ -816,6 +817,7 @@ namespace arma #include "armadillo_bits/spop_vectorise_meat.hpp" #include "armadillo_bits/spop_norm_meat.hpp" #include "armadillo_bits/spop_vecnorm_meat.hpp" + #include "armadillo_bits/spop_shift_meat.hpp" #include "armadillo_bits/spglue_plus_meat.hpp" #include "armadillo_bits/spglue_minus_meat.hpp" diff --git a/inst/include/armadillo_bits/arma_config.hpp b/inst/include/armadillo_bits/arma_config.hpp index 3670199e..139efa7d 100644 --- a/inst/include/armadillo_bits/arma_config.hpp +++ b/inst/include/armadillo_bits/arma_config.hpp @@ -181,6 +181,13 @@ struct arma_config #endif + #if defined(ARMA_HAVE_CXX23) + static constexpr bool cxx23 = true; + #else + static constexpr bool cxx23 = false; + #endif + + #if (!defined(ARMA_DONT_USE_STD_MUTEX)) static constexpr bool std_mutex = true; #else diff --git a/inst/include/armadillo_bits/arma_version.hpp b/inst/include/armadillo_bits/arma_version.hpp index a335bb3b..092454dc 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 12 -#define ARMA_VERSION_MINOR 6 -#define ARMA_VERSION_PATCH 7 -#define ARMA_VERSION_NAME "Cortisol Retox" +#define ARMA_VERSION_MINOR 8 +#define ARMA_VERSION_PATCH 0 +#define ARMA_VERSION_NAME "Cortisol Injector" diff --git a/inst/include/armadillo_bits/compiler_check.hpp b/inst/include/armadillo_bits/compiler_check.hpp index 8a653d24..b6fbd8a6 100644 --- a/inst/include/armadillo_bits/compiler_check.hpp +++ b/inst/include/armadillo_bits/compiler_check.hpp @@ -20,6 +20,7 @@ #undef ARMA_HAVE_CXX14 #undef ARMA_HAVE_CXX17 #undef ARMA_HAVE_CXX20 +#undef ARMA_HAVE_CXX23 #if (__cplusplus >= 201103L) #define ARMA_HAVE_CXX11 @@ -37,6 +38,10 @@ #define ARMA_HAVE_CXX20 #endif +#if (__cplusplus >= 202302L) + #define ARMA_HAVE_CXX23 +#endif + // MS really can't get its proverbial shit together #if defined(_MSVC_LANG) @@ -59,6 +64,11 @@ #define ARMA_HAVE_CXX20 #endif + #if (_MSVC_LANG >= 202302L) + #undef ARMA_HAVE_CXX23 + #define ARMA_HAVE_CXX23 + #endif + #endif diff --git a/inst/include/armadillo_bits/debug.hpp b/inst/include/armadillo_bits/debug.hpp index 7a0b95c2..b436c432 100644 --- a/inst/include/armadillo_bits/debug.hpp +++ b/inst/include/armadillo_bits/debug.hpp @@ -263,8 +263,7 @@ arma_print(const T1& x, const T2& y, const T3& z) // // arma_sigprint -//! print a message to the log stream with a preceding @ character. -//! by default the log stream is cout. +//! print a message to the cerr stream with a preceding @ character. //! used for printing the signature of a function //! (see the arma_extra_debug_sigprint macro) inline @@ -1424,6 +1423,7 @@ arma_assert_atlas_size(const T1& A, const T2& B) out << "@ arma_config::cxx14 = " << arma_config::cxx14 << '\n'; out << "@ arma_config::cxx17 = " << arma_config::cxx17 << '\n'; out << "@ arma_config::cxx20 = " << arma_config::cxx20 << '\n'; + out << "@ arma_config::cxx23 = " << arma_config::cxx23 << '\n'; out << "@ arma_config::std_mutex = " << arma_config::std_mutex << '\n'; out << "@ arma_config::posix = " << arma_config::posix << '\n'; out << "@ arma_config::openmp = " << arma_config::openmp << '\n'; diff --git a/inst/include/armadillo_bits/diskio_meat.hpp b/inst/include/armadillo_bits/diskio_meat.hpp index 4f716cc0..1658e596 100644 --- a/inst/include/armadillo_bits/diskio_meat.hpp +++ b/inst/include/armadillo_bits/diskio_meat.hpp @@ -915,11 +915,18 @@ diskio::save_csv_ascii(const Mat& x, std::ostream& f, const char separator) uword x_n_rows = x.n_rows; uword x_n_cols = x.n_cols; + const eT eT_int_lowest = eT(std::numeric_limits::lowest()); + const eT eT_int_max = eT(std::numeric_limits::max()); + for(uword row=0; row < x_n_rows; ++row) { for(uword col=0; col < x_n_cols; ++col) { - arma_ostream::raw_print_elem(f, x.at(row,col)); + const eT val = x.at(row,col); + + const bool is_real_int = (is_real::yes) && arma_isfinite(val) && (val > eT_int_lowest) && (val < eT_int_max) && (eT(int(val)) == val); + + (is_real_int) ? arma_ostream::raw_print_elem(f, int(val)) : arma_ostream::raw_print_elem(f, val); if( col < (x_n_cols-1) ) { f.put(separator); } } @@ -950,6 +957,9 @@ diskio::save_csv_ascii(const Mat< std::complex >& x, std::ostream& f, const c diskio::prepare_stream(f); + const T T_int_lowest = T(std::numeric_limits::lowest()); + const T T_int_max = T(std::numeric_limits::max()); + uword x_n_rows = x.n_rows; uword x_n_cols = x.n_cols; @@ -959,14 +969,20 @@ diskio::save_csv_ascii(const Mat< std::complex >& x, std::ostream& f, const c { const eT& val = x.at(row,col); - const T tmp_r = std::real(val); - const T tmp_i = std::imag(val); - const T tmp_i_abs = (tmp_i < T(0)) ? T(-tmp_i) : T(tmp_i); - const char tmp_sign = (tmp_i < T(0)) ? char('-') : char('+'); + const T val_r = std::real(val); + const T val_i = std::imag(val); + const T abs_i = (val_i < T(0)) ? T(-val_i) : T(val_i); + const char sgn_i = (val_i < T(0)) ? char('-') : char('+'); + + const bool val_r_is_real_int = (is_real::yes) && arma_isfinite(val_r) && (val_r > T_int_lowest) && (val_r < T_int_max) && (T(int(val_r)) == val_r); + const bool abs_i_is_real_int = (is_real::yes) && arma_isfinite(abs_i) && (abs_i < T_int_max) && (T(int(abs_i)) == abs_i); + + (val_r_is_real_int) ? arma_ostream::raw_print_elem(f, int(val_r)) : arma_ostream::raw_print_elem(f, val_r); + + f.put(sgn_i); + + (abs_i_is_real_int) ? arma_ostream::raw_print_elem(f, int(abs_i)) : arma_ostream::raw_print_elem(f, abs_i); - arma_ostream::raw_print_elem(f, tmp_r ); - f.put(tmp_sign); - arma_ostream::raw_print_elem(f, tmp_i_abs); f.put('i'); if( col < (x_n_cols-1) ) { f.put(separator); } @@ -1025,15 +1041,25 @@ diskio::save_coord_ascii(const Mat& x, std::ostream& f) diskio::prepare_stream(f); + const eT eT_zero = eT(0); + const eT eT_int_lowest = eT(std::numeric_limits::lowest()); + const eT eT_int_max = eT(std::numeric_limits::max()); + for(uword col=0; col < x.n_cols; ++col) for(uword row=0; row < x.n_rows; ++row) { const eT val = x.at(row,col); - if(val != eT(0)) - { - f << row << ' ' << col << ' ' << val << '\n'; - } + if(val == eT_zero) { continue; } + + f << row; f.put(' '); + f << col; f.put(' '); + + const bool is_real_int = (is_real::yes) && arma_isfinite(val) && (val > eT_int_lowest) && (val < eT_int_max) && (eT(int(val)) == val); + + (is_real_int) ? arma_ostream::raw_print_elem(f, int(val)) : arma_ostream::raw_print_elem(f, val); + + f.put('\n'); } // make sure it's possible to figure out the matrix size later @@ -1070,17 +1096,33 @@ diskio::save_coord_ascii(const Mat< std::complex >& x, std::ostream& f) diskio::prepare_stream(f); - const eT eT_zero = eT(0); + const eT eT_zero = eT(0); + const T T_int_lowest = T(std::numeric_limits::lowest()); + const T T_int_max = T(std::numeric_limits::max()); for(uword col=0; col < x.n_cols; ++col) for(uword row=0; row < x.n_rows; ++row) { const eT& val = x.at(row,col); - if(val != eT_zero) - { - f << row << ' ' << col << ' ' << val.real() << ' ' << val.imag() << '\n'; - } + if(val == eT_zero) { continue; } + + f << row; f.put(' '); + f << col; f.put(' '); + + const T val_r = std::real(val); + const T val_i = std::imag(val); + + const bool val_r_is_real_int = (is_real::yes) && arma_isfinite(val_r) && (val_r > T_int_lowest) && (val_r < T_int_max) && (T(int(val_r)) == val_r); + const bool val_i_is_real_int = (is_real::yes) && arma_isfinite(val_i) && (val_i > T_int_lowest) && (val_i < T_int_max) && (T(int(val_i)) == val_i); + + (val_r_is_real_int) ? arma_ostream::raw_print_elem(f, int(val_r)) : arma_ostream::raw_print_elem(f, val_r); + + f.put(' '); + + (val_i_is_real_int) ? arma_ostream::raw_print_elem(f, int(val_i)) : arma_ostream::raw_print_elem(f, val_i); + + f.put('\n'); } // make sure it's possible to figure out the matrix size later @@ -2904,7 +2946,9 @@ diskio::save_csv_ascii(const SpMat& x, std::ostream& f, const char separator uword x_n_rows = x.n_rows; uword x_n_cols = x.n_cols; - const eT eT_zero = eT(0); + const eT eT_zero = eT(0); + const eT eT_int_lowest = eT(std::numeric_limits::lowest()); + const eT eT_int_max = eT(std::numeric_limits::max()); for(uword row=0; row < x_n_rows; ++row) { @@ -2918,7 +2962,9 @@ diskio::save_csv_ascii(const SpMat& x, std::ostream& f, const char separator } else { - arma_ostream::raw_print_elem(f, val); + const bool is_real_int = (is_real::yes) && arma_isfinite(val) && (val > eT_int_lowest) && (val < eT_int_max) && (eT(int(val)) == val); + + (is_real_int) ? arma_ostream::raw_print_elem(f, int(val)) : arma_ostream::raw_print_elem(f, val); } if( col < (x_n_cols-1) ) { f.put(separator); } @@ -2998,14 +3044,24 @@ diskio::save_coord_ascii(const SpMat& x, std::ostream& f) diskio::prepare_stream(f); + const eT eT_int_lowest = eT(std::numeric_limits::lowest()); + const eT eT_int_max = eT(std::numeric_limits::max()); + typename SpMat::const_iterator iter = x.begin(); typename SpMat::const_iterator iter_end = x.end(); for(; iter != iter_end; ++iter) { + f << iter.row(); f.put(' '); + f << iter.col(); f.put(' '); + const eT val = (*iter); - f << iter.row() << ' ' << iter.col() << ' ' << val << '\n'; + const bool is_real_int = (is_real::yes) && arma_isfinite(val) && (val > eT_int_lowest) && (val < eT_int_max) && (eT(int(val)) == val); + + (is_real_int) ? arma_ostream::raw_print_elem(f, int(val)) : arma_ostream::raw_print_elem(f, val); + + f.put('\n'); } @@ -3044,14 +3100,32 @@ diskio::save_coord_ascii(const SpMat< std::complex >& x, std::ostream& f) diskio::prepare_stream(f); + const T T_int_lowest = T(std::numeric_limits::lowest()); + const T T_int_max = T(std::numeric_limits::max()); + typename SpMat::const_iterator iter = x.begin(); typename SpMat::const_iterator iter_end = x.end(); for(; iter != iter_end; ++iter) { + f << iter.row(); f.put(' '); + f << iter.col(); f.put(' '); + const eT val = (*iter); - f << iter.row() << ' ' << iter.col() << ' ' << val.real() << ' ' << val.imag() << '\n'; + const T val_r = std::real(val); + const T val_i = std::imag(val); + + const bool val_r_is_real_int = (is_real::yes) && arma_isfinite(val_r) && (val_r > T_int_lowest) && (val_r < T_int_max) && (T(int(val_r)) == val_r); + const bool val_i_is_real_int = (is_real::yes) && arma_isfinite(val_i) && (val_i > T_int_lowest) && (val_i < T_int_max) && (T(int(val_i)) == val_i); + + (val_r_is_real_int) ? arma_ostream::raw_print_elem(f, int(val_r)) : arma_ostream::raw_print_elem(f, val_r); + + f.put(' '); + + (val_i_is_real_int) ? arma_ostream::raw_print_elem(f, int(val_i)) : arma_ostream::raw_print_elem(f, val_i); + + f.put('\n'); } // make sure it's possible to figure out the matrix size later diff --git a/inst/include/armadillo_bits/eop_aux.hpp b/inst/include/armadillo_bits/eop_aux.hpp index 2b66ef2a..107ff2ae 100644 --- a/inst/include/armadillo_bits/eop_aux.hpp +++ b/inst/include/armadillo_bits/eop_aux.hpp @@ -116,6 +116,10 @@ class eop_aux template arma_inline static typename arma_real_only::result arma_abs (const eT x) { return std::abs(x); } template arma_inline static typename arma_real_only< T>::result arma_abs (const std::complex& x) { return std::abs(x); } + template arma_inline static typename arma_integral_only::result cbrt (const eT x) { return eT( std::cbrt(double(x)) ); } + template arma_inline static typename arma_real_only::result cbrt (const eT x) { return std::cbrt(x); } + template arma_inline static typename arma_cx_only::result cbrt (const eT& x) { arma_ignore(x); return eT(0); } + template arma_inline static typename arma_integral_only::result erf (const eT x) { return eT( std::erf(double(x)) ); } template arma_inline static typename arma_real_only::result erf (const eT x) { return std::erf(x); } template arma_inline static typename arma_cx_only::result erf (const eT& x) { arma_ignore(x); return eT(0); } diff --git a/inst/include/armadillo_bits/eop_core_bones.hpp b/inst/include/armadillo_bits/eop_core_bones.hpp index 8b9c75cc..24663504 100644 --- a/inst/include/armadillo_bits/eop_core_bones.hpp +++ b/inst/include/armadillo_bits/eop_core_bones.hpp @@ -98,6 +98,7 @@ class eop_ceil : public eop_core , public eo class eop_round : public eop_core , public eop_use_mp_false {}; class eop_trunc : public eop_core , public eop_use_mp_false {}; class eop_sign : public eop_core , public eop_use_mp_false {}; +class eop_cbrt : public eop_core , public eop_use_mp_true {}; class eop_erf : public eop_core , public eop_use_mp_true {}; class eop_erfc : public eop_core , public eop_use_mp_true {}; class eop_lgamma : public eop_core , public eop_use_mp_true {}; diff --git a/inst/include/armadillo_bits/eop_core_meat.hpp b/inst/include/armadillo_bits/eop_core_meat.hpp index 4bc0c7fe..330052f6 100644 --- a/inst/include/armadillo_bits/eop_core_meat.hpp +++ b/inst/include/armadillo_bits/eop_core_meat.hpp @@ -1142,6 +1142,9 @@ eop_core::process(const eT val, const eT ) { return eop_ template<> template arma_inline eT eop_core::process(const eT val, const eT ) { return arma_sign(val); } +template<> template arma_inline eT +eop_core::process(const eT val, const eT ) { return eop_aux::cbrt(val); } + template<> template arma_inline eT eop_core::process(const eT val, const eT ) { return eop_aux::erf(val); } diff --git a/inst/include/armadillo_bits/fn_conv_to.hpp b/inst/include/armadillo_bits/fn_conv_to.hpp index dbfc7fe2..4e2ed593 100644 --- a/inst/include/armadillo_bits/fn_conv_to.hpp +++ b/inst/include/armadillo_bits/fn_conv_to.hpp @@ -147,10 +147,15 @@ class conv_to< Mat > template inline static Mat from(const Base& in, const typename arma_cx_only::result* junk = nullptr); - template - inline static Mat from(const SpBase& in); + // + template + inline static Mat from(const SpBase& in, const typename arma_not_cx::result* junk = nullptr); + + template + inline static Mat from(const SpBase& in, const typename arma_cx_only::result* junk = nullptr); + // template inline static Mat from(const std::vector& in, const typename arma_not_cx::result* junk = nullptr); @@ -206,15 +211,59 @@ conv_to< Mat >::from(const Base& in, const typename arma_cx_o template -template +template arma_warn_unused inline Mat -conv_to< Mat >::from(const SpBase& in) +conv_to< Mat >::from(const SpBase& in, const typename arma_not_cx::result* junk) { arma_extra_debug_sigprint(); + arma_ignore(junk); + + const unwrap_spmat U(in.get_ref()); + const SpMat& X = U.M; + + Mat out(X.n_rows, X.n_cols, arma_zeros_indicator()); + + podarray tmp(X.n_nonzero); + + arrayops::convert( tmp.memptr(), X.values, X.n_nonzero ); - return Mat(in.get_ref()); + typename SpMat::const_iterator it = X.begin(); + typename SpMat::const_iterator it_end = X.end(); + + for(uword count=0; it != it_end; ++it, ++count) { out.at(it.row(), it.col()) = tmp[count]; } + + return out; + } + + + +template +template +arma_warn_unused +inline +Mat +conv_to< Mat >::from(const SpBase& in, const typename arma_cx_only::result* junk) + { + arma_extra_debug_sigprint(); + arma_ignore(junk); + + const unwrap_spmat U(in.get_ref()); + const SpMat& X = U.M; + + Mat out(X.n_rows, X.n_cols, arma_zeros_indicator()); + + podarray tmp(X.n_nonzero); + + arrayops::convert_cx( tmp.memptr(), X.values, X.n_nonzero ); + + typename SpMat::const_iterator it = X.begin(); + typename SpMat::const_iterator it_end = X.end(); + + for(uword count=0; it != it_end; ++it, ++count) { out.at(it.row(), it.col()) = tmp[count]; } + + return out; } @@ -279,7 +328,7 @@ class conv_to< Row > template inline static Row from(const Base& in, const typename arma_cx_only::result* junk = nullptr); - + // template inline static Row from(const std::vector& in, const typename arma_not_cx::result* junk = nullptr); @@ -398,7 +447,7 @@ class conv_to< Col > template inline static Col from(const Base& in, const typename arma_cx_only::result* junk = nullptr); - + // template inline static Col from(const std::vector& in, const typename arma_not_cx::result* junk = nullptr); @@ -517,8 +566,13 @@ class conv_to< SpMat > template inline static SpMat from(const SpBase& in, const typename arma_cx_only::result* junk = nullptr); - template - inline static SpMat from(const Base& in); + // + + template + inline static SpMat from(const Base& in, const typename arma_not_cx::result* junk = nullptr); + + template + inline static SpMat from(const Base& in, const typename arma_cx_only::result* junk = nullptr); }; @@ -572,15 +626,177 @@ conv_to< SpMat >::from(const SpBase& in, const typename arma_ template -template +template +arma_warn_unused +inline +SpMat +conv_to< SpMat >::from(const Base& in, const typename arma_not_cx::result* junk) + { + arma_extra_debug_sigprint(); + arma_ignore(junk); + + SpMat out; + + const quasi_unwrap U(in.get_ref()); + const Mat& X = U.M; + + if(is_same_type::yes) + { + const Mat& Y = reinterpret_cast&>(X); + + SpMat tmp(Y); + + out.steal_mem(tmp); + } + else + { + const uword X_n_rows = X.n_rows; + const uword X_n_cols = X.n_cols; + const uword X_n_elem = X.n_elem; + + const in_eT* X_mem = X.memptr(); + + uword X_nnz = 0; + + for(uword i=0; i < X_n_elem; ++i) { X_nnz += (X_mem[i] != in_eT(0)) ? uword(1) : uword(0); } + + podarray< in_eT> X_nonzeros(X_nnz); + podarray Y_nonzeros(X_nnz); + + for(uword i=0,count=0; i < X_n_elem; ++i) + { + const in_eT X_val = X_mem[i]; + + if(X_val != in_eT(0)) { X_nonzeros[count] = X_val; ++count; } + } + + arrayops::convert( Y_nonzeros.memptr(), X_nonzeros.memptr(), X_nnz ); + + if(X_nnz == 0) + { + out.set_size(X_n_rows, X.n_cols); + } + else + { + SpMat tmp(arma_reserve_indicator(), X_n_rows, X_n_cols, X_nnz); + + uword count = 0; + + for(uword c=0; c < X_n_cols; ++c) + for(uword r=0; r < X_n_rows; ++r) + { + const in_eT X_val = (*X_mem); ++X_mem; + + if(X_val != in_eT(0)) + { + access::rw(tmp.values[count]) = Y_nonzeros[count]; + access::rw(tmp.row_indices[count]) = r; + access::rw(tmp.col_ptrs[c + 1])++; + ++count; + } + } + + // Sum column counts to be column pointers. + for(uword c=1; c <= tmp.n_cols; ++c) + { + access::rw(tmp.col_ptrs[c]) += tmp.col_ptrs[c - 1]; + } + + tmp.remove_zeros(); // in case conversion resulted in an element equal to zero + + out.steal_mem(tmp); + } + } + + return out; + } + + + +template +template arma_warn_unused inline SpMat -conv_to< SpMat >::from(const Base& in) +conv_to< SpMat >::from(const Base& in, const typename arma_cx_only::result* junk) { arma_extra_debug_sigprint(); + arma_ignore(junk); - return SpMat(in.get_ref()); + SpMat out; + + const quasi_unwrap U(in.get_ref()); + const Mat& X = U.M; + + if(is_same_type::yes) + { + const Mat& Y = reinterpret_cast&>(X); + + SpMat tmp(Y); + + out.steal_mem(tmp); + } + else + { + const uword X_n_rows = X.n_rows; + const uword X_n_cols = X.n_cols; + const uword X_n_elem = X.n_elem; + + const in_eT* X_mem = X.memptr(); + + uword X_nnz = 0; + + for(uword i=0; i < X_n_elem; ++i) { X_nnz += (X_mem[i] != in_eT(0)) ? uword(1) : uword(0); } + + podarray< in_eT> X_nonzeros(X_nnz); + podarray Y_nonzeros(X_nnz); + + for(uword i=0,count=0; i < X_n_elem; ++i) + { + const in_eT X_val = X_mem[i]; + + if(X_val != in_eT(0)) { X_nonzeros[count] = X_val; ++count; } + } + + arrayops::convert_cx( Y_nonzeros.memptr(), X_nonzeros.memptr(), X_nnz ); + + if(X_nnz == 0) + { + out.set_size(X_n_rows, X.n_cols); + } + else + { + SpMat tmp(arma_reserve_indicator(), X_n_rows, X_n_cols, X_nnz); + + uword count = 0; + + for(uword c=0; c < X_n_cols; ++c) + for(uword r=0; r < X_n_rows; ++r) + { + const in_eT X_val = (*X_mem); ++X_mem; + + if(X_val != in_eT(0)) + { + access::rw(tmp.values[count]) = Y_nonzeros[count]; + access::rw(tmp.row_indices[count]) = r; + access::rw(tmp.col_ptrs[c + 1])++; + ++count; + } + } + + // Sum column counts to be column pointers. + for(uword c=1; c <= tmp.n_cols; ++c) + { + access::rw(tmp.col_ptrs[c]) += tmp.col_ptrs[c - 1]; + } + + tmp.remove_zeros(); // in case conversion resulted in an element equal to zero + + out.steal_mem(tmp); + } + } + + return out; } diff --git a/inst/include/armadillo_bits/fn_elem.hpp b/inst/include/armadillo_bits/fn_elem.hpp index 917537f4..f00f4fdb 100644 --- a/inst/include/armadillo_bits/fn_elem.hpp +++ b/inst/include/armadillo_bits/fn_elem.hpp @@ -677,6 +677,48 @@ sqrt(const SpBase& A) +// +// cbrt + +template +arma_warn_unused +arma_inline +typename enable_if2< (is_arma_type::value && is_cx::no), const eOp >::result +cbrt(const T1& A) + { + arma_extra_debug_sigprint(); + + return eOp(A); + } + + + +template +arma_warn_unused +arma_inline +typename enable_if2< is_cx::no, const eOpCube >::result +cbrt(const BaseCube& A) + { + arma_extra_debug_sigprint(); + + return eOpCube(A.get_ref()); + } + + + +template +arma_warn_unused +arma_inline +typename enable_if2< is_cx::no, const SpOp >::result +cbrt(const SpBase& A) + { + arma_extra_debug_sigprint(); + + return SpOp(A.get_ref()); + } + + + // // conj diff --git a/inst/include/armadillo_bits/fn_shift.hpp b/inst/include/armadillo_bits/fn_shift.hpp index d3de6a7c..879f5dca 100644 --- a/inst/include/armadillo_bits/fn_shift.hpp +++ b/inst/include/armadillo_bits/fn_shift.hpp @@ -115,4 +115,39 @@ shift +// + + + +template +arma_warn_unused +inline +SpMat +shift + ( + const SpBase& expr, + const sword N, + const uword dim = 0 + ) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + arma_debug_check( (dim > 1), "shift(): parameter 'dim' must be 0 or 1" ); + + const uword len = (N < 0) ? uword(-N) : uword(N); + const uword neg = (N < 0) ? uword( 1) : uword(0); + + unwrap_spmat U(expr.get_ref()); + + SpMat out; + + spop_shift::apply_noalias(out, U.M, len, neg, dim); + + return out; + } + + + //! @} diff --git a/inst/include/armadillo_bits/op_pinv_meat.hpp b/inst/include/armadillo_bits/op_pinv_meat.hpp index 326a0bef..f4dac14f 100644 --- a/inst/include/armadillo_bits/op_pinv_meat.hpp +++ b/inst/include/armadillo_bits/op_pinv_meat.hpp @@ -117,12 +117,17 @@ op_pinv::apply_direct(Mat& out, const Base::eval(expr.get_ref()); - sym_helper::analyse_matrix(is_approx_sym, is_approx_sympd, A); - - do_sym = ((is_cx::no) ? (is_approx_sym) : (is_approx_sym && is_approx_sympd)); + if(do_sym == false) + { + bool is_approx_sym = false; + bool is_approx_sympd = false; + + sym_helper::analyse_matrix(is_approx_sym, is_approx_sympd, A); + + do_sym = ((is_cx::no) ? (is_approx_sym) : (is_approx_sym && is_approx_sympd)); + } } if(do_sym) diff --git a/inst/include/armadillo_bits/op_rank_meat.hpp b/inst/include/armadillo_bits/op_rank_meat.hpp index ef00dd39..91c413e9 100644 --- a/inst/include/armadillo_bits/op_rank_meat.hpp +++ b/inst/include/armadillo_bits/op_rank_meat.hpp @@ -44,14 +44,21 @@ op_rank::apply(uword& out, const Base& expr, const ty bool do_sym = false; - if((arma_config::optimise_sym) && (auxlib::crippled_lapack(A) == false) && (A.n_rows >= (is_cx::yes ? uword(64) : uword(128)))) + const bool is_sym_size_ok = (A.n_rows == A.n_cols) && (A.n_rows > (is_cx::yes ? uword(20) : uword(40))); // for consistency with op_pinv + + if( (is_sym_size_ok) && (arma_config::optimise_sym) && (auxlib::crippled_lapack(A) == false) ) { - bool is_approx_sym = false; - bool is_approx_sympd = false; - - sym_helper::analyse_matrix(is_approx_sym, is_approx_sympd, A); + do_sym = is_sym_expr::eval(expr.get_ref()); - do_sym = (is_cx::no) ? (is_approx_sym) : (is_approx_sym && is_approx_sympd); + if(do_sym == false) + { + bool is_approx_sym = false; + bool is_approx_sympd = false; + + sym_helper::analyse_matrix(is_approx_sym, is_approx_sympd, A); + + do_sym = (is_cx::no) ? (is_approx_sym) : (is_approx_sym && is_approx_sympd); + } } if(do_sym) diff --git a/inst/include/armadillo_bits/op_shift_meat.hpp b/inst/include/armadillo_bits/op_shift_meat.hpp index b369b5d3..216290b5 100644 --- a/inst/include/armadillo_bits/op_shift_meat.hpp +++ b/inst/include/armadillo_bits/op_shift_meat.hpp @@ -64,6 +64,8 @@ op_shift::apply_noalias(Mat& out, const Mat& X, const uword len, const u arma_debug_check_bounds( ((dim == 0) && (len >= X.n_rows)), "shift(): shift amount out of bounds" ); arma_debug_check_bounds( ((dim == 1) && (len >= X.n_cols)), "shift(): shift amount out of bounds" ); + if(len == 0) { out = X; return; } + out.copy_size(X); const uword X_n_rows = X.n_rows; diff --git a/inst/include/armadillo_bits/sp_auxlib_meat.hpp b/inst/include/armadillo_bits/sp_auxlib_meat.hpp index dbfdf2d7..97df7e8c 100644 --- a/inst/include/armadillo_bits/sp_auxlib_meat.hpp +++ b/inst/include/armadillo_bits/sp_auxlib_meat.hpp @@ -1984,49 +1984,14 @@ sp_auxlib::run_aupd_plain case 1: { // We need to calculate the matrix-vector multiplication y = OP * x - // where x is of length n and starts at workd(ipntr(0)), and y is of - // length n and starts at workd(ipntr(1)). - - // // OLD METHOD - // - // // operator*(sp_mat, vec) doesn't properly put the result into the - // // right place so we'll just reimplement it here for now... - // - // // Set the output to point at the right memory. We have to subtract - // // one from FORTRAN pointers... - // Col out(workd.memptr() + ipntr(1) - 1, n, false /* don't copy */); - // // Set the input to point at the right memory. - // Col in(workd.memptr() + ipntr(0) - 1, n, false /* don't copy */); - // - // out.zeros(); - // - // T* out_mem = out.memptr(); - // const T* in_mem = in.memptr(); - // - // typename SpMat::const_iterator X_it = X.begin(); - // - // const uword X_nnz = X.n_nonzero; - // - // for(uword count=0; count < X_nnz; ++count, ++X_it) - // { - // const eT X_it_val = (*X_it); - // const uword X_it_row = X_it.row(); - // const uword X_it_col = X_it.col(); - // - // out_mem[X_it_row] += X_it_val * in_mem[X_it_col]; - // } - // - // // No need to modify memory further since it was all done in-place. - - - // NEW METHOD - // - // both operator*(rowvec, sp_mat) and operator*(sp_mat, colvec) can now write to an existing object + // where x is of length n and starts at workd(ipntr(0)), + // and y is of length n and starts at workd(ipntr(1)). + // We have to subtract one from FORTRAN pointers. Row out(workd.memptr() + ipntr(1) - 1, n, false, true); Row in(workd.memptr() + ipntr(0) - 1, n, false, true); - out = in * Xst; + out = in * Xst; // using transposed version break; } @@ -2259,14 +2224,12 @@ sp_auxlib::run_aupd_shiftinvert case 1: { // We need to calculate the matrix-vector multiplication y = OP * x - // where x is of length n and starts at workd(ipntr(0)), and y is of - // length n and starts at workd(ipntr(1)). + // where x is of length n and starts at workd(ipntr(0)), + // and y is of length n and starts at workd(ipntr(1)). + // We have to subtract one from FORTRAN pointers. - // Set the output to point at the right memory. We have to subtract - // one from FORTRAN pointers... - Col out(workd.memptr() + ipntr(1) - 1, n, false /* don't copy */); - // Set the input to point at the right memory. - Col in(workd.memptr() + ipntr(0) - 1, n, false /* don't copy */); + Col out(workd.memptr() + ipntr(1) - 1, n, false, true); + Col in(workd.memptr() + ipntr(0) - 1, n, false, true); // Consider getting the LU factorization from ZGSTRF, and then // solve the system L*U*out = in (possibly with permutation matrix?) diff --git a/inst/include/armadillo_bits/spop_misc_bones.hpp b/inst/include/armadillo_bits/spop_misc_bones.hpp index 42117f9b..ed62df20 100644 --- a/inst/include/armadillo_bits/spop_misc_bones.hpp +++ b/inst/include/armadillo_bits/spop_misc_bones.hpp @@ -64,6 +64,17 @@ class spop_sqrt +class spop_cbrt + : public traits_op_passthru + { + public: + + template + inline static void apply(SpMat& out, const SpOp& in); + }; + + + class spop_abs : public traits_op_passthru { diff --git a/inst/include/armadillo_bits/spop_misc_meat.hpp b/inst/include/armadillo_bits/spop_misc_meat.hpp index 1ef51cc4..78b8636a 100644 --- a/inst/include/armadillo_bits/spop_misc_meat.hpp +++ b/inst/include/armadillo_bits/spop_misc_meat.hpp @@ -146,6 +146,29 @@ spop_sqrt::apply(SpMat& out, const SpOp& i +namespace priv + { + struct functor_cbrt + { + template + arma_inline eT operator()(const eT val) const { return eop_aux::cbrt(val); } + }; + } + + + +template +inline +void +spop_cbrt::apply(SpMat& out, const SpOp& in) + { + arma_extra_debug_sigprint(); + + out.init_xform(in.m, priv::functor_cbrt()); + } + + + namespace priv { struct functor_abs @@ -328,8 +351,8 @@ spop_repelem::apply(SpMat& out, const SpOp 0) && (out_n_cols > 0) && (out_nnz > 0) ) { - umat locs(2, out_nnz, arma_nozeros_indicator()); - Col vals( out_nnz, arma_nozeros_indicator()); + Mat locs(2, out_nnz, arma_nozeros_indicator()); + Col vals( out_nnz, arma_nozeros_indicator()); uword* locs_mem = locs.memptr(); eT* vals_mem = vals.memptr(); diff --git a/inst/include/armadillo_bits/spop_shift_bones.hpp b/inst/include/armadillo_bits/spop_shift_bones.hpp new file mode 100644 index 00000000..4475c406 --- /dev/null +++ b/inst/include/armadillo_bits/spop_shift_bones.hpp @@ -0,0 +1,35 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// ------------------------------------------------------------------------ + + + +//! \addtogroup spop_shift +//! @{ + + + +class spop_shift + : public traits_op_default + { + public: + + template inline static void apply_noalias(SpMat& out, const SpMat& X, const uword len, const uword neg, const uword dim); + }; + + + +//! @} diff --git a/inst/include/armadillo_bits/spop_shift_meat.hpp b/inst/include/armadillo_bits/spop_shift_meat.hpp new file mode 100644 index 00000000..1fd14aa0 --- /dev/null +++ b/inst/include/armadillo_bits/spop_shift_meat.hpp @@ -0,0 +1,93 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// ------------------------------------------------------------------------ + + + +//! \addtogroup spop_shift +//! @{ + + + +template +inline +void +spop_shift::apply_noalias(SpMat& out, const SpMat& X, const uword len, const uword neg, const uword dim) + { + arma_extra_debug_sigprint(); + + arma_debug_check_bounds( ((dim == 0) && (len >= X.n_rows)), "shift(): shift amount out of bounds" ); + arma_debug_check_bounds( ((dim == 1) && (len >= X.n_cols)), "shift(): shift amount out of bounds" ); + + if(X.n_nonzero == 0) { out.zeros(X.n_rows, X.n_cols); return; } + + if(len == 0) { out = X; return; } + + typename SpMat::const_iterator it = X.begin(); + typename SpMat::const_iterator it_end = X.end(); + + Mat locs(2, X.n_nonzero, arma_nozeros_indicator()); + + uword* locs_mem = locs.memptr(); + + if(dim == 0) + { + const uword X_row_threshold = X.n_rows - len; + + for(; it != it_end; ++it) + { + const uword X_row = it.row(); + uword Y_row = 0; + + if(neg == 0) { Y_row = (X_row < X_row_threshold) ? (X_row + len) : (X_row - X_row_threshold); } + else if(neg == 1) { Y_row = (X_row >= len ) ? (X_row - len) : (X_row + X_row_threshold); } + + locs_mem[0] = Y_row; + locs_mem[1] = it.col(); + + locs_mem += 2; + } + } + else + if(dim == 1) + { + const uword X_col_threshold = X.n_cols - len; + + for(; it != it_end; ++it) + { + const uword X_col = it.col(); + uword Y_col = 0; + + if(neg == 0) { Y_col = (X_col < X_col_threshold) ? (X_col + len) : (X_col - X_col_threshold); } + else if(neg == 1) { Y_col = (X_col >= len ) ? (X_col - len) : (X_col + X_col_threshold); } + + locs_mem[0] = it.row(); + locs_mem[1] = Y_col; + + locs_mem += 2; + } + } + + const Col vals(const_cast(X.values), X.n_nonzero, false); + + SpMat Y(locs, vals, X.n_rows, X.n_cols, true, false); + + out.steal_mem(Y); + } + + + +//! @} diff --git a/inst/include/armadillo_bits/traits.hpp b/inst/include/armadillo_bits/traits.hpp index bde22009..e778fef8 100644 --- a/inst/include/armadillo_bits/traits.hpp +++ b/inst/include/armadillo_bits/traits.hpp @@ -1312,4 +1312,43 @@ struct has_nested_glue_traits }; + + +template +struct is_sym_expr + { + static constexpr bool eval(const T1&) { return false; } + }; + +template +struct is_sym_expr< Glue< Mat, Op, op_htrans>, glue_times > > + { + static + arma_inline + bool + eval(const Glue< Mat, Op, op_htrans>, glue_times >& expr) + { + const Mat& X = expr.A; + const Mat& Y = expr.B.m; + + return (&X == &Y); + } + }; + +template +struct is_sym_expr< Glue< Op, op_htrans>, Mat, glue_times > > + { + static + arma_inline + bool + eval(const Glue< Op, op_htrans>, Mat, glue_times >& expr) + { + const Mat& X = expr.A.m; + const Mat& Y = expr.B; + + return (&X == &Y); + } + }; + + //! @}