From 1795d9405d167a2fbad69866397d78bc69aa6bb0 Mon Sep 17 00:00:00 2001 From: conrad Date: Wed, 16 Mar 2022 14:24:37 +1000 Subject: [PATCH] move classes --- include/armadillo | 2 - include/armadillo_bits/fn_inv.hpp | 52 ++--- include/armadillo_bits/op_inv_gen_bones.hpp | 11 + include/armadillo_bits/op_inv_gen_meat.hpp | 94 ++++++++ include/armadillo_bits/op_inv_rcond_bones.hpp | 38 ---- include/armadillo_bits/op_inv_rcond_meat.hpp | 207 ------------------ include/armadillo_bits/op_inv_spd_bones.hpp | 11 + include/armadillo_bits/op_inv_spd_meat.hpp | 95 ++++++++ 8 files changed, 237 insertions(+), 273 deletions(-) delete mode 100644 include/armadillo_bits/op_inv_rcond_bones.hpp delete mode 100644 include/armadillo_bits/op_inv_rcond_meat.hpp diff --git a/include/armadillo b/include/armadillo index cae5573f..35df03de 100644 --- a/include/armadillo +++ b/include/armadillo @@ -214,7 +214,6 @@ namespace arma #include "armadillo_bits/op_log_det_bones.hpp" #include "armadillo_bits/op_inv_gen_bones.hpp" #include "armadillo_bits/op_inv_spd_bones.hpp" - #include "armadillo_bits/op_inv_rcond_bones.hpp" #include "armadillo_bits/op_htrans_bones.hpp" #include "armadillo_bits/op_max_bones.hpp" #include "armadillo_bits/op_min_bones.hpp" @@ -646,7 +645,6 @@ namespace arma #include "armadillo_bits/op_log_det_meat.hpp" #include "armadillo_bits/op_inv_gen_meat.hpp" #include "armadillo_bits/op_inv_spd_meat.hpp" - #include "armadillo_bits/op_inv_rcond_meat.hpp" #include "armadillo_bits/op_htrans_meat.hpp" #include "armadillo_bits/op_max_meat.hpp" #include "armadillo_bits/op_index_max_meat.hpp" diff --git a/include/armadillo_bits/fn_inv.hpp b/include/armadillo_bits/fn_inv.hpp index 6ff78d5c..4cb2f74c 100644 --- a/include/armadillo_bits/fn_inv.hpp +++ b/include/armadillo_bits/fn_inv.hpp @@ -61,31 +61,6 @@ inv -template -inline -typename enable_if2< is_supported_blas_type::value, bool >::result -inv - ( - Mat& out_inv, - typename T1::pod_type& out_rcond, - const Base& X - ) - { - arma_extra_debug_sigprint(); - - const bool status = op_inv_rcond::apply_direct_gen(out_inv, out_rcond, X.get_ref()); - - if(status == false) - { - out_inv.soft_reset(); - arma_debug_warn_level(3, "inv(): matrix is singular"); - } - - return status; - } - - - template arma_warn_unused arma_inline @@ -128,6 +103,31 @@ inv +template +inline +typename enable_if2< is_supported_blas_type::value, bool >::result +inv + ( + Mat& out_inv, + typename T1::pod_type& out_rcond, + const Base& X + ) + { + arma_extra_debug_sigprint(); + + const bool status = op_inv_gen_rcond::apply_direct(out_inv, out_rcond, X.get_ref()); + + if(status == false) + { + out_inv.soft_reset(); + arma_debug_warn_level(3, "inv(): matrix is singular"); + } + + return status; + } + + + // @@ -226,7 +226,7 @@ inv_sympd { arma_extra_debug_sigprint(); - const bool status = op_inv_rcond::apply_direct_spd(out_inv, out_rcond, X.get_ref()); + const bool status = op_inv_spd_rcond::apply_direct(out_inv, out_rcond, X.get_ref()); if(status == false) { diff --git a/include/armadillo_bits/op_inv_gen_bones.hpp b/include/armadillo_bits/op_inv_gen_bones.hpp index 4e817821..42d2345b 100644 --- a/include/armadillo_bits/op_inv_gen_bones.hpp +++ b/include/armadillo_bits/op_inv_gen_bones.hpp @@ -60,6 +60,17 @@ class op_inv_gen_full +class op_inv_gen_rcond + : public traits_op_default + { + public: + + template + inline static bool apply_direct(Mat& out_inv, typename T1::pod_type& out_rcond, const Base& expr); + }; + + + namespace inv_opts { struct opts diff --git a/include/armadillo_bits/op_inv_gen_meat.hpp b/include/armadillo_bits/op_inv_gen_meat.hpp index 7c84115a..c88d947c 100644 --- a/include/armadillo_bits/op_inv_gen_meat.hpp +++ b/include/armadillo_bits/op_inv_gen_meat.hpp @@ -295,4 +295,98 @@ op_inv_gen_full::apply_tiny(Mat& X) +template +inline +bool +op_inv_gen_rcond::apply_direct(Mat& out, typename T1::pod_type& out_rcond, const Base& expr) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + typedef typename T1::pod_type T; + + out = expr.get_ref(); + out_rcond = T(0); + + arma_debug_check( (out.is_square() == false), "inv(): given matrix must be square sized" ); + + const uword N = out.n_rows; + + if(is_op_diagmat::value || out.is_diagmat()) + { + arma_extra_debug_print("op_inv_gen_rcond: detected diagonal matrix"); + + eT* colmem = out.memptr(); + + T max_abs_src_val = T(0); + T max_abs_inv_val = T(0); + + for(uword i=0; i max_abs_src_val) ? abs_src_val : max_abs_src_val; + max_abs_inv_val = (abs_inv_val > max_abs_inv_val) ? abs_inv_val : max_abs_inv_val; + + colmem += N; + } + + out_rcond = T(1) / (max_abs_src_val * max_abs_inv_val); + + return true; + } + + const strip_trimat strip(expr.get_ref()); + + const bool is_triu_expr = strip.do_triu; + const bool is_tril_expr = strip.do_tril; + + const bool is_triu_mat = (is_triu_expr || is_tril_expr) ? false : ( trimat_helper::is_triu(out)); + const bool is_tril_mat = (is_triu_expr || is_tril_expr) ? false : ((is_triu_mat) ? false : trimat_helper::is_tril(out)); + + if(is_triu_expr || is_tril_expr || is_triu_mat || is_tril_mat) + { + return auxlib::inv_tr_rcond(out, out_rcond, ((is_triu_expr || is_triu_mat) ? uword(0) : uword(1))); + } + + #if defined(ARMA_OPTIMISE_SYMPD) + const bool try_sympd = (auxlib::crippled_lapack(out)) ? false : sympd_helper::guess_sympd(out); + #else + const bool try_sympd = false; + #endif + + if(try_sympd) + { + arma_extra_debug_print("op_inv_gen_rcond: attempting sympd optimisation"); + + Mat tmp = out; + + bool sympd_state = false; + + const bool status = auxlib::inv_sympd_rcond(tmp, sympd_state, out_rcond, T(-1)); + + if(status) { out.steal_mem(tmp); return true; } + + if((status == false) && (sympd_state == true)) { return false; } + + arma_extra_debug_print("op_inv_gen_rcond: sympd optimisation failed"); + + // fallthrough if optimisation failed + } + + return auxlib::inv_rcond(out, out_rcond); + } + + + //! @} diff --git a/include/armadillo_bits/op_inv_rcond_bones.hpp b/include/armadillo_bits/op_inv_rcond_bones.hpp deleted file mode 100644 index 0b0b11e1..00000000 --- a/include/armadillo_bits/op_inv_rcond_bones.hpp +++ /dev/null @@ -1,38 +0,0 @@ -// 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 op_inv_rcond -//! @{ - - - -class op_inv_rcond - : public traits_op_default - { - public: - - template - inline static bool apply_direct_gen(Mat& out_inv, typename T1::pod_type& out_rcond, const Base& expr); - - template - inline static bool apply_direct_spd(Mat& out_inv, typename T1::pod_type& out_rcond, const Base& expr); - }; - - - -//! @} diff --git a/include/armadillo_bits/op_inv_rcond_meat.hpp b/include/armadillo_bits/op_inv_rcond_meat.hpp deleted file mode 100644 index 1ce87095..00000000 --- a/include/armadillo_bits/op_inv_rcond_meat.hpp +++ /dev/null @@ -1,207 +0,0 @@ -// 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 op_inv_rcond -//! @{ - - - -template -inline -bool -op_inv_rcond::apply_direct_gen(Mat& out, typename T1::pod_type& out_rcond, const Base& expr) - { - arma_extra_debug_sigprint(); - - typedef typename T1::elem_type eT; - typedef typename T1::pod_type T; - - out = expr.get_ref(); - out_rcond = T(0); - - arma_debug_check( (out.is_square() == false), "inv(): given matrix must be square sized" ); - - const uword N = out.n_rows; - - if(is_op_diagmat::value || out.is_diagmat()) - { - arma_extra_debug_print("op_inv_rcond: detected diagonal matrix"); - - eT* colmem = out.memptr(); - - T max_abs_src_val = T(0); - T max_abs_inv_val = T(0); - - for(uword i=0; i max_abs_src_val) ? abs_src_val : max_abs_src_val; - max_abs_inv_val = (abs_inv_val > max_abs_inv_val) ? abs_inv_val : max_abs_inv_val; - - colmem += N; - } - - out_rcond = T(1) / (max_abs_src_val * max_abs_inv_val); - - return true; - } - - const strip_trimat strip(expr.get_ref()); - - const bool is_triu_expr = strip.do_triu; - const bool is_tril_expr = strip.do_tril; - - const bool is_triu_mat = (is_triu_expr || is_tril_expr) ? false : ( trimat_helper::is_triu(out)); - const bool is_tril_mat = (is_triu_expr || is_tril_expr) ? false : ((is_triu_mat) ? false : trimat_helper::is_tril(out)); - - if(is_triu_expr || is_tril_expr || is_triu_mat || is_tril_mat) - { - return auxlib::inv_tr_rcond(out, out_rcond, ((is_triu_expr || is_triu_mat) ? uword(0) : uword(1))); - } - - #if defined(ARMA_OPTIMISE_SYMPD) - const bool try_sympd = (auxlib::crippled_lapack(out)) ? false : sympd_helper::guess_sympd(out); - #else - const bool try_sympd = false; - #endif - - if(try_sympd) - { - arma_extra_debug_print("op_inv_rcond: attempting sympd optimisation"); - - Mat tmp = out; - - bool sympd_state = false; - - const bool status = auxlib::inv_sympd_rcond(tmp, sympd_state, out_rcond, T(-1)); - - if(status) { out.steal_mem(tmp); return true; } - - if((status == false) && (sympd_state == true)) { return false; } - - arma_extra_debug_print("op_inv_rcond: sympd optimisation failed"); - - // fallthrough if optimisation failed - } - - return auxlib::inv_rcond(out, out_rcond); - } - - - -template -inline -bool -op_inv_rcond::apply_direct_spd(Mat& out, typename T1::pod_type& out_rcond, const Base& expr) - { - arma_extra_debug_sigprint(); - - typedef typename T1::elem_type eT; - typedef typename T1::pod_type T; - - if(auxlib::crippled_lapack(out)) - { - out.soft_reset(); - out_rcond = T(0); - arma_stop_runtime_error("inv_sympd(): not available for complex matrices due to presence of incomplete LAPACK"); - return false; - } - - out = expr.get_ref(); - out_rcond = T(0); - - arma_debug_check( (out.is_square() == false), "inv_sympd(): given matrix must be square sized" ); - - if((arma_config::debug) && (auxlib::rudimentary_sym_check(out) == false)) - { - if(is_cx::no ) { arma_debug_warn_level(1, "inv_sympd(): given matrix is not symmetric"); } - if(is_cx::yes) { arma_debug_warn_level(1, "inv_sympd(): given matrix is not hermitian"); } - } - - const uword N = out.n_rows; - - if(is_cx::yes) - { - arma_extra_debug_print("op_inv_rcond: checking imaginary components of diagonal elements"); - - const T tol = T(100) * std::numeric_limits::epsilon(); // allow some leeway - - const eT* colmem = out.memptr(); - - for(uword i=0; i tol) { return false; } - - colmem += N; - } - } - - if(is_op_diagmat::value || out.is_diagmat()) - { - arma_extra_debug_print("op_inv_rcond: detected diagonal matrix"); - - eT* colmem = out.memptr(); - - T max_abs_src_val = T(0); - T max_abs_inv_val = T(0); - - for(uword i=0; i max_abs_src_val) ? abs_src_val : max_abs_src_val; - max_abs_inv_val = (abs_inv_val > max_abs_inv_val) ? abs_inv_val : max_abs_inv_val; - - colmem += N; - } - - out_rcond = T(1) / (max_abs_src_val * max_abs_inv_val); - - return true; - } - - return auxlib::inv_sympd_rcond(out, out_rcond, T(-1)); - } - - - -//! @} diff --git a/include/armadillo_bits/op_inv_spd_bones.hpp b/include/armadillo_bits/op_inv_spd_bones.hpp index ec2fdc6b..ac716922 100644 --- a/include/armadillo_bits/op_inv_spd_bones.hpp +++ b/include/armadillo_bits/op_inv_spd_bones.hpp @@ -52,4 +52,15 @@ class op_inv_spd_full +class op_inv_spd_rcond + : public traits_op_default + { + public: + + template + inline static bool apply_direct(Mat& out_inv, typename T1::pod_type& out_rcond, const Base& expr); + }; + + + //! @} diff --git a/include/armadillo_bits/op_inv_spd_meat.hpp b/include/armadillo_bits/op_inv_spd_meat.hpp index 535cb77f..f12349c2 100644 --- a/include/armadillo_bits/op_inv_spd_meat.hpp +++ b/include/armadillo_bits/op_inv_spd_meat.hpp @@ -228,4 +228,99 @@ op_inv_spd_full::apply_tiny(Mat& out) +// + + + +template +inline +bool +op_inv_spd_rcond::apply_direct(Mat& out, typename T1::pod_type& out_rcond, const Base& expr) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + typedef typename T1::pod_type T; + + if(auxlib::crippled_lapack(out)) + { + out.soft_reset(); + out_rcond = T(0); + arma_stop_runtime_error("inv_sympd(): not available for complex matrices due to presence of incomplete LAPACK"); + return false; + } + + out = expr.get_ref(); + out_rcond = T(0); + + arma_debug_check( (out.is_square() == false), "inv_sympd(): given matrix must be square sized" ); + + if((arma_config::debug) && (auxlib::rudimentary_sym_check(out) == false)) + { + if(is_cx::no ) { arma_debug_warn_level(1, "inv_sympd(): given matrix is not symmetric"); } + if(is_cx::yes) { arma_debug_warn_level(1, "inv_sympd(): given matrix is not hermitian"); } + } + + const uword N = out.n_rows; + + if(is_cx::yes) + { + arma_extra_debug_print("op_inv_spd_rcond: checking imaginary components of diagonal elements"); + + const T tol = T(100) * std::numeric_limits::epsilon(); // allow some leeway + + const eT* colmem = out.memptr(); + + for(uword i=0; i tol) { return false; } + + colmem += N; + } + } + + if(is_op_diagmat::value || out.is_diagmat()) + { + arma_extra_debug_print("op_inv_spd_rcond: detected diagonal matrix"); + + eT* colmem = out.memptr(); + + T max_abs_src_val = T(0); + T max_abs_inv_val = T(0); + + for(uword i=0; i max_abs_src_val) ? abs_src_val : max_abs_src_val; + max_abs_inv_val = (abs_inv_val > max_abs_inv_val) ? abs_inv_val : max_abs_inv_val; + + colmem += N; + } + + out_rcond = T(1) / (max_abs_src_val * max_abs_inv_val); + + return true; + } + + bool is_sympd_junk = false; + + return auxlib::inv_sympd_rcond(out, is_sympd_junk, out_rcond, T(-1)); + } + + + //! @}