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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_aligned.hpp
* DESCRIPTION: matrices aligned in memory
* COPYRIGHT  : (C) 2009  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_ALIGNED__HPP
#define __MMX__MATRIX_ALIGNED__HPP
#include <basix/int.hpp>
#include <numerix/simd.hpp>
#include <algebramix/vector_aligned.hpp>

namespace mmx {

/******************************************************************************
*  Variants for matrices aligned in memory
******************************************************************************/

template<typename V, typename W>
struct matrix_assume_aligned: public V {
  // All the columns of these matrices are assumed to be aligned
  // V is the variant to be used on base types
  // W is the variant to be used on simd types (see matrix_sse.hpp)
 
  typedef vector_assume_aligned<typename V::Vec,
                                typename W::Vec> Vec;
  typedef matrix_assume_aligned<typename V::Naive,
                                typename W::Naive> Naive;
  typedef matrix_assume_aligned<typename V::Positive,
                                typename W::Positive> Positive;
  typedef typename V::No_aligned No_aligned;
  typedef matrix_assume_aligned<typename V::No_thread,
                                typename W::No_thread> No_thread;
  typedef matrix_assume_aligned<typename V::No_scaled,
                                typename W::No_scaled> No_scaled;
};

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

template<typename Z, typename V, typename W>
struct implementation<vector_allocate,Z,matrix_assume_aligned<V,W> >:
  public implementation<vector_allocate,
                        typename matrix_assume_aligned<V,W>::Vec> {};

template<typename V, typename W>
struct matrix_aligned: public V {
  // These matrices are not necessarily aligned in memory
  typedef vector_aligned<typename V::Vec,
                         typename W::Vec> Vec;
  typedef matrix_aligned<typename V::Naive,
                         typename W::Naive> Naive;
  typedef matrix_aligned<typename V::Positive,
                         typename W::Positive> Positive;
  typedef typename V::No_aligned No_aligned;
  typedef matrix_aligned<typename V::No_thread,
                         typename W::No_thread> No_thread;
  typedef matrix_aligned<typename V::No_scaled,
                         typename W::No_scaled> No_scaled;
};

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

template<typename Z, typename V, typename W>
struct implementation<vector_allocate,Z,matrix_aligned<V,W> >:
  public implementation<vector_allocate,
                        typename matrix_aligned<V,W>::Vec> {};

/******************************************************************************
* Helpers to be specialized on aligned matrices hardware data
******************************************************************************/

template<typename V, typename W,
         typename Op, typename D, typename S1, typename S2>
struct mat_mul_aligned_helper {
  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 implementation<matrix_multiply,V> Mat;
    Mat::template mul<Op> (d, s1, s2, r, rr, l, ll, c, cc);
  }
};

/******************************************************************************
* Aligned matrix multiplication
******************************************************************************/

template<typename Z, typename V, typename W>
struct implementation<matrix_multiply_base,Z,matrix_assume_aligned<V,W> >:
  public implementation<matrix_linear,Z>
{
  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_mul_aligned_helper<V,W,Op,D,S1,S2>::
      mul (d, s1, s2, r, rr, l, ll, c, cc);
  }
}; // implementation<matrix_multiply_base,Z,matrix_assume_aligned<V,W> >

/******************************************************************************
* Matrix multiplication, not necessarily aligned
******************************************************************************/

template<typename Z, typename V, typename W>
struct implementation<matrix_multiply_base,Z,matrix_aligned<V,W> >:
  public implementation<matrix_linear,Z>
{
  typedef implementation<matrix_multiply,V> Mat;
  typedef implementation<matrix_multiply,matrix_assume_aligned<V,W> > AMat;
  typedef implementation<vector_allocate,
                         typename matrix_aligned<V,W>::Vec> Vec;

  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)
  {
    if (r == 0) return;
    if (mask_helper<D>::len != 16 || mask_helper<S1>::len != 16) {
      Mat::template mul<Op> (d, s1, s2, r, rr, l, ll, c, cc);
      return;
    }
    nat ar= Vec::template vec_floor_aligned_size<D> (r);
    nat sr= r - ar;
    nat arr= Vec::template vec_floor_aligned_size<D> (rr);    
    if (sr != 0) {
      if (ar != 0) mul<Op> (d, s1, s2, ar, rr, l, ll, c, cc);
      Mat::template mul<Op> (d  + Mat::index (ar, 0, rr, cc),
                             s1 + Mat::index (ar, 0, rr, ll), s2,
                             sr, rr, l, ll, c, cc);
      // NOTE: row times matrix can be further optimized for SIMD instructions
    }
    else if (! Vec::template vec_is_aligned (d) ||
             ! Vec::template vec_is_aligned (s1) ||
             arr != rr) {
      D* xd = mmx_new<D> (r * c);
      S1* xs1= mmx_new<S1> (r * l);
      S2* xs2= mmx_new<S2> (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);
      AMat::template mul<Op> (xd, xs1, xs2, r, r, l, l, c, c);
      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, r * c);
      mmx_delete<S1> (xs1, r * l);
      mmx_delete<S2> (xs2, l * c);
    }
    else
      AMat::template mul<Op> (d, s1, s2, r, rr, l, ll, c, c);
  }
  
}; // implementation<matrix_multiply_base,Z,matrix_aligned<V,W> >

} // namespace mmx
#endif // __MMX__MATRIX_ALIGNED__HPP