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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_unrolled.hpp
* DESCRIPTION: unrolling and subdividing 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_UNROLLED__HPP
#define __MMX__MATRIX_UNROLLED__HPP
#include <algebramix/vector_unrolled.hpp>
#include <algebramix/matrix_fixed.hpp>

namespace mmx {

/******************************************************************************
* Unrolled variant
******************************************************************************/

template<nat sz, typename V=matrix_naive>
struct matrix_unrolled: public V {
  typedef vector_unrolled<sz,typename V::Vec> Vec;
  typedef typename V::Naive Naive;
  typedef matrix_unrolled<sz,typename V::Positive> Positive;
  typedef matrix_unrolled<sz,typename V::No_aligned> No_aligned;
  typedef matrix_unrolled<sz,typename V::No_simd> No_simd;
  typedef matrix_unrolled<sz,typename V::No_thread> No_thread;
  typedef matrix_unrolled<sz,typename V::No_scaled> No_scaled;
};

template<nat sz, typename F, typename V, typename W>
struct implementation<F,V,matrix_unrolled<sz,W> >:
  public implementation<F,V,W> {};

/******************************************************************************
* Complete a few missing rows, columns, inner products using another algorithm
******************************************************************************/

template<typename Op, typename V, typename D, typename S1, typename S2> void
mul_complete (D* dest, const S1* src1, const S2* src2,
              nat r, nat rr, nat l, nat ll, nat c, nat cc,
              nat hr, nat hl, nat hc)
{
  typedef implementation<matrix_multiply,V> Mat;
  typedef typename Op::acc_op Acc;
  if (hr < r && hl != 0 && hc != 0)
    Mat::template mul<Op > (dest + Mat::index (hr, 0, rr, cc),
                            src1 + Mat::index (hr, 0, rr, ll),
                            src2,
                            r-hr, rr, hl, ll, hc, cc);
  if (hc < c && hl != 0)
    Mat::template mul<Op > (dest + Mat::index (0, hc, rr, cc),
                            src1,
                            src2 + Mat::index (0, hc, ll, cc),
                            r , rr, hl, ll, c-hc, cc);
  if (hl < l)
    Mat::template mul<Acc> (dest,
                            src1 + Mat::index (0, hl, rr, ll),
                            src2 + Mat::index (hl, 0, ll, cc),
                            r , rr, l-hl, ll, c , cc);
}

/******************************************************************************
* Using a triple loop using big unrolled blocks
******************************************************************************/

template<typename Op, nat ur, nat ul, nat uc, typename V,
         typename D, typename S1, typename S2> void
mul_unrolled (D* dest, const S1* src1, const S2* src2,
              nat r, nat rr, nat l, nat ll, nat c, nat cc)
{
  typedef implementation<matrix_multiply,V> Mat;
  typedef implementation<matrix_multiply_base,matrix_naive> NMat;
  typedef typename Op::acc_op Acc;
  nat nr= r/ur, nl= l/ul, nc= c/uc;
  if (nl == 0)
    NMat::template clr<Op> (dest, r, rr, c, cc);
  else
    for (nat ir=0; ir<nr; ir++)
      for (nat ic=0; ic<nc; ic++) {
        nat il=0;
        for (; il<1; il++)
          matrix_multiply_helper<Op,D,S1,S2,ur,ul,uc>::
            mul_stride (dest + Mat::index (ir*ur, ic*uc, rr, cc),
                        src1 + Mat::index (ir*ur, il*ul, rr, ll),
                        src2 + Mat::index (il*ul, ic*uc, ll, cc),
                        rr, ll);
        for (; il<nl; il++)
          matrix_multiply_helper<Acc,D,S1,S2,ur,ul,uc>::
            mul_stride (dest + Mat::index (ir*ur, ic*uc, rr, cc),
                        src1 + Mat::index (ir*ur, il*ul, rr, ll),
                        src2 + Mat::index (il*ul, ic*uc, ll, cc),
                        rr, ll);
      }
  mul_complete<Op,V> (dest, src1, src2, r, rr, l, ll, c, cc,
                      ur*nr, ul*nl, uc*nc);
}

/******************************************************************************
* Interface for unrolled multiplication
******************************************************************************/

template<nat sz, typename V, typename W>
struct implementation<matrix_multiply_base,V,matrix_unrolled<sz,W> >:
  public implementation<matrix_linear,V>
{
  const static nat ur= sz;
  const static nat ul= sz;
  const static nat uc= sz;

  template<typename Op, typename D, typename S1, typename S2>
  static inline void
  mul (D* dest, const S1* src1, const S2* src2,
       nat r, nat rr, nat l, nat ll, nat c, nat cc)
  {
    mul_unrolled<Op,ur,ul,uc,W> (dest, src1, src2, r, rr, l, ll, c, cc);
  }
}; // implementation<matrix_multiply,V,matrix_unrolled<sz,BV> >

} // namespace mmx
#endif //__MMX__MATRIX_UNROLLED__HPP