include/algebramix/matrix_strassen.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_strassen.hpp
* DESCRIPTION: Matrix multiplication using Strassen's algorithm
* 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_STRASSEN__HPP
#define __MMX__MATRIX_STRASSEN__HPP
#include <algebramix/matrix_unrolled.hpp>

namespace mmx {

/******************************************************************************
* Variant for Strassen multiplication
******************************************************************************/

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

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

/******************************************************************************
* Strassen multiplication (even-odd decomposition for numerical stability)
******************************************************************************/

template<typename V, typename W>
struct implementation<matrix_multiply_base,V,matrix_strassen<W> >:
  public implementation<matrix_linear,V>
{
  static const nat thr= 128; // use 512 for double >> TODO
  typedef implementation<vector_linear,W> Vec;
  typedef implementation<matrix_multiply,W> Mat;

  template<typename C> static inline void
  mat_load (C* d, nat r, nat c, const C* s, nat rr, nat cc) {
    nat j= c;
    for (; j!=0; j--, d += Mat::index(0,1,r,c), s += Mat::index(0,2,rr,cc)) {
      nat      i = r;
      C*       dd= d;
      const C* ss= s;
      for (; i!=0; i--, dd += Mat::index(1,0,r,c), ss += Mat::index(2,0,rr,cc))
        *dd= *ss;
    } 
  }

  template<typename Op, typename C> static inline void
  mat_save (C* d, nat rr, nat cc, const C* s, nat r, nat c) {
    typedef typename Op::nomul_op Set;
    nat j= c;
    for (; j!=0; j--, d += Mat::index(0,2,rr,cc), s += Mat::index(0,1,r,c)) {
      nat      i = r;
      C*       dd= d;
      const C* ss= s;
      for (; i!=0; i--, dd += Mat::index(2,0,rr,cc), ss += Mat::index(1,0,r,c))
        Set::set_op (*dd, *ss);
    } 
  }

  template<typename Op, typename D, typename S1, typename S2> static void
  mul (D* d, const S1* a, const S2* b,
       nat r, nat rr, nat l, nat ll, nat c, nat cc)
  {
    if (r < thr || l < thr || c < thr) {
      Mat::template mul<Op> (d, a, b, r, rr, l, ll, c, cc);
      return;
    }

    nat hr= r>>1, hl= l>>1, hc= c>>1, fr= hr<<1, fl= hl<<1, fc= hc<<1;
    
    nat sza= aligned_size<S1,W> (hr * hl);
    S1* a11= mmx_new<S1> (5 * sza);
    S1* a12= a11 + sza;
    S1* a21= a12 + sza;
    S1* a22= a21 + sza;
    S1* aaa= a22 + sza;

    nat szb= aligned_size<S2,W> (hl * hc);
    S2* b11= mmx_new<S2> (5 * szb);
    S2* b12= b11 + szb;
    S2* b21= b12 + szb;
    S2* b22= b21 + szb;
    S2* bbb= b22 + szb;

    nat szd= aligned_size<D,W> (hr * hc);
    D*  m1 = mmx_new<D> (11 * szd);
    D*  m2 = m1  + szd;
    D*  m3 = m2  + szd;
    D*  m4 = m3  + szd;
    D*  m5 = m4  + szd;
    D*  m6 = m5  + szd;
    D*  m7 = m6  + szd;
    D*  d11= m7  + szd;
    D*  d12= d11 + szd;
    D*  d21= d12 + szd;
    D*  d22= d21 + szd;

    mat_load (a11, hr, hl, a                            , rr, ll);
    mat_load (a12, hr, hl, a + Mat::index (0, 1, rr, ll), rr, ll);
    mat_load (a21, hr, hl, a + Mat::index (1, 0, rr, ll), rr, ll);
    mat_load (a22, hr, hl, a + Mat::index (1, 1, rr, ll), rr, ll);
    mat_load (b11, hr, hl, b                            , rr, ll);
    mat_load (b12, hr, hl, b + Mat::index (0, 1, rr, ll), rr, ll);
    mat_load (b21, hr, hl, b + Mat::index (1, 0, rr, ll), rr, ll);
    mat_load (b22, hr, hl, b + Mat::index (1, 1, rr, ll), rr, ll);

    Vec::add (aaa, a11, a22, hr * hl);
    Vec::add (bbb, b11, b22, hl * hc);
    mul<mul_op> (m1, aaa, bbb, hr, hr, hl, hl, hc, hc);
    Vec::add (aaa, a21, a22, hr * hl);
    mul<mul_op> (m2, aaa, b11, hr, hr, hl, hl, hc, hc);
    Vec::sub (bbb, b12, b22, hl * hc);
    mul<mul_op> (m3, a11, bbb, hr, hr, hl, hl, hc, hc);
    Vec::sub (bbb, b21, b11, hl * hc);
    mul<mul_op> (m4, a22, bbb, hr, hr, hl, hl, hc, hc);
    Vec::add (aaa, a11, a12, hr * hl);
    mul<mul_op> (m5, aaa, b22, hr, hr, hl, hl, hc, hc);
    Vec::sub (aaa, a21, a11, hr * hl);
    Vec::add (bbb, b11, b12, hl * hc);
    mul<mul_op> (m6, aaa, bbb, hr, hr, hl, hl, hc, hc);
    Vec::sub (aaa, a12, a22, hr * hl);
    Vec::add (bbb, b21, b22, hl * hc);
    mul<mul_op> (m7, aaa, bbb, hr, hr, hl, hl, hc, hc);

    Vec::add (d11, m1, m4, hr * hc);
    Vec::sub (d11, m5, hr * hc);
    Vec::add (d11, m7, hr * hc);
    Vec::add (d12, m3, m5, hr * hc);
    Vec::add (d21, m2, m4, hr * hc);
    Vec::sub (d22, m1, m2, hr * hc);
    Vec::add (d22, m3, hr * hc);
    Vec::add (d22, m6, hr * hc);

    mat_save<Op> (d                            , rr, cc, d11, hr, hc);
    mat_save<Op> (d + Mat::index (0, 1, rr, cc), rr, ll, d12, hr, hc);
    mat_save<Op> (d + Mat::index (1, 0, rr, cc), rr, ll, d21, hr, hc);
    mat_save<Op> (d + Mat::index (1, 1, rr, cc), rr, ll, d22, hr, hc);

    mmx_delete<S1> (a11, 5  * sza);
    mmx_delete<S2> (b11, 5  * szb);
    mmx_delete<D > (m1 , 11 * szd);

    mul_complete<Op,W> (d, a, b, r, rr, l, ll, c, cc, fr, fl, fc);
  }
}; // implementation<matrix_multiply_base,V,matrix_strassen<W> >

} // namespace mmx
#endif //__MMX__MATRIX_STRASSEN__HPP

---