max@0: // Copyright (C) 2008-2011 NICTA (www.nicta.com.au)
max@0: // Copyright (C) 2008-2011 Conrad Sanderson
max@0: // 
max@0: // This file is part of the Armadillo C++ library.
max@0: // It is provided without any warranty of fitness
max@0: // for any purpose. You can redistribute this file
max@0: // and/or modify it under the terms of the GNU
max@0: // Lesser General Public License (LGPL) as published
max@0: // by the Free Software Foundation, either version 3
max@0: // of the License or (at your option) any later version.
max@0: // (see http://www.opensource.org/licenses for more info)
max@0: 
max@0: 
max@0: //! \addtogroup fn_norm
max@0: //! @{
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: arma_hot
max@0: inline
max@0: typename T1::pod_type
max@0: arma_vec_norm_1(const Proxy<T1>& A)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::pod_type T;
max@0: 
max@0:   T acc = T(0);
max@0:     
max@0:   if(Proxy<T1>::prefer_at_accessor == false)
max@0:     {
max@0:     typename Proxy<T1>::ea_type P = A.get_ea();
max@0:     
max@0:     const uword N = A.get_n_elem();
max@0:     
max@0:     uword i,j;
max@0:     
max@0:     for(i=0, j=1; j<N; i+=2, j+=2)
max@0:       {
max@0:       acc += std::abs(P[i]);
max@0:       acc += std::abs(P[j]);
max@0:       }
max@0:     
max@0:     if(i < N)
max@0:       {
max@0:       acc += std::abs(P[i]);
max@0:       }
max@0:     }
max@0:   else
max@0:     {
max@0:     const uword n_rows = A.get_n_rows();
max@0:     const uword n_cols = A.get_n_cols();
max@0:     
max@0:     for(uword col=0; col<n_cols; ++col)
max@0:       {
max@0:       uword i,j;
max@0:       
max@0:       for(i=0, j=1; j<n_rows; i+=2, j+=2)
max@0:         {
max@0:         acc += std::abs(A.at(i,col));
max@0:         acc += std::abs(A.at(j,col));
max@0:         }
max@0:       
max@0:       if(i < n_rows)
max@0:         {
max@0:         acc += std::abs(A.at(i,col));
max@0:         }
max@0:       }
max@0:     }
max@0:     
max@0:   return acc;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: arma_hot
max@0: inline
max@0: typename T1::pod_type
max@0: arma_vec_norm_2(const Proxy<T1>& A, const typename arma_not_cx<typename T1::elem_type>::result* junk = 0)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   typedef typename T1::pod_type T;
max@0:   
max@0:   T acc = T(0);
max@0:   
max@0:   if(Proxy<T1>::prefer_at_accessor == false)
max@0:     {
max@0:     typename Proxy<T1>::ea_type P = A.get_ea();
max@0:     
max@0:     const uword N = A.get_n_elem();
max@0:     
max@0:     uword i,j;
max@0:     
max@0:     for(i=0, j=1; j<N; i+=2, j+=2)
max@0:       {
max@0:       const T tmp_i = P[i];
max@0:       const T tmp_j = P[j];
max@0:       
max@0:       acc += tmp_i * tmp_i;
max@0:       acc += tmp_j * tmp_j;
max@0:       }
max@0:     
max@0:     if(i < N)
max@0:       {
max@0:       const T tmp_i = P[i];
max@0:       
max@0:       acc += tmp_i * tmp_i;
max@0:       }
max@0:     }
max@0:   else
max@0:     {
max@0:     const uword n_rows = A.get_n_rows();
max@0:     const uword n_cols = A.get_n_cols();
max@0:     
max@0:     for(uword col=0; col<n_cols; ++col)
max@0:       {
max@0:       uword i,j;
max@0:       
max@0:       for(i=0, j=1; j<n_rows; i+=2, j+=2)
max@0:         {
max@0:         const T tmp_i = A.at(i,col);
max@0:         const T tmp_j = A.at(j,col);
max@0:         
max@0:         acc += tmp_i * tmp_i;
max@0:         acc += tmp_j * tmp_j;
max@0:         }
max@0:       
max@0:       if(i < n_rows)
max@0:         {
max@0:         const T tmp_i = A.at(i,col);
max@0:         
max@0:         acc += tmp_i * tmp_i;
max@0:         }
max@0:       }
max@0:     }
max@0:   
max@0:   return std::sqrt(acc);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: arma_hot
max@0: inline
max@0: typename T1::pod_type
max@0: arma_vec_norm_2(const Proxy<T1>& A, const typename arma_cx_only<typename T1::elem_type>::result* junk = 0)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   typedef typename T1::pod_type T;
max@0:   
max@0:   T acc = T(0);
max@0:   
max@0:   if(Proxy<T1>::prefer_at_accessor == false)
max@0:     {
max@0:     typename Proxy<T1>::ea_type P = A.get_ea();
max@0:     
max@0:     const uword N = A.get_n_elem();
max@0:     
max@0:     for(uword i=0; i<N; ++i)
max@0:       {
max@0:       const T tmp = std::abs(P[i]);
max@0:       acc += tmp*tmp;
max@0:       }
max@0:     }
max@0:   else
max@0:     {
max@0:     const uword n_rows = A.get_n_rows();
max@0:     const uword n_cols = A.get_n_cols();
max@0:     
max@0:     for(uword col=0; col<n_cols; ++col)
max@0:     for(uword row=0; row<n_rows; ++row)
max@0:       {
max@0:       const T tmp = std::abs(A.at(row,col));
max@0:       acc += tmp*tmp;
max@0:       }
max@0:     }
max@0:   
max@0:   return std::sqrt(acc);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: arma_hot
max@0: inline
max@0: typename T1::pod_type
max@0: arma_vec_norm_k(const Proxy<T1>& A, const int k)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::pod_type T;
max@0:   
max@0:   T acc = T(0);
max@0:   
max@0:   if(Proxy<T1>::prefer_at_accessor == false)
max@0:     {
max@0:     typename Proxy<T1>::ea_type P = A.get_ea();
max@0:     
max@0:     const uword N = A.get_n_elem();
max@0:     
max@0:     uword i,j;
max@0:     
max@0:     for(i=0, j=1; j<N; i+=2, j+=2)
max@0:       {
max@0:       acc += std::pow(std::abs(P[i]), k);
max@0:       acc += std::pow(std::abs(P[j]), k);
max@0:       }
max@0:     
max@0:     if(i < N)
max@0:       {
max@0:       acc += std::pow(std::abs(P[i]), k);
max@0:       }
max@0:     }
max@0:   else
max@0:     {
max@0:     const uword n_rows = A.get_n_rows();
max@0:     const uword n_cols = A.get_n_cols();
max@0:     
max@0:     for(uword col=0; col<n_cols; ++col)
max@0:     for(uword row=0; row<n_rows; ++row)
max@0:       {
max@0:       acc += std::pow(std::abs(A.at(row,col)), k);
max@0:       }
max@0:     }
max@0:   
max@0:   return std::pow(acc, T(1)/T(k));
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: arma_hot
max@0: inline
max@0: typename T1::pod_type
max@0: arma_vec_norm_max(const Proxy<T1>& A)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::pod_type       T;
max@0:   typedef typename Proxy<T1>::ea_type ea_type;
max@0:   
max@0:         ea_type P = A.get_ea();
max@0:   const uword     N = A.get_n_elem();
max@0:   
max@0:   T max_val = (N != 1) ? priv::most_neg<T>() : std::abs(P[0]);
max@0:   
max@0:   uword i,j;
max@0:   
max@0:   for(i=0, j=1; j<N; i+=2, j+=2)
max@0:     {
max@0:     const T tmp_i = std::abs(P[i]);
max@0:     const T tmp_j = std::abs(P[j]);
max@0:     
max@0:     if(max_val < tmp_i) { max_val = tmp_i; }
max@0:     if(max_val < tmp_j) { max_val = tmp_j; }
max@0:     }
max@0:   
max@0:   if(i < N)
max@0:     {
max@0:     const T tmp_i = std::abs(P[i]);
max@0:     
max@0:     if(max_val < tmp_i) { max_val = tmp_i; }
max@0:     }
max@0:   
max@0:   return max_val;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: arma_hot
max@0: inline
max@0: typename T1::pod_type
max@0: arma_vec_norm_min(const Proxy<T1>& A)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::pod_type       T;
max@0:   typedef typename Proxy<T1>::ea_type ea_type;
max@0:   
max@0:         ea_type P = A.get_ea();
max@0:   const uword     N = A.get_n_elem();
max@0:   
max@0:   T min_val = (N != 1) ? priv::most_pos<T>() : std::abs(P[0]);
max@0:   
max@0:   uword i,j;
max@0:   
max@0:   for(i=0, j=1; j<N; i+=2, j+=2)
max@0:     {
max@0:     const T tmp_i = std::abs(P[i]);
max@0:     const T tmp_j = std::abs(P[j]);
max@0:     
max@0:     if(min_val > tmp_i) { min_val = tmp_i; }
max@0:     if(min_val > tmp_j) { min_val = tmp_j; }
max@0:     }
max@0:   
max@0:   if(i < N)
max@0:     {
max@0:     const T tmp_i = std::abs(P[i]);
max@0:     
max@0:     if(min_val > tmp_i) { min_val = tmp_i; }
max@0:     }
max@0:   
max@0:   return min_val;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: inline
max@0: typename T1::pod_type
max@0: arma_mat_norm_1(const Proxy<T1>& A)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::elem_type eT;
max@0:   typedef typename T1::pod_type   T;
max@0:   
max@0:   const unwrap<typename Proxy<T1>::stored_type> tmp(A.Q);
max@0:   const Mat<eT>& X = tmp.M;
max@0:   
max@0:   // TODO: this can be sped up with a dedicated implementation
max@0:   return as_scalar( max( sum(abs(X)), 1) );
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: inline
max@0: typename T1::pod_type
max@0: arma_mat_norm_2(const Proxy<T1>& A)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::elem_type eT;
max@0:   typedef typename T1::pod_type   T;
max@0:   
max@0:   const unwrap<typename Proxy<T1>::stored_type> tmp(A.Q);
max@0:   const Mat<eT>& X = tmp.M;
max@0:   
max@0:   Col<T> S;
max@0:   svd(S, X);
max@0:   
max@0:   return (S.n_elem > 0) ? max(S) : T(0);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: inline
max@0: typename T1::pod_type
max@0: arma_mat_norm_inf(const Proxy<T1>& A)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::elem_type eT;
max@0:   typedef typename T1::pod_type   T;
max@0:   
max@0:   const unwrap<typename Proxy<T1>::stored_type> tmp(A.Q);
max@0:   const Mat<eT>& X = tmp.M;
max@0:   
max@0:   // TODO: this can be sped up with a dedicated implementation
max@0:   return as_scalar( max( sum(abs(X),1) ) );
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: inline
max@0: arma_warn_unused
max@0: typename T1::pod_type
max@0: norm
max@0:   (
max@0:   const Base<typename T1::elem_type,T1>& X,
max@0:   const uword k,
max@0:   const typename arma_float_or_cx_only<typename T1::elem_type>::result* junk = 0
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   typedef typename T1::elem_type eT;
max@0:   typedef typename T1::pod_type   T;
max@0:   
max@0:   const Proxy<T1> A(X.get_ref());
max@0:   
max@0:   if(A.get_n_elem() == 0)
max@0:     {
max@0:     return T(0);
max@0:     }
max@0:   
max@0:   const bool is_vec = (A.get_n_rows() == 1) || (A.get_n_cols() == 1);
max@0:   
max@0:   if(is_vec == true)
max@0:     {
max@0:     switch(k)
max@0:       {
max@0:       case 1:
max@0:         return arma_vec_norm_1(A);
max@0:         break;
max@0:       
max@0:       case 2:
max@0:         return arma_vec_norm_2(A);
max@0:         break;
max@0:       
max@0:       default:
max@0:         {
max@0:         arma_debug_check( (k == 0), "norm(): k must be greater than zero"   );
max@0:         return arma_vec_norm_k(A, int(k));
max@0:         }
max@0:       }
max@0:     }
max@0:   else
max@0:     {
max@0:     switch(k)
max@0:       {
max@0:       case 1:
max@0:         return arma_mat_norm_1(A);
max@0:         break;
max@0:       
max@0:       case 2:
max@0:         return arma_mat_norm_2(A);
max@0:         break;
max@0:       
max@0:       default:
max@0:         arma_stop("norm(): unsupported matrix norm type");
max@0:         return T(0);
max@0:       }
max@0:     }
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: inline
max@0: arma_warn_unused
max@0: typename T1::pod_type
max@0: norm
max@0:   (
max@0:   const Base<typename T1::elem_type,T1>& X,
max@0:   const char* method,
max@0:   const typename arma_float_or_cx_only<typename T1::elem_type>::result* junk = 0
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   typedef typename T1::elem_type eT;
max@0:   typedef typename T1::pod_type   T;
max@0:   
max@0:   const Proxy<T1> A(X.get_ref());
max@0:   
max@0:   if(A.get_n_elem() == 0)
max@0:     {
max@0:     return T(0);
max@0:     }
max@0:   
max@0:   const char sig    = method[0];
max@0:   const bool is_vec = (A.get_n_rows() == 1) || (A.get_n_cols() == 1);
max@0:   
max@0:   if(is_vec == true)
max@0:     {
max@0:     if( (sig == 'i') || (sig == 'I') || (sig == '+') )   // max norm
max@0:       {
max@0:       return arma_vec_norm_max(A);
max@0:       }
max@0:     else
max@0:     if(sig == '-')   // min norm
max@0:       {
max@0:       return arma_vec_norm_min(A);
max@0:       }
max@0:     else
max@0:     if( (sig == 'f') || (sig == 'F') )
max@0:       {
max@0:       return arma_vec_norm_2(A);
max@0:       }
max@0:     else
max@0:       {
max@0:       arma_stop("norm(): unsupported vector norm type");
max@0:       return T(0);
max@0:       }
max@0:     }
max@0:   else
max@0:     {
max@0:     if( (sig == 'i') || (sig == 'I') || (sig == '+') )   // inf norm
max@0:       {
max@0:       return arma_mat_norm_inf(A);
max@0:       }
max@0:     else
max@0:     if( (sig == 'f') || (sig == 'F') )
max@0:       {
max@0:       return arma_vec_norm_2(A);
max@0:       }
max@0:     else
max@0:       {
max@0:       arma_stop("norm(): unsupported matrix norm type");
max@0:       return T(0);
max@0:       }
max@0:     }
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! @}