phynexis-deps/armadillo-code
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

refactor rank() into op_rank

Browse Source
This commit is contained in:
conrad
2021-12-01 12:05:42 +10:00
parent 93abe2557d
commit 3a12b86552
4 changed files with 227 additions and 151 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/armadillo
+2
View File
@@ -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"
include/armadillo_bits/fn_rank.hpp
+9 -151
View File
@@ -21,177 +21,35 @@
template<typename eT>
arma_warn_unused
inline
typename enable_if2< is_supported_blas_type<eT>::value, uword >::result
rank_diagmat
(
const Mat<eT>& A,
typename get_pod_type<eT>::result tol = 0.0
)
{
arma_extra_debug_sigprint();
typedef typename get_pod_type<eT>::result T;
if(A.is_empty()) { return uword(0); }
const uword N = (std::min)(A.n_rows, A.n_cols);
podarray<T> diag_abs_vals(N);
T max_abs_Aii = T(0);
for(uword i=0; i<N; ++i)
{
const eT Aii = A.at(i,i);
const T abs_Aii = std::abs(Aii);
if(arma_isnan(Aii)) { arma_stop_runtime_error("rank(): detected NaN"); return uword(0); }
diag_abs_vals[i] = abs_Aii;
max_abs_Aii = (abs_Aii > 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<T>::epsilon(); }
uword count = 0;
for(uword i=0; i<N; ++i) { count += (diag_abs_vals[i] > tol) ? uword(1) : uword(0); }
return count;
}
template<typename T1>
arma_warn_unused
inline
typename enable_if2< is_supported_blas_type<typename T1::elem_type>::value, uword >::result
rank_sym
(
const Base<typename T1::elem_type,T1>& X,
typename T1::pod_type tol = 0.0
)
rank(const Base<typename T1::elem_type,T1>& 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<eT> 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<T1>::value || A.is_diagmat())
// {
// arma_extra_debug_print("rank(): detected diagonal matrix");
// // TODO: return rank_diagmat(A, tol);
// }
Col<T> 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<T>::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<typename T1>
arma_warn_unused
inline
uword
rank
(
const Base<typename T1::elem_type,T1>& X,
typename T1::pod_type tol = 0.0,
const typename arma_blas_type_only<typename T1::elem_type>::result* junk = nullptr
)
typename enable_if2< is_supported_blas_type<typename T1::elem_type>::value, bool >::result
rank(uword& out, const Base<typename T1::elem_type,T1>& 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<eT> A(X.get_ref());
if(A.is_empty()) { return uword(0); }
if(is_op_diagmat<T1>::value || A.is_diagmat())
{
arma_extra_debug_print("rank(): detected diagonal matrix");
return rank_diagmat(A, tol);
}
if(is_cx<eT>::no && A.is_symmetric())
{
arma_extra_debug_print("rank(): detected symmetric matrix");
return rank_sym(A, tol);
}
if(is_cx<eT>::yes && A.is_hermitian())
{
arma_extra_debug_print("rank(): detected hermitian matrix");
return rank_sym(A, tol);
}
Col<T> 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<T>::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);
}
include/armadillo_bits/op_rank_bones.hpp
+41
View File
@@ -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<typename T1> inline static bool apply(uword& out, const Base<typename T1::elem_type,T1>& expr, const typename T1::pod_type tol);
template<typename eT> inline static bool apply_gen(uword& out, Mat<eT>& A, typename get_pod_type<eT>::result tol);
template<typename eT> inline static bool apply_sym(uword& out, Mat<eT>& A, typename get_pod_type<eT>::result tol);
template<typename eT> inline static bool apply_diag(uword& out, Mat<eT>& A, typename get_pod_type<eT>::result tol);
};
//! @}
include/armadillo_bits/op_rank_meat.hpp
+175
View File
@@ -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<typename T1>
inline
bool
op_rank::apply(uword& out, const Base<typename T1::elem_type,T1>& expr, typename T1::pod_type tol)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT;
Mat<eT> A(expr.get_ref());
if(A.is_empty()) { out = uword(0); return true; }
if(is_op_diagmat<T1>::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<eT>::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<eT>::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<typename eT>
inline
bool
op_rank::apply_diag(uword& out, Mat<eT>& A, typename get_pod_type<eT>::result tol)
{
arma_extra_debug_sigprint();
typedef typename get_pod_type<eT>::result T;
const uword N = (std::min)(A.n_rows, A.n_cols);
podarray<T> diag_abs_vals(N);
T max_abs_Aii = T(0);
for(uword i=0; i<N; ++i)
{
const eT Aii = A.at(i,i);
const T abs_Aii = std::abs(Aii);
if(arma_isnan(Aii)) { out = uword(0); return false; }
diag_abs_vals[i] = abs_Aii;
max_abs_Aii = (abs_Aii > 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<T>::epsilon(); }
uword count = 0;
for(uword i=0; i<N; ++i) { count += (diag_abs_vals[i] > tol) ? uword(1) : uword(0); }
out = count;
return true;
}
template<typename eT>
inline
bool
op_rank::apply_sym(uword& out, Mat<eT>& A, typename get_pod_type<eT>::result tol)
{
arma_extra_debug_sigprint();
typedef typename get_pod_type<eT>::result T;
if(A.is_square() == false) { out = uword(0); return false; }
Col<T> 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<T>::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<typename eT>
inline
bool
op_rank::apply_gen(uword& out, Mat<eT>& A, typename get_pod_type<eT>::result tol)
{
arma_extra_debug_sigprint();
typedef typename get_pod_type<eT>::result T;
Col<T> 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<T>::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;
}
//! @}