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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : structured_block.hpp
* DESCRIPTION: block algorithms for structured matrices
* COPYRIGHT  : (C) 2009  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__STRUCTURED_BLOCK__HPP
#define __MMX__STRUCTURED_BLOCK__HPP
#include <algebramix/structured_naive.hpp>
#include <algebramix/matrix_block.hpp>

namespace mmx {

/******************************************************************************
* Multiplication of structured matrices
******************************************************************************/

struct matrix_structured_multiply {};

template<typename V>
struct implementation<matrix_structured_multiply,V,matrix_naive>:
  public implementation<matrix_multiply,V>
{
  typedef implementation<matrix_multiply,V> Mat;
  typedef implementation<vector_linear,V> Vec;

static inline
half_size (nat n) {
  return (n+1) >> 1;
}

template<typename Op,
         typename TD, typename T1, typename T2,
         typename D , typename S1, typename S2> static void
smul (D* d, const S1* s1, const S2* s2,
      nat r, nat rr, nat l, nat ll, nat c, nat cc);

template<typename TD, typename T1, typename T2>
struct structured_helper {
  template<typename Op, typename D, typename S1, typename S2> static 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 (l == 0) Mat::clr<Op> (d, r, rr, c, cc);
    else if (l == 1) Op::set_op (*d, *s1, *s2);
    else {
      nat hr= half_size (r);
      nat hl= half_size (l);
      nat hc= half_size (c);

      smul<Op, typename TD::ul, typename T1::ul, typename T2::ul>
        (d,
         s1,
         s2,
         hr, rr, hl, ll, hc, cc);
      smul<Acc, typename TD::ul, typename T1::ur, typename T2::ll>
        (d,
         s1 + Mat::index (0, hl, rr, ll),
         s2 + Mat::index (hl, 0, ll, cc),
         hr, rr, ll - hl, ll, hc, cc);
      smul<Op, typename TD::ur, typename T1::ul, typename T2::ur>
        (d + Mat::index (0, hc, rr, cc),
         s1,
         s2 + Mat::index (0, hc, ll, cc),
         hr, rr, hl, ll, cc - hc, cc);
      smul<Acc, typename TD::ur, typename T1::ur, typename T2::lr>
        (d + Mat::index (0, hc, rr, cc),
         s1 + Mat::index (0, hl, rr, ll),
         s2 + Mat::index (hl, hc, ll, cc),
         hr, rr, ll - hl, ll, cc - hc, cc);

      smul<Op, typename TD::ll, typename T1::lr, typename T2::ll>
        (d + Mat::index (hr, 0, rr, cc),
         s1 + Mat::index (hr, hl, rr, ll),
         s2 + Mat::index (hl, 0, ll, cc),
         rr - hr, rr, ll - hl, ll, hc, cc);
      smul<Acc, typename TD::ll, typename T1::ll, typename T2::ul>
        (d + Mat::index (hr, 0, rr, cc),
         s1 + Mat::index (hr, 0, rr, ll),
         s2,
         rr - hr, rr, hl, ll, hc, cc);
      smul<Op, typename TD::lr, typename T1::lr, typename T2::lr>
        (d + Mat::index (hr, hc, rr, cc),
         s1 + Mat::index (hr, hl, rr, ll),
         s2 + Mat::index (hl, hc, ll, cc),
         rr - hr, rr, ll - hl, ll, cc - hc, cc);
      smul<Acc, typename TD::lr, typename T1::ll, typename T2::ur>
        (d + Mat::index (hr, hc, rr, cc),
         s1 + Mat::index (hr, 0, rr, ll),
         s2 + Mat::index (0, hc, ll, cc),
         rr - hr, rr, hl, ll, cc - hc, cc);
    }
  }
};

struct noop_helper {
  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) {
    (void) d; (void) s1; (void) s2;
    (void) r; (void) rr; (void) l; (void) ll; (void) c; (void) cc;
  }
};

template<typename T1, typename T2>
struct structured_helper<zero_matrix,T1,T2>: public noop_helper {};

template<typename TD, typename T2>
struct structured_helper<TD,zero_matrix,T2>: public noop_helper {};

template<typename TD, typename T1>
struct structured_helper<TD,T1,zero_matrix>: public noop_helper {};

template<typename TD>
struct structured_helper<TD,zero_matrix,zero_matrix>: public noop_helper {};

template<typename T1>
struct structured_helper<zero_matrix,T1,zero_matrix>: public noop_helper {};

template<typename T2>
struct structured_helper<zero_matrix,zero_matrix,T2>: public noop_helper {};

template<>
struct structured_helper<zero_matrix,zero_matrix,zero_matrix>:
  public noop_helper {};

template<typename Op,
         typename TD, typename T1, typename T2,
         typename D, typename S1, typename S2> static inline void
smul (D* d, const S1* s1, const S2* s2,
      nat r, nat rr, nat l, nat ll, nat c, nat cc) {
  structured_helper<TD,T1,T2>::template mul<Op>
    (d, s1, s2, r, rr, l, ll, c, cc);
}

}; // implementation<matrix_structured,V,matrix_naive>

} // namespace mmx
#endif //__MMX__STRUCTURED_BLOCK__HPP