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

move sparse matrix norms to separate files

Browse Source
This commit is contained in:
conrad
2021-03-25 15:22:49 +10:00
parent 14d0b35966
commit be701c8e1f
6 changed files with 152 additions and 109 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/armadillo
+2
View File
@@ -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"
include/armadillo_bits/fn_norm.hpp
+3 -3
View File
@@ -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') )
include/armadillo_bits/op_norm_bones.hpp
-12
View File
@@ -23,8 +23,6 @@ class op_norm
{
public:
// norms for dense vectors and matrices
template<typename T1> arma_hot inline static typename T1::pod_type vec_norm_1(const Proxy<T1>& P, const typename arma_not_cx<typename T1::elem_type>::result* junk = nullptr);
template<typename T1> arma_hot inline static typename T1::pod_type vec_norm_1(const Proxy<T1>& P, const typename arma_cx_only<typename T1::elem_type>::result* junk = nullptr);
template<typename eT> arma_hot inline static eT vec_norm_1_direct_std(const Mat<eT>& X);
@@ -45,16 +43,6 @@ class op_norm
template<typename eT> inline static typename get_pod_type<eT>::result mat_norm_2(const Mat<eT>& X);
template<typename eT> inline static typename get_pod_type<eT>::result mat_norm_inf(const Mat<eT>& X);
// norms for sparse matrices
template<typename T1> inline static typename T1::pod_type mat_norm_1(const SpProxy<T1>& P);
template<typename T1> inline static typename T1::pod_type mat_norm_2(const SpProxy<T1>& P, const typename arma_real_only<typename T1::elem_type>::result* junk = nullptr);
template<typename T1> inline static typename T1::pod_type mat_norm_2(const SpProxy<T1>& P, const typename arma_cx_only<typename T1::elem_type>::result* junk = nullptr);
template<typename T1> inline static typename T1::pod_type mat_norm_inf(const SpProxy<T1>& P);
};
include/armadillo_bits/op_norm_meat.hpp
-94
View File
@@ -912,98 +912,4 @@ op_norm::mat_norm_inf(const Mat<eT>& X)
//
// norms for sparse matrices
template<typename T1>
inline
typename T1::pod_type
op_norm::mat_norm_1(const SpProxy<T1>& 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<typename T1>
inline
typename T1::pod_type
op_norm::mat_norm_2(const SpProxy<T1>& P, const typename arma_real_only<typename T1::elem_type>::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<typename SpProxy<T1>::stored_type> tmp(P.Q);
const SpMat<eT>& A = tmp.M;
const SpMat<eT> B = trans(A);
const SpMat<eT> C = (A.n_rows <= A.n_cols) ? (A*B) : (B*A);
Col<T> eigval;
eigs_sym(eigval, C, 1);
return (eigval.n_elem > 0) ? std::sqrt(eigval[0]) : T(0);
}
template<typename T1>
inline
typename T1::pod_type
op_norm::mat_norm_2(const SpProxy<T1>& P, const typename arma_cx_only<typename T1::elem_type>::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<typename SpProxy<T1>::stored_type> tmp(P.Q);
const SpMat<eT>& A = tmp.M;
const SpMat<eT> B = trans(A);
const SpMat<eT> C = (A.n_rows <= A.n_cols) ? (A*B) : (B*A);
Col<eT> eigval;
eigs_gen(eigval, C, 1);
return (eigval.n_elem > 0) ? std::sqrt(std::real(eigval[0])) : T(0);
}
template<typename T1>
inline
typename T1::pod_type
op_norm::mat_norm_inf(const SpProxy<T1>& 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) );
}
//! @}
include/armadillo_bits/spop_norm_bones.hpp
+36
View File
@@ -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<typename T1> inline static typename T1::pod_type mat_norm_1(const SpProxy<T1>& P);
template<typename T1> inline static typename T1::pod_type mat_norm_2(const SpProxy<T1>& P, const typename arma_real_only<typename T1::elem_type>::result* junk = nullptr);
template<typename T1> inline static typename T1::pod_type mat_norm_2(const SpProxy<T1>& P, const typename arma_cx_only<typename T1::elem_type>::result* junk = nullptr);
template<typename T1> inline static typename T1::pod_type mat_norm_inf(const SpProxy<T1>& P);
};
//! @}
include/armadillo_bits/spop_norm_meat.hpp
+111
View File
@@ -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<typename T1>
inline
typename T1::pod_type
spop_norm::mat_norm_1(const SpProxy<T1>& 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<typename T1>
inline
typename T1::pod_type
spop_norm::mat_norm_2(const SpProxy<T1>& P, const typename arma_real_only<typename T1::elem_type>::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<typename SpProxy<T1>::stored_type> tmp(P.Q);
const SpMat<eT>& A = tmp.M;
const SpMat<eT> B = trans(A);
const SpMat<eT> C = (A.n_rows <= A.n_cols) ? (A*B) : (B*A);
Col<T> eigval;
eigs_sym(eigval, C, 1);
return (eigval.n_elem > 0) ? std::sqrt(eigval[0]) : T(0);
}
template<typename T1>
inline
typename T1::pod_type
spop_norm::mat_norm_2(const SpProxy<T1>& P, const typename arma_cx_only<typename T1::elem_type>::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<typename SpProxy<T1>::stored_type> tmp(P.Q);
const SpMat<eT>& A = tmp.M;
const SpMat<eT> B = trans(A);
const SpMat<eT> C = (A.n_rows <= A.n_cols) ? (A*B) : (B*A);
Col<eT> eigval;
eigs_gen(eigval, C, 1);
return (eigval.n_elem > 0) ? std::sqrt(std::real(eigval[0])) : T(0);
}
template<typename T1>
inline
typename T1::pod_type
spop_norm::mat_norm_inf(const SpProxy<T1>& 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) );
}
//! @}