diff --git a/include/armadillo b/include/armadillo index 325357e6..0edff89a 100644 --- a/include/armadillo +++ b/include/armadillo @@ -278,6 +278,7 @@ namespace arma #include "armadillo_bits/op_sp_plus_bones.hpp" #include "armadillo_bits/op_sp_minus_bones.hpp" #include "armadillo_bits/op_powmat_bones.hpp" + #include "armadillo_bits/op_rank_bones.hpp" #include "armadillo_bits/glue_times_bones.hpp" #include "armadillo_bits/glue_mixed_bones.hpp" @@ -707,6 +708,7 @@ namespace arma #include "armadillo_bits/op_sp_plus_meat.hpp" #include "armadillo_bits/op_sp_minus_meat.hpp" #include "armadillo_bits/op_powmat_meat.hpp" + #include "armadillo_bits/op_rank_meat.hpp" #include "armadillo_bits/glue_times_meat.hpp" #include "armadillo_bits/glue_mixed_meat.hpp" diff --git a/include/armadillo_bits/fn_rank.hpp b/include/armadillo_bits/fn_rank.hpp index f4453def..7701a049 100644 --- a/include/armadillo_bits/fn_rank.hpp +++ b/include/armadillo_bits/fn_rank.hpp @@ -21,177 +21,35 @@ -template -arma_warn_unused -inline -typename enable_if2< is_supported_blas_type::value, uword >::result -rank_diagmat - ( - const Mat& A, - typename get_pod_type::result tol = 0.0 - ) - { - arma_extra_debug_sigprint(); - - typedef typename get_pod_type::result T; - - if(A.is_empty()) { return uword(0); } - - const uword N = (std::min)(A.n_rows, A.n_cols); - - podarray diag_abs_vals(N); - - T max_abs_Aii = T(0); - - for(uword i=0; i max_abs_Aii) ? abs_Aii : max_abs_Aii; - } - - if(tol == T(0)) { tol = (std::max)(A.n_rows, A.n_cols) * max_abs_Aii * std::numeric_limits::epsilon(); } - - uword count = 0; - - for(uword i=0; i tol) ? uword(1) : uword(0); } - - return count; - } - - - template arma_warn_unused inline typename enable_if2< is_supported_blas_type::value, uword >::result -rank_sym - ( - const Base& X, - typename T1::pod_type tol = 0.0 - ) +rank(const Base& expr, const typename T1::pod_type tol = 0) { arma_extra_debug_sigprint(); - typedef typename T1::elem_type eT; - typedef typename T1::pod_type T; + uword out = uword(0); - Mat A(X.get_ref()); + const bool status = op_rank::apply(out, expr.get_ref(), tol); - if(A.is_square() == false) - { - arma_stop_runtime_error("rank_sym(): given matrix must be square sized"); - return uword(0); - } + if(status == false) { arma_stop_runtime_error("rank(): failed"); return uword(0); } - if(A.is_empty()) { return uword(0); } - - // if(is_op_diagmat::value || A.is_diagmat()) - // { - // arma_extra_debug_print("rank(): detected diagonal matrix"); - // // TODO: return rank_diagmat(A, tol); - // } - - Col v; - - const bool status = auxlib::eig_sym(v, A); - - if(status == false) - { - arma_stop_runtime_error("rank_sym(): eig_sym failed"); - return uword(0); - } - - const uword v_n_elem = v.n_elem; - const T* v_mem = v.memptr(); - - // set tolerance to default if it hasn't been specified - if( (tol == T(0)) && (v_n_elem > 0) ) - { - tol = (std::max)(A.n_rows, A.n_cols) * v_mem[v_n_elem-1] * std::numeric_limits::epsilon(); - } - - uword count = 0; - - for(uword i=0; i < v_n_elem; ++i) { count += (v_mem[i] > tol) ? uword(1) : uword(0); } - - return count; + return out; } template -arma_warn_unused inline -uword -rank - ( - const Base& X, - typename T1::pod_type tol = 0.0, - const typename arma_blas_type_only::result* junk = nullptr - ) +typename enable_if2< is_supported_blas_type::value, bool >::result +rank(uword& out, const Base& expr, const typename T1::pod_type tol = 0) { arma_extra_debug_sigprint(); - arma_ignore(junk); - typedef typename T1::elem_type eT; - typedef typename T1::pod_type T; + out = uword(0); - Mat A(X.get_ref()); - - if(A.is_empty()) { return uword(0); } - - if(is_op_diagmat::value || A.is_diagmat()) - { - arma_extra_debug_print("rank(): detected diagonal matrix"); - - return rank_diagmat(A, tol); - } - - if(is_cx::no && A.is_symmetric()) - { - arma_extra_debug_print("rank(): detected symmetric matrix"); - - return rank_sym(A, tol); - } - - if(is_cx::yes && A.is_hermitian()) - { - arma_extra_debug_print("rank(): detected hermitian matrix"); - - return rank_sym(A, tol); - } - - Col s; - - const bool status = auxlib::svd_dc(s, A); - - if(status == false) - { - arma_stop_runtime_error("rank(): svd failed"); - - return uword(0); - } - - const uword s_n_elem = s.n_elem; - const T* s_mem = s.memptr(); - - // set tolerance to default if it hasn't been specified - if( (tol == T(0)) && (s_n_elem > 0) ) - { - tol = (std::max)(A.n_rows, A.n_cols) * s_mem[0] * std::numeric_limits::epsilon(); - } - - uword count = 0; - - for(uword i=0; i < s_n_elem; ++i) { count += (s_mem[i] > tol) ? uword(1) : uword(0); } - - return count; + return op_rank::apply(out, expr.get_ref(), tol); } diff --git a/include/armadillo_bits/op_rank_bones.hpp b/include/armadillo_bits/op_rank_bones.hpp new file mode 100644 index 00000000..f0c4a072 --- /dev/null +++ b/include/armadillo_bits/op_rank_bones.hpp @@ -0,0 +1,41 @@ +// 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_rank +//! @{ + + + +class op_rank + : public traits_op_default + { + public: + + template inline static bool apply(uword& out, const Base& expr, const typename T1::pod_type tol); + + template inline static bool apply_gen(uword& out, Mat& A, typename get_pod_type::result tol); + + template inline static bool apply_sym(uword& out, Mat& A, typename get_pod_type::result tol); + + template inline static bool apply_diag(uword& out, Mat& A, typename get_pod_type::result tol); + }; + + + +//! @} diff --git a/include/armadillo_bits/op_rank_meat.hpp b/include/armadillo_bits/op_rank_meat.hpp new file mode 100644 index 00000000..428bf1d1 --- /dev/null +++ b/include/armadillo_bits/op_rank_meat.hpp @@ -0,0 +1,175 @@ +// 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_rank +//! @{ + + + +template +inline +bool +op_rank::apply(uword& out, const Base& expr, typename T1::pod_type tol) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + Mat A(expr.get_ref()); + + if(A.is_empty()) { out = uword(0); return true; } + + if(is_op_diagmat::value || A.is_diagmat()) + { + arma_extra_debug_print("op_rank::apply(): detected diagonal matrix"); + + return op_rank::apply_diag(out, A, tol); + } + + if(is_cx::no && A.is_symmetric()) + { + arma_extra_debug_print("op_rank::apply(): detected symmetric matrix"); + + return op_rank::apply_sym(out, A, tol); + } + + if(is_cx::yes && A.is_hermitian()) + { + arma_extra_debug_print("op_rank::apply(): detected hermitian matrix"); + + return op_rank::apply_sym(out, A, tol); + } + + return op_rank::apply_gen(out, A, tol); + } + + + +template +inline +bool +op_rank::apply_diag(uword& out, Mat& A, typename get_pod_type::result tol) + { + arma_extra_debug_sigprint(); + + typedef typename get_pod_type::result T; + + const uword N = (std::min)(A.n_rows, A.n_cols); + + podarray diag_abs_vals(N); + + T max_abs_Aii = T(0); + + for(uword i=0; i max_abs_Aii) ? abs_Aii : max_abs_Aii; + } + + // set tolerance to default if it hasn't been specified + if(tol == T(0)) { tol = (std::max)(A.n_rows, A.n_cols) * max_abs_Aii * std::numeric_limits::epsilon(); } + + uword count = 0; + + for(uword i=0; i tol) ? uword(1) : uword(0); } + + out = count; + + return true; + } + + + +template +inline +bool +op_rank::apply_sym(uword& out, Mat& A, typename get_pod_type::result tol) + { + arma_extra_debug_sigprint(); + + typedef typename get_pod_type::result T; + + if(A.is_square() == false) { out = uword(0); return false; } + + Col v; + + const bool status = auxlib::eig_sym(v, A); + + if(status == false) { out = uword(0); return false; } + + const uword v_n_elem = v.n_elem; + const T* v_mem = v.memptr(); + + if(v_n_elem == 0) { out = uword(0); return true; } + + // set tolerance to default if it hasn't been specified + if(tol == T(0)) { tol = (std::max)(A.n_rows, A.n_cols) * v_mem[v_n_elem-1] * std::numeric_limits::epsilon(); } + + uword count = 0; + + for(uword i=0; i < v_n_elem; ++i) { count += (v_mem[i] > tol) ? uword(1) : uword(0); } + + out = count; + + return true; + } + + + +template +inline +bool +op_rank::apply_gen(uword& out, Mat& A, typename get_pod_type::result tol) + { + arma_extra_debug_sigprint(); + + typedef typename get_pod_type::result T; + + Col s; + + const bool status = auxlib::svd_dc(s, A); + + if(status == false) { out = uword(0); return false; } + + const uword s_n_elem = s.n_elem; + const T* s_mem = s.memptr(); + + if(s_n_elem == 0) { out = uword(0); return true; } + + // set tolerance to default if it hasn't been specified + if(tol == T(0)) { tol = (std::max)(A.n_rows, A.n_cols) * s_mem[0] * std::numeric_limits::epsilon(); } + + uword count = 0; + + for(uword i=0; i < s_n_elem; ++i) { count += (s_mem[i] > tol) ? uword(1) : uword(0); } + + out = count; + + return true; + } + + + +//! @}