eop_core_meat.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 eop_core
00018 //! @{
00019 
00020 
00021 
00022 template<typename eop_type>
00023 template<typename T1>
00024 arma_hot
00025 arma_inline
00026 typename T1::elem_type
00027 eop_core<eop_type>::get_elem(const eOp<T1, eop_type>& x, const u32 i)
00028   {
00029   typedef typename T1::elem_type eT;
00030   
00031        if(is_generator<eop_type>::value                == true) { return eop_aux::generate<eT,eop_type>();                       }
00032   else if(is_same_type<eop_type, eop_ones_diag>::value == true) { return ((i % x.P.n_rows) == (i / x.P.n_rows)) ? eT(1) : eT(0); }
00033   else                                                          { return eop_core<eop_type>::process(x, x.P[i]);                 }
00034   }
00035 
00036 
00037 
00038 template<typename eop_type>
00039 template<typename T1>
00040 arma_hot
00041 arma_inline
00042 typename T1::elem_type
00043 eop_core<eop_type>::get_elem(const eOp<T1, eop_type>& x, const u32 row, const u32 col)
00044   {
00045   typedef typename T1::elem_type eT;
00046   
00047        if(is_generator<eop_type>::value                == true) { return eop_aux::generate<eT,eop_type>();                 }
00048   else if(is_same_type<eop_type, eop_ones_diag>::value == true) { return (row == col) ? eT(1) : eT(0);                     }
00049   else                                                          { return eop_core<eop_type>::process(x, x.P.at(row, col)); }
00050   }
00051 
00052 
00053 
00054 template<typename eop_type>
00055 template<typename T1>
00056 arma_hot
00057 arma_inline
00058 typename T1::elem_type
00059 eop_core<eop_type>::process(const eOp<T1, eop_type>& x, const typename T1::elem_type val)
00060   {
00061   typedef typename T1::elem_type eT;
00062   
00063   // the optimiser will keep only one return statement
00064   
00065        if(is_same_type<eop_type, eop_neg              >::value == true) { return -val;                     }
00066   else if(is_same_type<eop_type, eop_scalar_plus      >::value == true) { return val + x.aux;              }
00067   else if(is_same_type<eop_type, eop_scalar_minus_pre >::value == true) { return x.aux - val;              }
00068   else if(is_same_type<eop_type, eop_scalar_minus_post>::value == true) { return val - x.aux;              }
00069   else if(is_same_type<eop_type, eop_scalar_times     >::value == true) { return val * x.aux;              }
00070   else if(is_same_type<eop_type, eop_scalar_div_pre   >::value == true) { return x.aux / val;              }
00071   else if(is_same_type<eop_type, eop_scalar_div_post  >::value == true) { return val / x.aux;              }
00072   else if(is_same_type<eop_type, eop_square           >::value == true) { return val*val;                  }
00073   else if(is_same_type<eop_type, eop_sqrt             >::value == true) { return eop_aux::sqrt(val);       }
00074   else if(is_same_type<eop_type, eop_log10            >::value == true) { return eop_aux::log10(val);      }
00075   else if(is_same_type<eop_type, eop_log              >::value == true) { return eop_aux::log(val);        }
00076   else if(is_same_type<eop_type, eop_trunc_log        >::value == true) { return    arma::trunc_log(val);  }
00077   else if(is_same_type<eop_type, eop_exp              >::value == true) { return eop_aux::exp(val);        }
00078   else if(is_same_type<eop_type, eop_trunc_exp        >::value == true) { return    arma::trunc_exp(val);  }
00079   else if(is_same_type<eop_type, eop_cos              >::value == true) { return eop_aux::cos(val);        }
00080   else if(is_same_type<eop_type, eop_sin              >::value == true) { return eop_aux::sin(val);        }
00081   else if(is_same_type<eop_type, eop_tan              >::value == true) { return eop_aux::tan(val);        }
00082   else if(is_same_type<eop_type, eop_acos             >::value == true) { return eop_aux::acos(val);       }
00083   else if(is_same_type<eop_type, eop_asin             >::value == true) { return eop_aux::asin(val);       }
00084   else if(is_same_type<eop_type, eop_atan             >::value == true) { return eop_aux::atan(val);       }
00085   else if(is_same_type<eop_type, eop_cosh             >::value == true) { return eop_aux::cosh(val);       }
00086   else if(is_same_type<eop_type, eop_sinh             >::value == true) { return eop_aux::sinh(val);       }
00087   else if(is_same_type<eop_type, eop_tanh             >::value == true) { return eop_aux::tanh(val);       }
00088   else if(is_same_type<eop_type, eop_acosh            >::value == true) { return eop_aux::acosh(val);      }
00089   else if(is_same_type<eop_type, eop_asinh            >::value == true) { return eop_aux::asinh(val);      }
00090   else if(is_same_type<eop_type, eop_atanh            >::value == true) { return eop_aux::atanh(val);      }
00091   else if(is_same_type<eop_type, eop_eps              >::value == true) { return eop_aux::direct_eps(val); }
00092   else if(is_same_type<eop_type, eop_abs              >::value == true) { return eop_aux::arma_abs(val);   }
00093   else if(is_same_type<eop_type, eop_conj             >::value == true) { return eop_aux::conj(val);       }
00094   else if(is_same_type<eop_type, eop_pow              >::value == true) { return eop_aux::pow(val, x.aux); }
00095   else if(is_same_type<eop_type, eop_pow_int          >::value == true)
00096     {
00097     const int exponent = (x.aux_u32_b == 0) ? int(x.aux_u32_a) : -int(x.aux_u32_a);
00098     
00099     return eop_aux::pow_int(val, exponent);
00100     }
00101   else
00102     {
00103     arma_stop("eop_core::process(): unhandled eop_type");
00104     return eT(0);
00105     }
00106   }
00107 
00108 
00109 
00110 template<typename eop_type>
00111 template<typename T1>
00112 arma_hot
00113 arma_inline
00114 void
00115 eop_core<eop_type>::apply(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00116   {
00117   arma_extra_debug_sigprint();
00118   
00119   if(is_Mat<T1>::value == true)
00120     {
00121     eop_core<eop_type>::apply_unwrap(out, x);
00122     }
00123   else
00124     {
00125     eop_core<eop_type>::apply_proxy(out, x);
00126     }
00127   }
00128 
00129 
00130 
00131 template<typename eop_type>
00132 template<typename T1>
00133 arma_hot
00134 inline
00135 void
00136 eop_core<eop_type>::apply_proxy(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00137   {
00138   arma_extra_debug_sigprint();
00139   
00140   // eop_type::apply_proxy() function is not allowed to unwrap things
00141   // (in order to get the input into a common format).
00142   // the proxy class is already providing objects with element access
00143   
00144   typedef typename T1::elem_type eT;
00145   
00146   const Proxy<T1>& P = x.P;
00147   
00148   out.set_size(P.n_rows, P.n_cols);
00149   
00150         eT* out_mem = out.memptr();
00151   const u32 n_elem  = P.n_elem;
00152     
00153   if(is_generator<eop_type>::value == true)
00154     {
00155     for(u32 i=0; i<n_elem; ++i)
00156       {
00157       out_mem[i] = eop_aux::generate<eT,eop_type>();
00158       }
00159     }
00160   else
00161     {
00162     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00163       {
00164       for(u32 col=0; col<P.n_rows; ++col)
00165         {
00166         for(u32 row=0; row<col; ++row)          { out.at(row,col) = eT(0); }
00167           
00168         out.at(col,col) = eT(1);
00169         
00170         for(u32 row=col+1; row<P.n_rows; ++row) { out.at(row,col) = eT(0); }
00171         }
00172       }
00173     else
00174       {
00175       for(u32 i=0; i<n_elem; ++i)
00176         {
00177         out_mem[i] = eop_core<eop_type>::process(x, P[i]);
00178         }
00179       }
00180     }
00181   }
00182 
00183 
00184 
00185 template<typename eop_type>
00186 template<typename T1>
00187 arma_hot
00188 inline
00189 void
00190 eop_core<eop_type>::apply_unwrap(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00191   {
00192   arma_extra_debug_sigprint();
00193   
00194   typedef typename T1::elem_type eT;
00195   
00196   const Proxy<T1>& P = x.P;
00197   
00198 //   cout << "*** P.n_rows = " << P.n_rows << endl;
00199 //   cout << "*** P.n_cols = " << P.n_cols << endl;
00200   
00201   out.set_size(P.n_rows, P.n_cols);
00202   
00203         eT* out_mem = out.memptr();
00204   const u32 n_elem  = P.n_elem;
00205     
00206   const unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
00207   const Mat<eT>& A = tmp.M;
00208   
00209   if(is_generator<eop_type>::value == true)
00210     {
00211     for(u32 i=0; i<n_elem; ++i)
00212       {
00213       out_mem[i] = eop_aux::generate<eT,eop_type>();
00214       }
00215     }
00216   else
00217     {
00218     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00219       {
00220       for(u32 col=0; col<P.n_rows; ++col)
00221         {
00222         for(u32 row=0; row<col; ++row)          { out.at(row,col) = eT(0); }
00223           
00224         out.at(col,col) = eT(1);
00225         
00226         for(u32 row=col+1; row<P.n_rows; ++row) { out.at(row,col) = eT(0); }
00227         }
00228       }
00229     else
00230       {
00231       const eT* A_mem = A.memptr();
00232       
00233       for(u32 i=0; i<n_elem; ++i)
00234         {
00235         out_mem[i] = eop_core<eop_type>::process(x, A_mem[i]);
00236         }
00237       }
00238     }
00239   }
00240 
00241 
00242 
00243 template<typename eop_type>
00244 template<typename T1>
00245 arma_hot
00246 inline
00247 void
00248 eop_core<eop_type>::apply_inplace_plus(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00249   {
00250   arma_extra_debug_sigprint();
00251   
00252   typedef typename T1::elem_type eT;
00253   
00254   const Proxy<T1>& P = x.P;
00255   
00256   arma_debug_assert_same_size(out.n_rows, out.n_cols, P.n_rows, P.n_cols, "matrix addition");
00257   
00258         eT* out_mem = out.memptr();
00259   const u32 n_elem  = P.n_elem;
00260     
00261   if(is_generator<eop_type>::value == true)
00262     {
00263     for(u32 i=0; i<n_elem; ++i)
00264       {
00265       out_mem[i] += eop_aux::generate<eT,eop_type>();
00266       }
00267     }
00268   else
00269     {
00270     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00271       {
00272       for(u32 row=0; row<P.n_rows; ++row)
00273         {
00274         out.at(row,row) += eT(1);
00275         }
00276       }
00277     else
00278       {
00279       for(u32 i=0; i<n_elem; ++i)
00280         {
00281         out_mem[i] += eop_core<eop_type>::process(x, P[i]);
00282         }
00283       }
00284     }
00285   }
00286 
00287 
00288 
00289 template<typename eop_type>
00290 template<typename T1>
00291 arma_hot
00292 inline
00293 void
00294 eop_core<eop_type>::apply_inplace_minus(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00295   {
00296   arma_extra_debug_sigprint();
00297   
00298   typedef typename T1::elem_type eT;
00299   
00300   const Proxy<T1>& P = x.P;
00301   
00302   arma_debug_assert_same_size(out.n_rows, out.n_cols, P.n_rows, P.n_cols, "matrix subtraction");
00303   
00304         eT* out_mem = out.memptr();
00305   const u32 n_elem  = P.n_elem;
00306     
00307   if(is_generator<eop_type>::value == true)
00308     {
00309     for(u32 i=0; i<n_elem; ++i)
00310       {
00311       out_mem[i] -= eop_aux::generate<eT,eop_type>();
00312       }
00313     }
00314   else
00315     {
00316     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00317       {
00318       for(u32 row=0; row<P.n_rows; ++row)
00319         {
00320         out.at(row,row) -= eT(1);
00321         }
00322       }
00323     else
00324       {
00325       for(u32 i=0; i<n_elem; ++i)
00326         {
00327         out_mem[i] -= eop_core<eop_type>::process(x, P[i]);
00328         }
00329       }
00330     }
00331   }
00332 
00333 
00334 
00335 template<typename eop_type>
00336 template<typename T1>
00337 arma_hot
00338 inline
00339 void
00340 eop_core<eop_type>::apply_inplace_schur(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00341   {
00342   arma_extra_debug_sigprint();
00343   
00344   typedef typename T1::elem_type eT;
00345   
00346   const Proxy<T1>& P = x.P;
00347   
00348   arma_debug_assert_same_size(out.n_rows, out.n_cols, P.n_rows, P.n_cols, "element-wise matrix multiplication");
00349   
00350         eT* out_mem = out.memptr();
00351   const u32 n_elem  = P.n_elem;
00352   
00353   if(is_generator<eop_type>::value == true)
00354     {
00355     for(u32 i=0; i<n_elem; ++i)
00356       {
00357       out_mem[i] *= eop_aux::generate<eT,eop_type>();
00358       }
00359     }
00360   else
00361     {
00362     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00363       {
00364       for(u32 col=0; col<P.n_rows; ++col)
00365         {
00366         for(u32 row=0;     row<col;      ++row) { out.at(row,col) = eT(0); }
00367         for(u32 row=col+1; row<P.n_rows; ++row) { out.at(row,col) = eT(0); }
00368         }
00369       }
00370     else
00371       {
00372       for(u32 i=0; i<n_elem; ++i)
00373         {
00374         out_mem[i] *= eop_core<eop_type>::process(x, P[i]);
00375         }
00376       }
00377     }
00378   }
00379 
00380 
00381 
00382 template<typename eop_type>
00383 template<typename T1>
00384 arma_hot
00385 inline
00386 void
00387 eop_core<eop_type>::apply_inplace_div(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00388   {
00389   arma_extra_debug_sigprint();
00390   
00391   typedef typename T1::elem_type eT;
00392   
00393   const Proxy<T1>& P = x.P;
00394   
00395   arma_debug_assert_same_size(out.n_rows, out.n_cols, P.n_rows, P.n_cols, "element-wise matrix division");
00396   
00397         eT* out_mem = out.memptr();
00398   const u32 n_elem  = P.n_elem;
00399   
00400   if(is_generator<eop_type>::value == true)
00401     {
00402     for(u32 i=0; i<n_elem; ++i)
00403       {
00404       out_mem[i] /= eop_aux::generate<eT,eop_type>();
00405       }
00406     }
00407   else
00408     {
00409     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00410       {
00411       for(u32 col=0; col<P.n_rows; ++col)
00412         {
00413         for(u32 row=0;     row<col;      ++row) { out.at(row,col) /= eT(0); }
00414         for(u32 row=col+1; row<P.n_rows; ++row) { out.at(row,col) /= eT(0); }
00415         }
00416       }
00417     else
00418       {
00419       for(u32 i=0; i<n_elem; ++i)
00420         {
00421         out_mem[i] /= eop_core<eop_type>::process(x, P[i]);
00422         }
00423       }
00424     }
00425   }
00426 
00427 
00428 
00429 //! @}