From be701c8e1f9672040f2c8da8aaa6b61aafc96ba6 Mon Sep 17 00:00:00 2001 From: conrad Date: Thu, 25 Mar 2021 15:22:49 +1000 Subject: [PATCH] move sparse matrix norms to separate files --- include/armadillo | 2 + include/armadillo_bits/fn_norm.hpp | 6 +- include/armadillo_bits/op_norm_bones.hpp | 12 --- include/armadillo_bits/op_norm_meat.hpp | 94 ----------------- include/armadillo_bits/spop_norm_bones.hpp | 36 +++++++ include/armadillo_bits/spop_norm_meat.hpp | 111 +++++++++++++++++++++ 6 files changed, 152 insertions(+), 109 deletions(-) create mode 100644 include/armadillo_bits/spop_norm_bones.hpp create mode 100644 include/armadillo_bits/spop_norm_meat.hpp diff --git a/include/armadillo b/include/armadillo index b1e25ce9..97840b5f 100644 --- a/include/armadillo +++ b/include/armadillo @@ -319,6 +319,7 @@ namespace arma #include "armadillo_bits/spop_reverse_bones.hpp" #include "armadillo_bits/spop_repmat_bones.hpp" #include "armadillo_bits/spop_vectorise_bones.hpp" + #include "armadillo_bits/spop_norm_bones.hpp" #include "armadillo_bits/spglue_plus_bones.hpp" #include "armadillo_bits/spglue_minus_bones.hpp" @@ -746,6 +747,7 @@ namespace arma #include "armadillo_bits/spop_reverse_meat.hpp" #include "armadillo_bits/spop_repmat_meat.hpp" #include "armadillo_bits/spop_vectorise_meat.hpp" + #include "armadillo_bits/spop_norm_meat.hpp" #include "armadillo_bits/spglue_plus_meat.hpp" #include "armadillo_bits/spglue_minus_meat.hpp" diff --git a/include/armadillo_bits/fn_norm.hpp b/include/armadillo_bits/fn_norm.hpp index c4b3d94b..20e84b53 100644 --- a/include/armadillo_bits/fn_norm.hpp +++ b/include/armadillo_bits/fn_norm.hpp @@ -164,8 +164,8 @@ norm } else { - if(k == uword(1)) { return op_norm::mat_norm_1(P); } - if(k == uword(2)) { return op_norm::mat_norm_2(P); } + if(k == uword(1)) { return spop_norm::mat_norm_1(P); } + if(k == uword(2)) { return spop_norm::mat_norm_2(P); } arma_stop_logic_error("norm(): unsupported or unimplemented norm type for sparse matrices"); } @@ -234,7 +234,7 @@ norm { if( (sig == 'i') || (sig == 'I') || (sig == '+') ) // inf norm { - return op_norm::mat_norm_inf(P); + return spop_norm::mat_norm_inf(P); } else if( (sig == 'f') || (sig == 'F') ) diff --git a/include/armadillo_bits/op_norm_bones.hpp b/include/armadillo_bits/op_norm_bones.hpp index a131bdd4..cc778e5e 100644 --- a/include/armadillo_bits/op_norm_bones.hpp +++ b/include/armadillo_bits/op_norm_bones.hpp @@ -23,8 +23,6 @@ class op_norm { public: - // norms for dense vectors and matrices - template arma_hot inline static typename T1::pod_type vec_norm_1(const Proxy& P, const typename arma_not_cx::result* junk = nullptr); template arma_hot inline static typename T1::pod_type vec_norm_1(const Proxy& P, const typename arma_cx_only::result* junk = nullptr); template arma_hot inline static eT vec_norm_1_direct_std(const Mat& X); @@ -45,16 +43,6 @@ class op_norm template inline static typename get_pod_type::result mat_norm_2(const Mat& X); template inline static typename get_pod_type::result mat_norm_inf(const Mat& X); - - - // norms for sparse matrices - - template inline static typename T1::pod_type mat_norm_1(const SpProxy& P); - - template inline static typename T1::pod_type mat_norm_2(const SpProxy& P, const typename arma_real_only::result* junk = nullptr); - template inline static typename T1::pod_type mat_norm_2(const SpProxy& P, const typename arma_cx_only::result* junk = nullptr); - - template inline static typename T1::pod_type mat_norm_inf(const SpProxy& P); }; diff --git a/include/armadillo_bits/op_norm_meat.hpp b/include/armadillo_bits/op_norm_meat.hpp index 2a15c782..95654669 100644 --- a/include/armadillo_bits/op_norm_meat.hpp +++ b/include/armadillo_bits/op_norm_meat.hpp @@ -912,98 +912,4 @@ op_norm::mat_norm_inf(const Mat& X) -// -// norms for sparse matrices - - - -template -inline -typename T1::pod_type -op_norm::mat_norm_1(const SpProxy& P) - { - arma_extra_debug_sigprint(); - - // TODO: this can be sped up with a dedicated implementation - return as_scalar( max( sum(abs(P.Q), 0), 1) ); - } - - - -template -inline -typename T1::pod_type -op_norm::mat_norm_2(const SpProxy& P, const typename arma_real_only::result* junk) - { - arma_extra_debug_sigprint(); - arma_ignore(junk); - - // norm = sqrt( largest eigenvalue of (A^H)*A ), where ^H is the conjugate transpose - // http://math.stackexchange.com/questions/4368/computing-the-largest-eigenvalue-of-a-very-large-sparse-matrix - - typedef typename T1::elem_type eT; - typedef typename T1::pod_type T; - - const unwrap_spmat::stored_type> tmp(P.Q); - - const SpMat& A = tmp.M; - const SpMat B = trans(A); - - const SpMat C = (A.n_rows <= A.n_cols) ? (A*B) : (B*A); - - Col eigval; - eigs_sym(eigval, C, 1); - - return (eigval.n_elem > 0) ? std::sqrt(eigval[0]) : T(0); - } - - - -template -inline -typename T1::pod_type -op_norm::mat_norm_2(const SpProxy& P, const typename arma_cx_only::result* junk) - { - arma_extra_debug_sigprint(); - arma_ignore(junk); - - typedef typename T1::elem_type eT; - typedef typename T1::pod_type T; - - // we're calling eigs_gen(), which currently requires ARPACK - #if !defined(ARMA_USE_ARPACK) - { - arma_stop_logic_error("norm(): use of ARPACK must be enabled for norm of complex matrices"); - return T(0); - } - #endif - - const unwrap_spmat::stored_type> tmp(P.Q); - - const SpMat& A = tmp.M; - const SpMat B = trans(A); - - const SpMat C = (A.n_rows <= A.n_cols) ? (A*B) : (B*A); - - Col eigval; - eigs_gen(eigval, C, 1); - - return (eigval.n_elem > 0) ? std::sqrt(std::real(eigval[0])) : T(0); - } - - - -template -inline -typename T1::pod_type -op_norm::mat_norm_inf(const SpProxy& P) - { - arma_extra_debug_sigprint(); - - // TODO: this can be sped up with a dedicated implementation - return as_scalar( max( sum(abs(P.Q), 1), 0) ); - } - - - //! @} diff --git a/include/armadillo_bits/spop_norm_bones.hpp b/include/armadillo_bits/spop_norm_bones.hpp new file mode 100644 index 00000000..84d0bf18 --- /dev/null +++ b/include/armadillo_bits/spop_norm_bones.hpp @@ -0,0 +1,36 @@ +// 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_norm +//! @{ + + +class spop_norm + : public traits_op_default + { + public: + + template inline static typename T1::pod_type mat_norm_1(const SpProxy& P); + + template inline static typename T1::pod_type mat_norm_2(const SpProxy& P, const typename arma_real_only::result* junk = nullptr); + template inline static typename T1::pod_type mat_norm_2(const SpProxy& P, const typename arma_cx_only::result* junk = nullptr); + + template inline static typename T1::pod_type mat_norm_inf(const SpProxy& P); + }; + + + +//! @} diff --git a/include/armadillo_bits/spop_norm_meat.hpp b/include/armadillo_bits/spop_norm_meat.hpp new file mode 100644 index 00000000..cce8fe95 --- /dev/null +++ b/include/armadillo_bits/spop_norm_meat.hpp @@ -0,0 +1,111 @@ +// 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_norm +//! @{ + + + +template +inline +typename T1::pod_type +spop_norm::mat_norm_1(const SpProxy& P) + { + arma_extra_debug_sigprint(); + + // TODO: this can be sped up with a dedicated implementation + return as_scalar( max( sum(abs(P.Q), 0), 1) ); + } + + + +template +inline +typename T1::pod_type +spop_norm::mat_norm_2(const SpProxy& P, const typename arma_real_only::result* junk) + { + arma_extra_debug_sigprint(); + arma_ignore(junk); + + // norm = sqrt( largest eigenvalue of (A^H)*A ), where ^H is the conjugate transpose + // http://math.stackexchange.com/questions/4368/computing-the-largest-eigenvalue-of-a-very-large-sparse-matrix + + typedef typename T1::elem_type eT; + typedef typename T1::pod_type T; + + const unwrap_spmat::stored_type> tmp(P.Q); + + const SpMat& A = tmp.M; + const SpMat B = trans(A); + + const SpMat C = (A.n_rows <= A.n_cols) ? (A*B) : (B*A); + + Col eigval; + eigs_sym(eigval, C, 1); + + return (eigval.n_elem > 0) ? std::sqrt(eigval[0]) : T(0); + } + + + +template +inline +typename T1::pod_type +spop_norm::mat_norm_2(const SpProxy& P, const typename arma_cx_only::result* junk) + { + arma_extra_debug_sigprint(); + arma_ignore(junk); + + typedef typename T1::elem_type eT; + typedef typename T1::pod_type T; + + // we're calling eigs_gen(), which currently requires ARPACK + #if !defined(ARMA_USE_ARPACK) + { + arma_stop_logic_error("norm(): use of ARPACK must be enabled for norm of complex matrices"); + return T(0); + } + #endif + + const unwrap_spmat::stored_type> tmp(P.Q); + + const SpMat& A = tmp.M; + const SpMat B = trans(A); + + const SpMat C = (A.n_rows <= A.n_cols) ? (A*B) : (B*A); + + Col eigval; + eigs_gen(eigval, C, 1); + + return (eigval.n_elem > 0) ? std::sqrt(std::real(eigval[0])) : T(0); + } + + + +template +inline +typename T1::pod_type +spop_norm::mat_norm_inf(const SpProxy& P) + { + arma_extra_debug_sigprint(); + + // TODO: this can be sped up with a dedicated implementation + return as_scalar( max( sum(abs(P.Q), 1), 0) ); + } + + + +//! @}