fn_trace.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2010 NICTA and the authors listed below
00002 // http://nicta.com.au
00003 // 
00004 // Authors:
00005 // - Conrad Sanderson (conradsand at ieee dot org)
00006 // 
00007 // This file is part of the Armadillo C++ library.
00008 // It is provided without any warranty of fitness
00009 // for any purpose. You can redistribute this file
00010 // and/or modify it under the terms of the GNU
00011 // Lesser General Public License (LGPL) as published
00012 // by the Free Software Foundation, either version 3
00013 // of the License or (at your option) any later version.
00014 // (see http://www.opensource.org/licenses for more info)
00015 
00016 
00017 //! \addtogroup fn_trace
00018 //! @{
00019 
00020 
00021 //! Immediate trace (sum of diagonal elements) of a square dense matrix
00022 template<typename T1>
00023 inline
00024 arma_warn_unused
00025 typename T1::elem_type
00026 trace(const Base<typename T1::elem_type,T1>& X)
00027   {
00028   arma_extra_debug_sigprint();
00029   
00030   typedef typename T1::elem_type eT;
00031   
00032   const Proxy<T1> A(X.get_ref());
00033 
00034   arma_debug_check( (A.n_rows != A.n_cols), "trace(): matrix must be square sized" );
00035   
00036   const u32 N   = A.n_rows;
00037         eT  val = eT(0);
00038   
00039   for(u32 i=0; i<N; ++i)
00040     {
00041     val += A.at(i,i);
00042     }
00043   
00044   return val;
00045   }
00046 
00047 
00048 
00049 template<typename T1>
00050 inline
00051 arma_warn_unused
00052 typename T1::elem_type
00053 trace(const Op<T1, op_diagmat>& X)
00054   {
00055   arma_extra_debug_sigprint();
00056   
00057   typedef typename T1::elem_type eT;
00058   
00059   const diagmat_proxy<T1> A(X.m);
00060   
00061   const u32 N = A.n_elem;
00062   
00063   eT val = eT(0);
00064   
00065   for(u32 i=0; i<N; ++i)
00066     {
00067     val += A[i];
00068     }
00069   
00070   return val;
00071   }
00072 
00073 
00074 //! speedup for trace(A*B), where the result of A*B is a square sized matrix
00075 template<typename T1, typename T2>
00076 inline
00077 arma_warn_unused
00078 typename T1::elem_type
00079 trace(const Glue<T1, T2, glue_times>& X)
00080   {
00081   arma_extra_debug_sigprint();
00082   
00083   typedef typename T1::elem_type eT;
00084   
00085   const unwrap<T1> tmp1(X.A);
00086   const unwrap<T2> tmp2(X.B);
00087   
00088   const Mat<eT>& A = tmp1.M;
00089   const Mat<eT>& B = tmp2.M;
00090   
00091   arma_debug_assert_mul_size(A, B, "matrix multiply");
00092   
00093   arma_debug_check( (A.n_rows != B.n_cols), "trace(): matrix must be square sized" );
00094   
00095   const u32 N1  = A.n_rows;
00096   const u32 N2  = A.n_cols;
00097         eT  val = eT(0);
00098   
00099   for(u32 i=0; i<N1; ++i)
00100     {
00101     const eT* B_colmem = B.colptr(i);
00102           eT  acc      = eT(0);
00103     
00104     for(u32 j=0; j<N2; ++j)
00105       {
00106       acc += A.at(i,j) * B_colmem[j];
00107       }
00108     
00109     val += acc;
00110     }
00111   
00112   return val;
00113   }
00114 
00115 
00116 
00117 //! @}