fn_accu.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2009 NICTA
00002 // 
00003 // Authors:
00004 // - Conrad Sanderson (conradsand at ieee dot org)
00005 // 
00006 // This file is part of the Armadillo C++ library.
00007 // It is provided without any warranty of fitness
00008 // for any purpose. You can redistribute this file
00009 // and/or modify it under the terms of the GNU
00010 // Lesser General Public License (LGPL) as published
00011 // by the Free Software Foundation, either version 3
00012 // of the License or (at your option) any later version.
00013 // (see http://www.opensource.org/licenses for more info)
00014 
00015 
00016 //! \addtogroup fn_accu
00017 //! @{
00018 
00019 
00020 
00021 //! accumulate the elements of a matrix
00022 template<typename T1>
00023 inline
00024 typename T1::elem_type
00025 accu(const Base<typename T1::elem_type,T1>& X)
00026   {
00027   arma_extra_debug_sigprint();
00028   
00029   typedef typename T1::elem_type eT;
00030   
00031   const unwrap<T1> tmp(X.get_ref());
00032   const Mat<eT>& A = tmp.M;
00033   
00034   const u32 A_n_elem = A.n_elem;
00035   const eT* A_mem    = A.mem;
00036   
00037   eT val = eT(0);
00038   
00039   for(u32 i=0; i<A_n_elem; ++i)
00040     {
00041     val += A_mem[i];
00042     }
00043   
00044   return val;
00045   }
00046 
00047 
00048 
00049 //! sum of values along the main diagonal
00050 template<typename T1>
00051 inline
00052 typename T1::elem_type
00053 accu(const Op<T1, op_diagmat>& X)
00054   {
00055   arma_extra_debug_sigprint();
00056   
00057   typedef typename T1::elem_type eT;
00058   
00059   const unwrap<T1> tmp(X.m);
00060   const Mat<eT>& A = tmp.M;
00061   
00062   arma_debug_check( !A.is_square(), "accu(): sum of diagonal values of a non-square matrix requested" );
00063   
00064   eT acc = eT(0);
00065   
00066   for(u32 i=0; i<A.n_rows; ++i)
00067     {
00068     acc += A.at(i,i);
00069     }
00070   
00071   return acc;
00072   }
00073 
00074 
00075 
00076 template<typename eT>
00077 inline
00078 eT
00079 accu(const Op<Mat<eT>, op_diagmat_vec>& X)
00080   {
00081   arma_extra_debug_sigprint();
00082   
00083   const Mat<eT>& A = X.m;
00084   arma_debug_check( !A.is_vec(), "accu(): internal error: expected a vector" );
00085   
00086   return accu(A);
00087   }
00088 
00089 
00090 
00091 //! sum of squares
00092 template<typename T1>
00093 inline
00094 typename T1::elem_type
00095 accu(const Op<T1, op_square>& in)
00096   {
00097   arma_extra_debug_sigprint();
00098   
00099   typedef typename T1::elem_type eT;
00100   
00101   const unwrap<T1> tmp(in.m);
00102   const Mat<eT>& A = tmp.M;
00103   
00104   const u32 A_n_elem = A.n_elem;
00105   const eT* A_mem    = A.mem;
00106   
00107   eT acc = eT(0);
00108   
00109   for(u32 i=0; i<A_n_elem; ++i)
00110     {
00111     const eT val = A_mem[i];
00112     acc += val*val;
00113     }
00114   
00115   return acc;
00116   }
00117 
00118 
00119 
00120 //! sum of square roots
00121 template<typename T1>
00122 inline
00123 typename T1::elem_type
00124 accu(const Op<T1, op_sqrt>& in)
00125   {
00126   arma_extra_debug_sigprint();
00127   
00128   typedef typename T1::elem_type eT;
00129   
00130   const unwrap<T1> tmp(in.m);
00131   const Mat<eT>& A = tmp.M;
00132   
00133   const u32 A_n_elem = A.n_elem;
00134   const eT* A_mem    = A.mem;
00135   
00136   eT acc = eT(0);
00137   for(u32 i=0; i<A_n_elem; ++i)
00138     {
00139     acc += std::sqrt(A_mem[i]);
00140     }
00141   
00142   return acc;
00143   }
00144 
00145 
00146 
00147 //! sum of squares of differences
00148 template<typename T1, typename T2>
00149 inline
00150 typename T1::elem_type
00151 accu(const Op< Glue<T1,T2, glue_minus>, op_square>& in)
00152   {
00153   arma_extra_debug_sigprint();
00154   
00155   typedef typename T1::elem_type eT;
00156   
00157   const unwrap<T1> tmp1(in.m.A);
00158   const unwrap<T2> tmp2(in.m.B);
00159   
00160   const Mat<eT>& A = tmp1.M;
00161   const Mat<eT>& B = tmp2.M;
00162   
00163   arma_debug_assert_same_size(A,B, "accu()");
00164   
00165   const u32 n_elem = A.n_elem;
00166   
00167   eT acc = eT(0);
00168   for(u32 i=0; i<n_elem; ++i)
00169     {
00170     const eT val = A.mem[i] - B.mem[i];
00171     acc += val*val;
00172     }
00173   
00174   return acc;
00175   }
00176 
00177 
00178 
00179 //! accumulate the elements of a subview (submatrix)
00180 template<typename eT>
00181 inline
00182 eT
00183 accu(const subview<eT>& X)
00184   {
00185   arma_extra_debug_sigprint();  
00186   
00187   eT val = eT(0);
00188   for(u32 col=0; col<X.n_cols; ++col)
00189     {
00190     const eT* coldata = X.colptr(col);
00191     
00192     for(u32 row=0; row<X.n_rows; ++row)
00193       {
00194       val += coldata[row];
00195       }
00196     
00197     }
00198   
00199   return val;
00200   }
00201 
00202 
00203 
00204 //! accumulate the elements of a diagview
00205 template<typename eT>
00206 inline
00207 eT
00208 accu(const diagview<eT>& X)
00209   {
00210   arma_extra_debug_sigprint();  
00211   
00212   const u32 n_elem = X.n_elem;
00213   eT val = eT(0);
00214   
00215   for(u32 i=0; i<n_elem; ++i)
00216     {
00217     val += X[i];
00218     }
00219   
00220   return val;
00221   }
00222 
00223 
00224 
00225 //! accumulate the result of A % B, where % is the Schur product (element-wise multiplication)
00226 template<typename eT>
00227 inline
00228 eT
00229 accu_schur(const Mat<eT>& A, const Mat<eT>& B)
00230   {
00231   arma_extra_debug_sigprint();
00232   
00233   arma_debug_assert_same_size(A,B, "accu()");
00234  
00235   const eT* const A_mem = A.mem;
00236   const eT* const B_mem = B.mem;
00237   
00238   const u32 n_elem = A.n_elem;
00239   eT val = eT(0);
00240   
00241   for(u32 i=0; i<n_elem; ++i)
00242     {
00243     val += A_mem[i] * B_mem[i];
00244     }
00245   
00246   return val;
00247   }
00248 
00249 
00250 
00251 //! accumulate the result of A % B, where % is the Schur product (element-wise multiplication)
00252 template<typename eT>
00253 inline
00254 eT
00255 accu(const Glue<Mat<eT>,Mat<eT>,glue_schur>& X)
00256   {
00257   return accu_schur(X.A, X.B);
00258   }
00259 
00260 
00261 
00262 //! accumulate the result of A % B % C, where % is the Schur product (element-wise multiplication)
00263 template<typename eT>
00264 inline
00265 eT
00266 accu(const Glue<Glue<Mat<eT>,Mat<eT>,glue_schur>,Mat<eT>,glue_schur>& X)
00267   {
00268   arma_extra_debug_sigprint();
00269   
00270   const Mat<eT>& A = X.A.A;
00271   const Mat<eT>& B = X.A.B;
00272   const Mat<eT>& C = X.B;
00273   
00274   arma_debug_assert_same_size(A,B, "accu()");
00275   arma_debug_assert_same_size(A,C, "accu()");
00276   
00277   const eT* const A_mem = A.mem;
00278   const eT* const B_mem = B.mem;
00279   const eT* const C_mem = C.mem;
00280   
00281   const u32 n_elem = A.n_elem;
00282   eT val = eT(0);
00283   
00284   for(u32 i=0; i<n_elem; ++i)
00285     {
00286     val += A_mem[i] * B_mem[i] * C_mem[i];
00287     }
00288     
00289   return val;
00290   }
00291 
00292 
00293 
00294 //! \brief
00295 //! accumulate the result of T1 % T2 
00296 //! where % is the Schur product (element-wise multiplication),
00297 //! while T1 and T2 can be 'mat', 'rowvec', 'colvec', 'Op', 'Glue'
00298     
00299 template<typename T1, typename T2>
00300 inline
00301 typename T1::elem_type
00302 accu(const Glue<T1,T2,glue_schur>& X)
00303   {
00304   arma_extra_debug_sigprint();
00305 
00306   isnt_same_type<typename T1::elem_type, typename T2::elem_type>::check();
00307 
00308   typedef typename T1::elem_type eT;
00309 
00310   const u32 N_mat = 1 + depth_lhs< glue_schur, Glue<T1,T2,glue_schur> >::num;
00311   arma_extra_debug_print( arma_boost::format("N_mat = %d") % N_mat );
00312 
00313   if(N_mat == 2)
00314     {
00315     return accu_schur(Mat<eT>(X.A), Mat<eT>(X.B));
00316     }
00317   else
00318     {
00319     const Mat<eT>* ptrs[N_mat];
00320     bool            del[N_mat];
00321   
00322     mat_ptrs<glue_schur, Glue<T1,T2,glue_schur> >::get_ptrs(ptrs, del, X);
00323   
00324     for(u32 i=0; i<N_mat; ++i)  arma_extra_debug_print( arma_boost::format("ptrs[%d] = %x") % i % ptrs[i] );
00325     for(u32 i=0; i<N_mat; ++i)  arma_extra_debug_print( arma_boost::format(" del[%d] = %d") % i %  del[i] );
00326   
00327     const Mat<eT>& tmp_mat = *(ptrs[0]);
00328     
00329     for(u32 i=1; i<N_mat; ++i)
00330       {
00331       arma_debug_assert_same_size(tmp_mat, *(ptrs[i]), "accu()");
00332       }
00333     
00334     // const u32 n_rows = ptrs[0]->n_rows;
00335     // const u32 n_cols = ptrs[0]->n_cols;
00336     
00337     eT val = eT(0);
00338     
00339     const u32 n_elem = ptrs[0]->n_elem;
00340     
00341     for(u32 j=0; j<n_elem; ++j)
00342       {
00343       eT tmp = ptrs[0]->mem[j];
00344     
00345       for(u32 i=1; i<N_mat; ++i)
00346         {
00347         tmp *= ptrs[i]->mem[j];
00348         }
00349     
00350       val += tmp;
00351       }
00352     
00353     
00354     for(u32 i=0; i<N_mat; ++i)
00355       {
00356       if(del[i] == true)
00357         {
00358         arma_extra_debug_print( arma_boost::format("delete mat_ptr[%d]") % i );
00359         delete ptrs[i];
00360         }
00361       }
00362     
00363     return val;    
00364     }
00365   }
00366 
00367 
00368 //! accumulate the result of submatrix % matrix, where % is the Schur product (element-wise multiplication)
00369 template<typename eT>
00370 inline
00371 eT
00372 accu(const Glue<subview<eT>,Mat<eT>,glue_schur>& X)
00373   {
00374   arma_extra_debug_sigprint();
00375   
00376   arma_debug_assert_same_size(X.A, X.B, "accu()");
00377   
00378   const Mat<eT>& A = X.A.m;
00379   const Mat<eT>& B = X.B;
00380   
00381   const u32 A_sub_n_rows = X.A.n_rows;
00382   const u32 A_sub_n_cols = X.A.n_cols;
00383   
00384   const u32 A_aux_row1 = X.A.aux_row1;
00385   const u32 A_aux_col1 = X.A.aux_col1;
00386   
00387   
00388   eT val = eT(0);
00389     
00390   for(u32 col = 0; col<A_sub_n_cols; ++col)
00391     {
00392     const u32 col_mod = A_aux_col1 + col;
00393     
00394     for(u32 row = 0; row<A_sub_n_rows; ++row)
00395       {
00396       const u32 row_mod = A_aux_row1 + row;
00397       
00398       val += A.at(row_mod, col_mod) * B.at(row,col);
00399       }
00400     
00401     }
00402   
00403   return val;
00404   }
00405 
00406 
00407 
00408 //! accumulate the result of matrix % submatrix, where % is the Schur product (element-wise multiplication)
00409 template<typename eT>
00410 inline
00411 eT
00412 accu(const Glue<Mat<eT>,subview<eT>,glue_schur>& X)
00413   {
00414   arma_extra_debug_sigprint();
00415   
00416   arma_debug_assert_same_size(X.A, X.B, "accu()");
00417   
00418   const Mat<eT>& A = X.A;
00419   const Mat<eT>& B = X.B.m;
00420   
00421   // const u32 B_sub_n_rows = X.B.n_rows;
00422   // const u32 B_sub_n_cols = X.B.n_cols;
00423   
00424   const u32 B_aux_row1 = X.B.aux_row1;
00425   const u32 B_aux_col1 = X.B.aux_col1;
00426   
00427   
00428   eT val = eT(0);
00429     
00430   for(u32 col = 0; col<A.n_cols; ++col)
00431     {
00432     const u32 col_mod = B_aux_col1 + col;
00433     
00434     for(u32 row = 0; row<A.n_rows; ++row)
00435       {
00436       const u32 row_mod = B_aux_row1 + row;
00437       
00438       val += A.at(row, col) * B.at(row_mod, col_mod);
00439       }
00440     
00441     }
00442   
00443   return val;
00444   }
00445 
00446 
00447 
00448 //! accumulate the result of submatrix % submatrix, where % is the Schur product (element-wise multiplication)
00449 template<typename eT>
00450 inline
00451 eT
00452 accu(const Glue<subview<eT>,subview<eT>,glue_schur>& X)
00453   {
00454   arma_extra_debug_sigprint();
00455   
00456   arma_debug_assert_same_size(X.A, X.B, "accu()");
00457   
00458   const Mat<eT>& A = X.A.m;
00459   const Mat<eT>& B = X.B.m;
00460   
00461   const u32 A_sub_n_rows = X.A.n_rows;
00462   const u32 A_sub_n_cols = X.A.n_cols;
00463   
00464   // const u32 B_sub_n_rows = X.B.n_rows;
00465   // const u32 B_sub_n_cols = X.B.n_cols;
00466   
00467   const u32 A_aux_row1 = X.A.aux_row1;
00468   const u32 A_aux_col1 = X.A.aux_col1;
00469   
00470   const u32 B_aux_row1 = X.B.aux_row1;
00471   const u32 B_aux_col1 = X.B.aux_col1;
00472   
00473   
00474   eT val = eT(0);
00475     
00476   for(u32 col = 0; col<A_sub_n_cols; ++col)
00477     {
00478     const u32 A_col_mod = A_aux_col1 + col;
00479     const u32 B_col_mod = B_aux_col1 + col;
00480     
00481     for(u32 row = 0; row<A_sub_n_rows; ++row)
00482       {
00483       const u32 A_row_mod = A_aux_row1 + row;
00484       const u32 B_row_mod = B_aux_row1 + row;
00485       
00486       val += A.at(A_row_mod, A_col_mod) * B.at(B_row_mod, B_col_mod);
00487       }
00488     
00489     }
00490   
00491   return val;
00492   }
00493 
00494 
00495 
00496 //! @}