glue_plus_meat.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 glue_plus
00017 //! @{
00018 
00019 
00020 
00021 template<typename eT>
00022 inline
00023 void
00024 glue_plus::apply(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B)
00025   {
00026   arma_extra_debug_sigprint();
00027   
00028   arma_debug_assert_same_size(A, B, "matrix addition");
00029   
00030   // no aliasing problem
00031   out.set_size(A.n_rows, A.n_cols);
00032   
00033         eT* out_mem = out.memptr();
00034   const eT* A_mem   = A.mem;
00035   const eT* B_mem   = B.mem;
00036     
00037   const u32 n_elem  = out.n_elem;
00038   
00039   for(u32 i=0; i<n_elem; ++i)
00040     {
00041     out_mem[i] = A_mem[i] + B_mem[i];
00042     }
00043     
00044   }
00045 
00046 
00047 
00048 template<typename eT>
00049 inline
00050 void
00051 glue_plus::apply(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B, const Mat<eT>& C)
00052   {
00053   arma_extra_debug_sigprint();
00054   
00055   arma_debug_assert_same_size(A, B, "matrix addition");
00056   arma_debug_assert_same_size(A, C, "matrix addition");
00057   
00058   // no aliasing problem
00059   out.set_size(A.n_rows, A.n_cols);
00060     
00061         eT* out_mem = out.memptr();
00062   const eT* A_mem   = A.mem;
00063   const eT* B_mem   = B.mem;
00064   const eT* C_mem   = C.mem;
00065   
00066   const u32 n_elem  = A.n_elem;
00067   
00068   for(u32 i=0; i<n_elem; ++i)
00069     {
00070     out_mem[i] = A_mem[i] + B_mem[i] + C_mem[i];
00071     }
00072     
00073   }
00074 
00075 
00076 
00077 template<typename eT>
00078 inline
00079 void
00080 glue_plus::apply(Mat<eT>& out, const Glue<Mat<eT>,Mat<eT>,glue_plus>& X)
00081   {
00082   glue_plus::apply(out, X.A, X.B);
00083   }
00084 
00085 
00086 
00087 template<typename eT>
00088 inline
00089 void
00090 glue_plus::apply(Mat<eT>& out, const Glue< Glue<Mat<eT>,Mat<eT>,glue_plus>, Mat<eT>, glue_plus>& X)
00091   {
00092   glue_plus::apply(out, X.A.A, X.A.B, X.B);
00093   }
00094 
00095 
00096 
00097 template<typename T1, typename T2>
00098 inline
00099 void
00100 glue_plus::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_plus>& X)
00101   {
00102   arma_extra_debug_sigprint();
00103   
00104   typedef typename T1::elem_type eT;
00105   
00106   const u32 N_mat = 1 + depth_lhs< glue_plus, Glue<T1,T2,glue_plus> >::num;
00107   arma_extra_debug_print( arma_boost::format("N_mat = %d") % N_mat );
00108 
00109   if(N_mat == 2)
00110     {
00111     if(is_glue_times<T1>::value == true)
00112       {
00113       out = X.B;
00114       glue_plus::apply_inplace(out, X.A);
00115       }
00116     else
00117     if(is_glue_times<T2>::value == true)
00118       {
00119       out = X.A;
00120       glue_plus::apply_inplace(out, X.B);
00121       }
00122     else
00123       {
00124       const unwrap<T1> tmp1(X.A);
00125       const unwrap<T2> tmp2(X.B);
00126       
00127       glue_plus::apply(out, tmp1.M, tmp2.M);
00128       }
00129     }
00130   else
00131     {
00132     const Mat<eT>* ptrs[N_mat];
00133     bool            del[N_mat];
00134 
00135     mat_ptrs<glue_plus, Glue<T1,T2,glue_plus> >::get_ptrs(ptrs, del, X);
00136 
00137     for(u32 i=0; i<N_mat; ++i)  arma_extra_debug_print( arma_boost::format("ptrs[%d] = %x") % i % ptrs[i] );
00138     for(u32 i=0; i<N_mat; ++i)  arma_extra_debug_print( arma_boost::format(" del[%d] = %d") % i %  del[i] );
00139 
00140     const u32 n_rows = ptrs[0]->n_rows;
00141     const u32 n_cols = ptrs[0]->n_cols;
00142   
00143     const Mat<eT>& tmp_mat = *(ptrs[0]);
00144     
00145     for(u32 i=1; i<N_mat; ++i)
00146       {
00147       arma_debug_assert_same_size(tmp_mat, *(ptrs[i]), "matrix addition");
00148       }
00149   
00150   
00151     // no aliasing problem
00152     out.set_size(n_rows, n_cols);
00153     
00154     const u32 n_elem = ptrs[0]->n_elem;
00155     
00156     for(u32 j=0; j<n_elem; ++j)
00157       {
00158       eT acc = ptrs[0]->mem[j];
00159     
00160       for(u32 i=1; i < N_mat; ++i)
00161         {
00162         acc += ptrs[i]->mem[j];
00163         }
00164     
00165       out[j] = acc;
00166       }
00167     
00168     
00169     for(u32 i=0; i<N_mat; ++i)
00170       {
00171       if(del[i] == true)
00172         {
00173         arma_extra_debug_print( arma_boost::format("delete mat_ptr[%d]") % i );
00174         delete ptrs[i];
00175         }
00176       }
00177     }
00178   }
00179 
00180 
00181 
00182 // possible aliasing cases:
00183 // Q = Q + Q.row(0)  -> no problem  (aliasing has no effect or incompatible matrix dimensions).
00184 //                      however, the only time the above will work is when Q has the same dimensions as Q.row(0),
00185 //                      meaning that doing this addition operation is pretty silly
00186 // Q = Q + R.row(0)  -> no problem
00187 // Q = R + Q.row(0)  -> output Q is set to size of R, which may destroy input Q
00188 //
00189 // strategy:
00190 // if the matrix from the second argument is an alias of the output matrix,
00191 // make a proper matrix out of the second argument
00192 
00193 template<typename eT>
00194 inline
00195 void
00196 glue_plus::apply(Mat<eT>& out, const Glue<Mat<eT>, subview<eT>, glue_plus>& X)
00197   {
00198   arma_extra_debug_sigprint();
00199   
00200   const Mat<eT>& orig_A = X.A;
00201   const Mat<eT>& orig_B = X.B.m;
00202   
00203   if( &out != &orig_B )
00204     {
00205     //const u32 sub_B_n_rows = X.B.n_rows;
00206     //const u32 sub_B_n_cols = X.B.n_cols;
00207     
00208     arma_debug_assert_same_size(X.A, X.B, "matrix addition");
00209       
00210     
00211     out.set_size(orig_A.n_rows, orig_A.n_cols);
00212     
00213     for(u32 col = 0; col<orig_A.n_cols; ++col)
00214       {
00215       const u32 B_col_mod = X.B.aux_col1 + col;
00216       
00217       for(u32 row = 0; row<orig_A.n_rows; ++row)
00218         {
00219         const u32 B_row_mod = X.B.aux_row1 + row;
00220         
00221         out.at(row,col) =  orig_A.at(row, col) + orig_B.at(B_row_mod, B_col_mod);
00222         }
00223       
00224       }
00225     
00226     }
00227   else
00228     {
00229     const Mat<eT> processed_B(X.B);  // create a matrix out of subview
00230     glue_plus::apply(out, orig_A, processed_B);
00231     }
00232      
00233   }
00234 
00235 
00236 // possible aliasing cases:
00237 // Q = Q.row(0) + Q  -> no problem (aliasing has no effect or incompatible matrix dimensions)
00238 // Q = Q.row(0) + R  -> problem (output Q is set to size of Q.row(0) which may destroy input Q)
00239 // Q = R.row(0) + Q  -> problem (output Q is set to size of R.row(0) which may destroy input Q)
00240 
00241 template<typename eT>
00242 inline
00243 void
00244 glue_plus::apply(Mat<eT>& out, const Glue< subview<eT>, Mat<eT>, glue_plus>& X)
00245   {
00246   arma_extra_debug_sigprint();
00247   
00248   const Mat<eT>& orig_A = X.A.m;
00249   
00250   const unwrap_check< Mat<eT> > tmp(X.B, out);
00251   const Mat<eT>& orig_B = tmp.M;
00252   
00253   if( &out != &orig_A )
00254     {
00255     const u32 sub_A_n_rows = X.A.n_rows;
00256     const u32 sub_A_n_cols = X.A.n_cols;
00257     
00258     arma_debug_assert_same_size(X.A, X.B, "matrix addition");
00259       
00260     out.set_size(sub_A_n_rows, sub_A_n_cols);
00261     
00262     for(u32 col = 0; col<sub_A_n_cols; ++col)
00263       {
00264       const u32 A_col_mod = X.A.aux_col1 + col;
00265       
00266       for(u32 row = 0; row<sub_A_n_rows; ++row)
00267         {
00268         const u32 A_row_mod = X.A.aux_row1 + row;
00269         
00270         out.at(row,col) =  orig_A.at(A_row_mod, A_col_mod) + orig_B.at(row, col);
00271         }
00272       
00273       }
00274     }
00275   else
00276     {
00277     const Mat<eT> processed_A(X.A);
00278     glue_plus::apply(out, processed_A, orig_B);
00279     }
00280   
00281   
00282   }
00283 
00284 
00285 // possible aliasing cases:
00286 // Q = Q.row(0) + Q.row(0)  -> input Q is destroyed unless Q.row(0) has the same size as Q 
00287 // Q = Q.row(0) + R.row(0)  -> input Q is destroyed unless Q.row(0) has the same size as Q 
00288 // Q = R.row(0) + Q.row(0)  -> input Q is destroyed
00289 
00290 template<typename eT>
00291 inline
00292 void
00293 glue_plus::apply(Mat<eT>& out, const Glue< subview<eT>, subview<eT>, glue_plus>& X)
00294   {
00295   arma_extra_debug_sigprint();
00296   
00297   const Mat<eT>& orig_A = X.A.m;
00298   const Mat<eT>& orig_B = X.B.m;
00299   
00300   if( (&out != &orig_A) && (&out != &orig_B) )
00301     {
00302     const u32 sub_A_n_rows = X.A.n_rows;
00303     const u32 sub_A_n_cols = X.A.n_cols;
00304     
00305     const u32 sub_B_n_rows = X.B.n_rows;
00306     const u32 sub_B_n_cols = X.B.n_cols;
00307     
00308     arma_debug_assert_same_size(X.A, X.B, "matrix addition");
00309       
00310     out.set_size(sub_A_n_rows, sub_A_n_cols);
00311     
00312     for(u32 col = 0; col<sub_A_n_cols; ++col)
00313       {
00314       const u32 A_col_mod = X.A.aux_col1 + col;
00315       const u32 B_col_mod = X.B.aux_col1 + col;
00316       
00317       for(u32 row = 0; row<sub_A_n_rows; ++row)
00318         {
00319         const u32 A_row_mod = X.A.aux_row1 + row;
00320         const u32 B_row_mod = X.B.aux_row1 + row;
00321         
00322         out.at(row,col) =  orig_A.at(A_row_mod, A_col_mod) + orig_B.at(B_row_mod, B_col_mod);
00323         }
00324       
00325       }
00326     }
00327   else
00328     {
00329     const Mat<eT> processed_A(X.A);
00330     const Mat<eT> processed_B(X.B);
00331     
00332     glue_plus::apply(out, processed_A, processed_B);
00333     }
00334   }
00335 
00336 
00337 
00338 template<typename eT>
00339 inline void glue_plus::apply_inplace(Mat<eT>& out, const Mat<eT>& B)
00340   {
00341   arma_extra_debug_sigprint();
00342   
00343   arma_debug_assert_same_size(out, B, "matrix addition");
00344   
00345   
00346         eT* out_mem = out.memptr();
00347   const eT* B_mem   = B.mem;
00348   
00349   const u32 n_elem  = B.n_elem;
00350   
00351   for(u32 i=0; i<n_elem; ++i)
00352     {
00353     out_mem[i] += B_mem[i];
00354     }
00355   
00356   }
00357 
00358 
00359 
00360 template<typename T1, typename op_type>
00361 inline
00362 void
00363 glue_plus::apply_inplace(Mat<typename T1::elem_type>& out, const Op<T1, op_type>& X)
00364   {
00365   arma_extra_debug_sigprint();
00366   
00367   typedef typename T1::elem_type eT;
00368   
00369   const Mat<eT> tmp(X);
00370   glue_plus::apply(out, out, tmp);
00371   }
00372   
00373 
00374 
00375 template<typename T1>
00376 inline
00377 void
00378 glue_plus::apply_inplace(Mat<typename T1::elem_type>& out, const Op<T1, op_square>& X)
00379   {
00380   arma_extra_debug_sigprint();
00381   
00382   typedef typename T1::elem_type eT;
00383   
00384   const unwrap<T1> tmp(X.m);
00385   const Mat<eT>& B = tmp.M;
00386   
00387   arma_debug_assert_same_size(out, B, "matrix addition");
00388     
00389         eT* out_mem = out.memptr();
00390   const eT* B_mem   = B.mem;
00391   
00392   const u32 n_elem  = out.n_elem;
00393   
00394   for(u32 i=0; i<n_elem; ++i)
00395     {
00396     const eT tmp = B_mem[i];
00397     out_mem[i] += tmp*tmp;
00398     }
00399   }
00400 
00401 
00402 
00403 template<typename T1, typename T2, typename glue_type>
00404 inline
00405 void
00406 glue_plus::apply_inplace(Mat<typename T1::elem_type>& out, const Glue<T1, T2, glue_type>& X)
00407   {
00408   arma_extra_debug_sigprint();
00409     
00410   out = X + out;
00411   }
00412 
00413 
00414 
00415 template<typename T1, typename T2>
00416 inline
00417 void
00418 glue_plus::apply_inplace(Mat<typename T1::elem_type>& out, const Glue<T1, T2, glue_times>& X)
00419   {
00420   arma_extra_debug_sigprint();
00421   
00422   typedef typename T1::elem_type eT;
00423   
00424   const unwrap_check<T1> tmp1(X.A, out);
00425   const unwrap_check<T2> tmp2(X.B, out);
00426   
00427   const Mat<eT>& A = tmp1.M;
00428   const Mat<eT>& B = tmp2.M;
00429   
00430   arma_debug_assert_mul_size(A, B, "matrix multiplication");
00431   arma_debug_assert_same_size(out.n_rows, out.n_cols, A.n_rows, B.n_cols, "matrix addition");
00432   
00433   gemm<false,false,false,true>::apply(out, A, B, eT(1), eT(1));
00434   }
00435 
00436 
00437 
00438 //
00439 //
00440 //
00441 
00442 
00443 
00444 template<typename T1, typename T2>
00445 inline
00446 void
00447 glue_plus::apply
00448   (
00449   Mat<typename T1::elem_type>& out,
00450   const Glue< Glue< T1, Col<typename T1::elem_type>, glue_times_vec>, T2, glue_plus>& in
00451   )
00452   {
00453   arma_extra_debug_sigprint();
00454   
00455   typedef typename T1::elem_type eT;
00456   
00457   const unwrap<T1> tmp1(in.A.A);
00458   const unwrap<T2> tmp2(in.B);
00459   
00460   const Mat<eT>& A = tmp1.M;
00461   const Col<eT>& B = in.A.B;
00462   const Mat<eT>& C = tmp2.M;
00463   
00464   arma_debug_assert_mul_size(A, B, "matrix multiplication");
00465   arma_debug_assert_same_size(A.n_rows, B.n_cols, C.n_rows, C.n_cols, "matrix addition");
00466   
00467   if( (&out != &A) && (&out != &B) )
00468     {
00469     out = C;
00470     gemv<false,false,true>::apply(out.memptr(), A, B.mem, eT(1), eT(1));
00471     }
00472   else
00473     {
00474     const unwrap_check< Mat<eT>    > tmp1(A,out);
00475     const unwrap_check< Col<eT> > tmp2(B,out);
00476     
00477     const Mat<eT>& A_safe = tmp1.M;
00478     const Col<eT>& B_safe = tmp2.M;
00479     
00480     out = C;
00481     gemv<false,false,true>::apply(out.memptr(), A_safe, B_safe.mem, eT(1), eT(1));
00482     }
00483   
00484   }
00485 
00486 
00487 
00488 template<typename T1, typename T2>
00489 inline
00490 void
00491 glue_plus::apply
00492   (
00493   Mat<typename T1::elem_type>& out,
00494   const Glue< Glue< Row<typename T1::elem_type>, T1, glue_times_vec>, T2, glue_plus>& in
00495   )
00496   {
00497   arma_extra_debug_sigprint();
00498   
00499   typedef typename T1::elem_type eT;
00500   
00501   const unwrap<T1> tmp1(in.A.B);
00502   const unwrap<T2> tmp2(in.B);
00503   
00504   const Row<eT>& A = in.A.A;
00505   const Mat<eT>& B = tmp1.M;
00506   const Mat<eT>& C = tmp2.M;
00507   
00508   arma_debug_assert_mul_size(A, B, "matrix multiplication");
00509   arma_debug_assert_same_size(A.n_rows, B.n_cols, C.n_rows, C.n_cols, "matrix addition");
00510   
00511   if( (&out != &A) && (&out != &B) )
00512     {
00513     out = C;
00514     gemv<true,false,true>::apply(out.memptr(), B, A.mem, eT(1), eT(1));
00515     }
00516   else
00517     {
00518     const unwrap_check< Mat<eT> > tmp1(A,out);
00519     const unwrap_check< Mat<eT> > tmp2(B,out);
00520     
00521     const Mat<eT>& A_safe = tmp1.M;
00522     const Mat<eT>& B_safe = tmp2.M;
00523     
00524     out = C;
00525     gemv<true,false,true>::apply(out.memptr(), B_safe, A_safe.mem, eT(1), eT(1));
00526     }
00527   
00528   }
00529 
00530 
00531 
00532 //
00533 // matrix addition with different element types
00534 
00535 template<typename eT1, typename eT2>
00536 inline
00537 void
00538 glue_plus::apply_mixed(Mat<typename promote_type<eT1,eT2>::result>& out, const Mat<eT1>& X, const Mat<eT2>& Y)
00539   {
00540   arma_extra_debug_sigprint();
00541   
00542   typedef typename promote_type<eT1,eT2>::result out_eT;
00543   
00544   arma_debug_assert_same_size(X,Y, "matrix addition");  
00545   
00546   out.set_size(X.n_rows, X.n_cols);
00547   
00548         out_eT* out_mem = out.memptr();
00549   const eT1*    X_mem   = X.mem;
00550   const eT2*    Y_mem   = Y.mem;
00551   
00552   const u32 n_elem = out.n_elem;
00553   
00554   for(u32 i=0; i<n_elem; ++i)
00555     {
00556     out_mem[i] = upgrade_val<eT1,eT2>::apply(X_mem[i]) + upgrade_val<eT1,eT2>::apply(Y_mem[i]);
00557     }
00558   }
00559 
00560 
00561 
00562 //
00563 // glue_plus_diag
00564 
00565 
00566 template<typename T1, typename T2>
00567 inline
00568 void
00569 glue_plus_diag::apply(Mat<typename T1::elem_type>& out, const T1& A_orig, const Op<T2,op_diagmat>& B_orig)
00570   {
00571   arma_extra_debug_sigprint();
00572   
00573   isnt_same_type<typename T1::elem_type, typename T2::elem_type>::check();
00574   
00575   typedef typename T1::elem_type eT;
00576   
00577   const unwrap<T1> tmp1(A_orig);
00578   const unwrap<T2> tmp2(B_orig.m);
00579   
00580   const Mat<eT>& A = tmp1.M;  
00581   const Mat<eT>& B = tmp2.M;
00582   
00583   arma_debug_check( !B.is_square(), "glue_plus_diag::apply(): matrices must be square" );
00584   arma_debug_assert_same_size(A, B, "matrix addition");
00585 
00586   
00587   // no aliasing problem
00588   out.set_size(A.n_rows, A.n_cols);
00589   
00590   for(u32 col=0; col<A.n_cols; ++col)
00591     {
00592     for(u32 row=0; row<A.n_rows; ++row)
00593       {
00594       if(col != row)
00595         {
00596         out.at(row,col) = A.at(row,col);
00597         }
00598       else
00599         {
00600         out.at(row,col) = A.at(row,col) + B.at(row,col);
00601         }
00602       }
00603     }
00604   
00605   }
00606 
00607 
00608 
00609 template<typename T1, typename T2>
00610 inline
00611 void
00612 glue_plus_diag::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_diagmat>& A_orig, const Op<T2,op_diagmat>& B_orig)
00613   {
00614   arma_extra_debug_sigprint();
00615   
00616   isnt_same_type<typename T1::elem_type, typename T2::elem_type>::check();
00617   
00618   const unwrap<T1> tmp1(A_orig.m);
00619   const unwrap<T2> tmp2(B_orig.m);
00620   
00621   typedef typename T1::elem_type eT;
00622   
00623   const Mat<eT>& A = tmp1.M;
00624   const Mat<eT>& B = tmp2.M;
00625   
00626   arma_debug_check( !A.is_square(), "glue_plus_diag::apply(): matrices must be square" );
00627   arma_debug_assert_same_size(A, B, "matrix addition");
00628   
00629   
00630   if( (&out != &A) && (&out != &B) )
00631     {
00632     out.zeros(A.n_rows, A.n_cols);
00633     
00634     for(u32 i=0; i<A.n_rows; ++i)
00635       {
00636       out.at(i,i) = A.at(i,i) + B.at(i,i);
00637       }
00638     }
00639   else
00640     {
00641     out.set_size(A.n_rows, A.n_cols);
00642   
00643     for(u32 col=0; col<A.n_cols; ++col)
00644       {
00645       for(u32 row=0; row<A.n_rows; ++row)
00646         {
00647         if(col != row)
00648           {
00649           out.at(row,col) = 0.0;
00650           }
00651         else
00652           {
00653           out.at(row,col) = A.at(row,col) + B.at(row,col);
00654           }
00655         }
00656       }
00657     }
00658   
00659   }
00660 
00661 
00662 
00663 template<typename T1, typename T2>
00664 inline
00665 void
00666 glue_plus_diag::apply(Mat<typename T1::elem_type>& out, const Glue<T1, Op<T2,op_diagmat>, glue_plus_diag>& X)
00667   {
00668   glue_plus_diag::apply(out, X.A, X.B);
00669   }
00670 
00671 
00672 
00673 template<typename T1, typename T2>
00674 inline
00675 void
00676 glue_plus_diag::apply(Mat<typename T1::elem_type>& out, const Glue<Op<T1,op_diagmat>, T2, glue_plus_diag>& X)
00677   {
00678   glue_plus_diag::apply(out, X.B, X.A);  // NOTE: arguments are swapped
00679   }
00680 
00681 
00682 
00683 template<typename T1, typename T2>
00684 inline
00685 void
00686 glue_plus_diag::apply(Mat<typename T1::elem_type>& out, const Glue<Op<T1,op_diagmat>, Op<T2,op_diagmat>, glue_plus_diag>& X)
00687   {
00688   glue_plus_diag::apply(out, X.A, X.B);
00689   }
00690 
00691 
00692 
00693 //! @}