/Users/mourrain/Devel/mmx/algebramix/include/algebramix/matrix_balanced.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_balanced.hpp
* DESCRIPTION: degree balancing
* COPYRIGHT  : (C) 2008  Gregoire Lecerf
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX__MATRIX_BALANCED__HPP
#define __MMX__MATRIX_BALANCED__HPP
#include <basix/operators.hpp>
#include <algebramix/matrix_naive.hpp>

namespace mmx {

/******************************************************************************
* Balanced algorithms for matrixs
******************************************************************************/

template<typename V>
struct matrix_balanced: public V {
  typedef typename V::Vec Vec;
  typedef typename V::Naive Naive;
  typedef matrix_balanced<typename V::Positive> Positive;
  typedef matrix_balanced<typename V::No_simd> No_simd;
  typedef matrix_balanced<typename V::No_thread> No_thread;
  typedef matrix_balanced<typename V::No_scaled> No_scaled;
};

template<typename F, typename V, typename W>
struct implementation<F,V,matrix_balanced<W> >:
  public implementation<F,V,W> {};

/******************************************************************************
* Threshold
******************************************************************************/

template<typename V>
struct matrix_balanced_threshold {};

template<typename W>
struct threshold_helper<integer,matrix_balanced_threshold<W> > {
  typedef fixed_value<nat,2> impl;
};

/******************************************************************************
* Multiplication
******************************************************************************/

template<typename V, typename W>
struct implementation<matrix_multiply,V,matrix_balanced<W> >:
  public implementation<matrix_multiply_base,V>
{
  typedef matrix_multiply_threshold<matrix_balanced<W> > Th;
  typedef matrix_balanced_threshold<matrix_balanced<W> > BTh;
  typedef implementation<matrix_multiply,W> Mat;
  typedef implementation<vector_linear, typename W::Vec> Vec;

template<typename C>
static inline nat mat_size (const C* s, nat r, nat c) {
  nat sz= 0;
  for (nat i= 0; i < r * c; i++) sz= max (sz, N(s[i]));
  return sz; }

template<typename C>
static inline void mat_rshift (C* d, const C* s, nat r, nat c, int h) {
  for (nat i= 0; i < r * c; i++) d[i]= rshiftz (s[i], h); }

template<typename C>
static inline void mat_co_rshift (C* d, const C* s, nat r, nat c, int h) {
  for (nat i= 0; i < r * c; i++) d[i]= co_rshiftz (s[i], h); }

template<typename C>
static inline void mat_lshift (C* d, const C* s, nat r, nat c, int h) {
  for (nat i= 0; i < r * c; i++) d[i]= lshiftz (s[i], h); }

template<typename Op, typename D, typename S1, typename S2>
static inline void
mul (D* d, const S1* s1, const S2* s2,
     nat r, nat rr, nat l, nat ll, nat c, nat cc) {
  Mat::template mul<Op> (d, s1, s2, r, rr, l, ll, c, cc); }

template<typename D, typename S1, typename S2>
static inline void
mul (D* d, const S1* s1, const S2* s2,
     nat r, nat l, nat c) {
  Mat::template mul<mul_op> (d, s1, s2, r, r, l, l, c, c); }

template<typename C>
static void
mul (C* d, const C* s1, const C* s2,
     nat r, nat l, nat c) {
  const nat th= Threshold(C,Th);
  const nat balance_threshold= Threshold(C,BTh);
  nat sz1= mat_size (s1, r, l);
  nat sz2= mat_size (s2, l, c);
  if (sz1 == 0 || sz2 == 0) 
    Mat::clear (d, r, c);
  else if (sz1 <= th && sz2 <= th) 
    Mat::template mul (d, s1, s2, r, l, c);
  else if (balance_threshold * sz1 < sz2) {
    int h= (int) sz2 >> 1;
    nat len_aux= aligned_size<C,V> (l, c);
    C* aux= mmx_new<C> (len_aux);
    nat len_tmp= aligned_size<C,V> (r, c);
    C* tmp= mmx_new<C> (len_tmp);
    mat_rshift (aux, s2, l, c, h);
    mul (tmp, s1, aux, r, l, c);
    mat_lshift (d, tmp, r, c, h);
    mat_co_rshift (aux, s2, l, c, h);
    mul (tmp, s1, aux, r, l, c);
    Vec::add (d, tmp, r * c);
    mmx_delete<C> (aux, len_aux);
    mmx_delete<C> (tmp, len_tmp);
  }
  else if (balance_threshold * sz2 < sz1) {
    nat h= (int) sz1 >> 1;
    nat len_aux= aligned_size<C,V> (r, l);
    C* aux= mmx_new<C> (len_aux);
    nat len_tmp= aligned_size<C,V> (r, c);
    C* tmp= mmx_new<C> (len_tmp);
    mat_rshift (aux, s1, r, l, h);
    mul (tmp, aux, s2, r, l, c);
    mat_lshift (d, tmp, r, c, h);
    mat_co_rshift (aux, s1, r, l, h);
    mul (tmp, aux, s2, r, l, c);
    Vec::add (d, tmp, r * c);
    mmx_delete<C> (aux, len_aux);
    mmx_delete<C> (tmp, len_tmp);
  }
  else {
    Mat::template mul (d, s1, s2, r, l, c); } }

}; // implementation<matrix_multiply,V,matrix_balanced<W,b> >

} // namespace mmx
#endif //__MMX__MATRIX_BALANCED__HPP