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 std::sqrt(val);           }
00074   else if(is_same_type<eop_type, eop_log10            >::value == true) { return std::log10(val);          }
00075   else if(is_same_type<eop_type, eop_log              >::value == true) { return std::log(val);            }
00076   else if(is_same_type<eop_type, eop_trunc_log        >::value == true) { return eop_aux::trunc_log(val);  }
00077   else if(is_same_type<eop_type, eop_exp              >::value == true) { return std::exp(val);            }
00078   else if(is_same_type<eop_type, eop_trunc_exp        >::value == true) { return eop_aux::trunc_exp(val);  }
00079   else if(is_same_type<eop_type, eop_cos              >::value == true) { return std::cos(val);            }
00080   else if(is_same_type<eop_type, eop_sin              >::value == true) { return std::sin(val);            }
00081   else if(is_same_type<eop_type, eop_tan              >::value == true) { return std::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 std::cosh(val);           }
00086   else if(is_same_type<eop_type, eop_sinh             >::value == true) { return std::sinh(val);           }
00087   else if(is_same_type<eop_type, eop_tanh             >::value == true) { return std::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   out.set_size(P.n_rows, P.n_cols);
00199   
00200         eT* out_mem = out.memptr();
00201   const u32 n_elem  = P.n_elem;
00202     
00203   const unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
00204   const Mat<eT>& A = tmp.M;
00205   
00206   if(is_generator<eop_type>::value == true)
00207     {
00208     for(u32 i=0; i<n_elem; ++i)
00209       {
00210       out_mem[i] = eop_aux::generate<eT,eop_type>();
00211       }
00212     }
00213   else
00214     {
00215     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00216       {
00217       for(u32 col=0; col<P.n_rows; ++col)
00218         {
00219         for(u32 row=0; row<col; ++row)          { out.at(row,col) = eT(0); }
00220           
00221         out.at(col,col) = eT(1);
00222         
00223         for(u32 row=col+1; row<P.n_rows; ++row) { out.at(row,col) = eT(0); }
00224         }
00225       }
00226     else
00227       {
00228       const eT* A_mem = A.memptr();
00229       
00230       for(u32 i=0; i<n_elem; ++i)
00231         {
00232         out_mem[i] = eop_core<eop_type>::process(x, A_mem[i]);
00233         }
00234       }
00235     }
00236   }
00237 
00238 
00239 
00240 template<typename eop_type>
00241 template<typename T1>
00242 arma_hot
00243 inline
00244 void
00245 eop_core<eop_type>::apply_inplace_plus(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00246   {
00247   arma_extra_debug_sigprint();
00248   
00249   typedef typename T1::elem_type eT;
00250   
00251   const Proxy<T1>& P = x.P;
00252   
00253   arma_debug_assert_same_size(out.n_rows, out.n_cols, P.n_rows, P.n_cols, "matrix addition");
00254   
00255         eT* out_mem = out.memptr();
00256   const u32 n_elem  = P.n_elem;
00257     
00258   if(is_generator<eop_type>::value == true)
00259     {
00260     for(u32 i=0; i<n_elem; ++i)
00261       {
00262       out_mem[i] += eop_aux::generate<eT,eop_type>();
00263       }
00264     }
00265   else
00266     {
00267     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00268       {
00269       for(u32 row=0; row<P.n_rows; ++row)
00270         {
00271         out.at(row,row) += eT(1);
00272         }
00273       }
00274     else
00275       {
00276       for(u32 i=0; i<n_elem; ++i)
00277         {
00278         out_mem[i] += eop_core<eop_type>::process(x, P[i]);
00279         }
00280       }
00281     }
00282   }
00283 
00284 
00285 
00286 template<typename eop_type>
00287 template<typename T1>
00288 arma_hot
00289 inline
00290 void
00291 eop_core<eop_type>::apply_inplace_minus(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00292   {
00293   arma_extra_debug_sigprint();
00294   
00295   typedef typename T1::elem_type eT;
00296   
00297   const Proxy<T1>& P = x.P;
00298   
00299   arma_debug_assert_same_size(out.n_rows, out.n_cols, P.n_rows, P.n_cols, "matrix subtraction");
00300   
00301         eT* out_mem = out.memptr();
00302   const u32 n_elem  = P.n_elem;
00303     
00304   if(is_generator<eop_type>::value == true)
00305     {
00306     for(u32 i=0; i<n_elem; ++i)
00307       {
00308       out_mem[i] -= eop_aux::generate<eT,eop_type>();
00309       }
00310     }
00311   else
00312     {
00313     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00314       {
00315       for(u32 row=0; row<P.n_rows; ++row)
00316         {
00317         out.at(row,row) -= eT(1);
00318         }
00319       }
00320     else
00321       {
00322       for(u32 i=0; i<n_elem; ++i)
00323         {
00324         out_mem[i] -= eop_core<eop_type>::process(x, P[i]);
00325         }
00326       }
00327     }
00328   }
00329 
00330 
00331 
00332 template<typename eop_type>
00333 template<typename T1>
00334 arma_hot
00335 inline
00336 void
00337 eop_core<eop_type>::apply_inplace_schur(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00338   {
00339   arma_extra_debug_sigprint();
00340   
00341   typedef typename T1::elem_type eT;
00342   
00343   const Proxy<T1>& P = x.P;
00344   
00345   arma_debug_assert_same_size(out.n_rows, out.n_cols, P.n_rows, P.n_cols, "element-wise matrix multiplication");
00346   
00347         eT* out_mem = out.memptr();
00348   const u32 n_elem  = P.n_elem;
00349   
00350   if(is_generator<eop_type>::value == true)
00351     {
00352     for(u32 i=0; i<n_elem; ++i)
00353       {
00354       out_mem[i] *= eop_aux::generate<eT,eop_type>();
00355       }
00356     }
00357   else
00358     {
00359     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00360       {
00361       for(u32 col=0; col<P.n_rows; ++col)
00362         {
00363         for(u32 row=0;     row<col;      ++row) { out.at(row,col) = eT(0); }
00364         for(u32 row=col+1; row<P.n_rows; ++row) { out.at(row,col) = eT(0); }
00365         }
00366       }
00367     else
00368       {
00369       for(u32 i=0; i<n_elem; ++i)
00370         {
00371         out_mem[i] *= eop_core<eop_type>::process(x, P[i]);
00372         }
00373       }
00374     }
00375   }
00376 
00377 
00378 
00379 template<typename eop_type>
00380 template<typename T1>
00381 arma_hot
00382 inline
00383 void
00384 eop_core<eop_type>::apply_inplace_div(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
00385   {
00386   arma_extra_debug_sigprint();
00387   
00388   typedef typename T1::elem_type eT;
00389   
00390   const Proxy<T1>& P = x.P;
00391   
00392   arma_debug_assert_same_size(out.n_rows, out.n_cols, P.n_rows, P.n_cols, "element-wise matrix division");
00393   
00394         eT* out_mem = out.memptr();
00395   const u32 n_elem  = P.n_elem;
00396   
00397   if(is_generator<eop_type>::value == true)
00398     {
00399     for(u32 i=0; i<n_elem; ++i)
00400       {
00401       out_mem[i] /= eop_aux::generate<eT,eop_type>();
00402       }
00403     }
00404   else
00405     {
00406     if(is_same_type<eop_type, eop_ones_diag>::value == true)
00407       {
00408       for(u32 col=0; col<P.n_rows; ++col)
00409         {
00410         for(u32 row=0;     row<col;      ++row) { out.at(row,col) /= eT(0); }
00411         for(u32 row=col+1; row<P.n_rows; ++row) { out.at(row,col) /= eT(0); }
00412         }
00413       }
00414     else
00415       {
00416       for(u32 i=0; i<n_elem; ++i)
00417         {
00418         out_mem[i] /= eop_core<eop_type>::process(x, P[i]);
00419         }
00420       }
00421     }
00422   }
00423 
00424 
00425 
00426 //! @}