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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_threads.hpp
* DESCRIPTION: multi-threaded matrix-multiplication
* COPYRIGHT  : (C) 2007  Joris van der Hoeven
*******************************************************************************
* 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_THREADS__HPP
#define __MMX__MATRIX_THREADS__HPP
#include <algebramix/matrix_unrolled.hpp>
#include <basix/threads.hpp>

namespace mmx {

/******************************************************************************
* Multi-threaded matrix variant
******************************************************************************/

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

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

/******************************************************************************
* Multi-threaded matrix-multiplication
******************************************************************************/

#ifdef BASIX_ENABLE_THREADS

template<typename V, typename W>
struct implementation<matrix_multiply_base,V,matrix_threads<W> >:
  public implementation<matrix_linear,V>
{
  static const nat thr= 1024;
  static const nat sz = 4;
  typedef implementation<matrix_multiply,V,W> Mat;

template<typename Op, typename D, typename S1, typename S2>
struct multiply_task_rep: public task_rep {
  D* d; const S1* s1; const S2* s2;
  nat r; nat rr; nat l; nat ll; nat c; nat cc;
  //D *xd; S1 *xs1; S2 *xs2;
public:
  inline multiply_task_rep (D* d2, const S1* s1b, const S2* s2b,
                            nat r2, nat rr2, nat l2, nat ll2, nat c2, nat cc2):
    d (d2), s1 (s1b), s2 (s2b),
    r (r2), rr (rr2), l (l2), ll (ll2), c (c2), cc (cc2)
  {
    /*
    xd = mmx_new<D > (aligned_size<C,V> (r * c));
    xs1= mmx_new<S1> (aligned_size<C,V> (r * l));
    xs2= mmx_new<S2> (aligned_size<C,V> (l * c));
    Mat::template mat_unary_stride<id_op>
      (xs1, Mat::index (1, 0, r , l ), Mat::index (0, 1, r , l ),
       s1 , Mat::index (1, 0, rr, ll), Mat::index (0, 1, rr, ll), r, l);
    Mat::template mat_unary_stride<id_op>
      (xs2, Mat::index (1, 0, l , c ), Mat::index (0, 1, l , c ),
       s2 , Mat::index (1, 0, ll, cc), Mat::index (0, 1, ll, cc), l, c);
    */
  }
  inline ~multiply_task_rep () {
    /*
    Mat::template mat_unary_stride<typename Op::nomul_op>
      (d , Mat::index (1, 0, rr, cc), Mat::index (0, 1, rr, cc),
       xd, Mat::index (1, 0, r , c ), Mat::index (0, 1, r , c ), r, c);
    mmx_delete<D > (xd , aligned_size<C,V> (r * c));
    mmx_delete<S1> (xs1, aligned_size<C,V> (r * l));
    mmx_delete<S2> (xs2, aligned_size<C,V> (l * c));
    */
  }
  void execute () {
    //Mat::template mul<mul_op> (xd, xs1, xs2, r, r, l, l, c, c);
    Mat::template mul<mul_op> (d, s1, s2, r, rr, l, ll, c, cc);
  }
};

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)
{
  typedef typename Op::acc_op Acc;
  if (r * c < thr) Mat::template mul<Op> (d, s1, s2, r, rr, l, ll, c, cc);
  else {
    nat tr= r, nr= 1, tc= c, nc= 1, tt= threads_number;
    while (tt != 1) {
      if ((tt & 1) == 0) {
        if (tr > tc) { tr= (tr+1) >> 1; nr <<= 1; }
        else { tc= (tc+1) >> 1; nc <<= 1; }
        tt >>= 1;
      }
      else {
        if (tr > tc) { tr= (tr+tt-1) / tt; nr *= tt; }
        else { tc= (tc+tt-1) / tt; nc *= tt; }
        tt= 1;
      }
    }
    tr= sz * ((tr + sz - 1) / sz);
    tc= sz * ((tc + sz - 1) / sz);

    task tasks[nr*nc];
    for (nat ir=0; ir<nr; ir++)
      for (nat ic=0; ic<nc; ic++) {
        nat r1= ir * tr, r2= min (r1 + tr, r);
        nat c1= ic * tc, c2= min (c1 + tc, c);
        if (r1 < r && c1 < c) {
          D*        td = d  + Mat::index (r1, c1, rr, cc);
          const S1* ts1= s1 + Mat::index (r1, 0 , rr, ll);
          const S2* ts2= s2 + Mat::index (0 , c1, ll, cc);
          tasks[Mat::index (ir, ic, nr, nc)]=
            new multiply_task_rep<Op,D,S1,S2>
                  (td, ts1, ts2, r2-r1, rr, l, ll, c2-c1, cc);
        }
      }
    threads_execute (tasks, nr*nc);
  }
}

}; // implementation<matrix_multiply_base,V,matrix_threads<W> >

#endif // BASIX_ENABLE_THREADS

} // namespace mmx
#endif //__MMX__MATRIX_THREADS__HPP