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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_fixed.hpp
* DESCRIPTION: inlined low-level matrix operations of a fixed size
* COPYRIGHT  : (C) 2007  Joris van der Hoeven and 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_FIXED__HPP
#define __MMX__MATRIX_FIXED__HPP
#include <algebramix/matrix_naive.hpp>

namespace mmx {

/******************************************************************************
* Inlined matrix multiplication for fixed sizes
******************************************************************************/

template<typename Op, nat tp> struct div_type {};

template<typename Op, typename D, typename S1, typename S2,
         nat r, nat l, nat c>
struct matrix_multiply_helper {
  static const nat tp= (r*c==1 || (r<l && c<l && l>8))? 3: (r<=c? 2: 1);

  template<typename rr, typename cc> static inline void
  mul (D* d, const S1* m1, const S2* m2) {
    matrix_multiply_helper <div_type<Op,tp>, D, S1, S2, r, l, c>::
      template mul<rr,cc> (d, m1, m2);
  }

  static inline void
  mul_stride (D* d, const S1* m1, const S2* m2, nat rr, nat ll) {
    matrix_multiply_helper <div_type<Op,tp>, D, S1, S2, r, l, c>::
      mul_stride (d, m1, m2, rr, ll);
  }
};

template<typename Op, typename D, typename S1, typename S2,
         nat r, nat l, nat c>
struct matrix_multiply_helper<div_type<Op,1>, D, S1, S2, r, l, c> {
  // Recurse on a smaller number of rows r
  static const nat r1= (r>>1), r2= r-r1;

  template<nat rr, nat ll> static inline void
  mul (D* d, const S1* m1, const S2* m2) {
    {
      matrix_multiply_helper <Op, D, S1, S2, r1, l, c>::
        template mul<rr,ll> (d, m1, m2);
      matrix_multiply_helper <Op, D, S1, S2, r2, l, c>::
        template mul<rr,ll> (d + r1, m1 + r1, m2);
    }
  }

  static inline void
  mul_stride (D* d, const S1* m1, const S2* m2, nat rr, nat ll) {
    {
      matrix_multiply_helper <Op, D, S1, S2, r1, l, c>::
	mul_stride (d, m1, m2, rr, ll);
      matrix_multiply_helper <Op, D, S1, S2, r2, l, c>::
	mul_stride (d + r1, m1 + r1, m2, rr, ll);
    }
  }
};

template<typename Op, typename D, typename S1, typename S2,
         nat r, nat l, nat c>
struct matrix_multiply_helper<div_type<Op,2>, D, S1, S2, r, l, c> {
  // Recurse on a smaller number of columns c
  static const nat c1= (c>>1), c2= c-c1;

  template<nat rr, nat ll> static inline void
  mul (D* d, const S1* m1, const S2* m2) {
    {
      matrix_multiply_helper <Op, D, S1, S2, r, l, c1>::
        template mul<rr,ll> (d, m1, m2);
      matrix_multiply_helper <Op, D, S1, S2, r, l, c2>::
        template mul<rr,ll> (d + c1 * rr, m1, m2 + c1 * ll);
    }
  }

  static inline void
  mul_stride (D* d, const S1* m1, const S2* m2, nat rr, nat ll) {
    {
      matrix_multiply_helper <Op, D, S1, S2, r, l, c1>::
	mul_stride (d, m1, m2, rr, ll);
      matrix_multiply_helper <Op, D, S1, S2, r, l, c2>::
	mul_stride (d + c1 * rr, m1, m2 + c1 * ll, rr, ll);
    }
  }
};

template<typename Op, typename D, typename S1, typename S2,
         nat r, nat l, nat c>
struct matrix_multiply_helper<div_type<Op,3>, D, S1, S2, r, l, c> {
  // Recurse on a smaller length l of the inner loop
  static const nat l1= (l>>1), l2= l-l1;
  typedef typename Op::acc_op Acc;

  template<nat rr, nat ll> static inline void
  mul (D* d, const S1* m1, const S2* m2) {
    {
      matrix_multiply_helper <Op , D, S1, S2, r, l1, c>::
        template mul<rr,ll> (d, m1, m2);
      matrix_multiply_helper <Acc, D, S1, S2, r, l2, c>::
        template mul<rr,ll> (d, m1 + l1 * rr, m2 + l1);
    }
  }

  static inline void
  mul_stride (D* d, const S1* m1, const S2* m2, nat rr, nat ll) {
    {
      matrix_multiply_helper <Op , D, S1, S2, r, l1, c>::
	mul_stride (d, m1, m2, rr, ll);
      matrix_multiply_helper <Acc, D, S1, S2, r, l2, c>::
	mul_stride (d, m1 + l1 * rr, m2 + l1, rr, ll);
    }
  }
};

template<typename Op, typename D, typename S1, typename S2>
struct matrix_multiply_helper<Op, D, S1, S2, 1, 1, 1> {
  template<nat rr, nat ll> static inline void
  mul (D* d, const S1* m1, const S2* m2) {
    Op::set_op (*d, *m1, *m2);
  }

  static inline void
  mul_stride (D* d, const S1* m1, const S2* m2, nat rr, nat ll) {
    (void) rr; (void) ll;
    Op::set_op (*d, *m1, *m2);
  }
};

} // namespace mmx
#endif //__MMX__MATRIX_FIXED__HPP