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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : structured_naive.hpp
* DESCRIPTION: naive 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_NAIVE__HPP
#define __MMX__STRUCTURED_NAIVE__HPP
#include <algebramix/matrix_naive.hpp>

namespace mmx {

/******************************************************************************
* Supported matrix types
******************************************************************************/

struct full_matrix {
  typedef full_matrix ul;
  typedef full_matrix ur;
  typedef full_matrix ll;
  typedef full_matrix lr;
  static inline nat first_row (nat c, nat rs) {
    (void) c; (void) rs; return 0; }
  static inline nat first_col (nat r, nat cs) {
    (void) r; (void) cs; return 0; }
  static inline nat last_row  (nat c, nat rs) {
    (void) c; return rs; }
  static inline nat last_col  (nat r, nat cs) {
    (void) r; return cs; }
};

struct zero_matrix {
  typedef zero_matrix ul;
  typedef zero_matrix ur;
  typedef zero_matrix ll;
  typedef zero_matrix lr;
  static inline nat first_row (nat c, nat rs) {
    (void) c; return rs; }
  static inline nat first_col (nat r, nat cs) {
    (void) r; return cs; }
  static inline nat last_row  (nat c, nat rs) {
    (void) c; return rs; }
  static inline nat last_col  (nat r, nat cs) {
    (void) r; return cs; }
};

struct diagonal_matrix {
  typedef diagonal_matrix ul;
  typedef zero_matrix ur;
  typedef zero_matrix ll;
  typedef diagonal_matrix lr;
  static inline nat first_row (nat c, nat rs) {
    (void) rs; return c; }
  static inline nat first_col (nat r, nat cs) {
    (void) cs; return r; }
  static inline nat last_row  (nat c, nat rs) {
    (void) rs; return c+1; }
  static inline nat last_col  (nat r, nat cs) {
    (void) cs; return r+1; }
};

struct upper_triangular_matrix {
  typedef upper_triangular_matrix ul;
  typedef full_matrix ur;
  typedef zero_matrix ll;
  typedef upper_triangular_matrix lr;
  static inline nat first_row (nat c, nat rs) {
    (void) c; (void) rs; return 0; }
  static inline nat first_col (nat r, nat cs) {
    (void) cs; return r; }
  static inline nat last_row  (nat c, nat rs) {
    (void) rs; return c+1; }
  static inline nat last_col  (nat r, nat cs) {
    (void) r; return cs; }
};

struct lower_triangular {
  typedef lower_triangular_matrix ul;
  typedef zero_matrix ur;
  typedef full_matrix ll;
  typedef lower_triangular_matrix lr;
  static inline nat first_row (nat c, nat rs) {
    (void) rs; return c; }
  static inline nat first_col (nat r, nat cs) {
    (void) r; (void) cs; return 0; }
  static inline nat last_row  (nat c, nat rs) {
    (void) c; return rs; }
  static inline nat last_col  (nat r, nat cs) {
    (void) cs; return r+1; }
};

/******************************************************************************
* 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;

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 rs, nat rr, nat ls, nat ll, nat cs, nat cc) {
  typedef typename Op::acc_op Acc;
  typedef typename Op::clr_op Clr;
  const S1 *s1r, *s1c;
  const S2 *s2r, *s2c;
  nat drp = Mat::index (1, 0, rr, cc), dcp = Mat::index (0, 1, rr, cc);
  nat s1p= Mat::index (0, 1, rr, ll);
  nat s2p= Mat::index (1, 0, ll, cc);
  for (nat ri=0; ri<rs; ri++) {
    nat c1= TD::first_col (ri, cs);
    nat c2= TD::last_col  (ri, cs);
    for (nat ci=0; ci<cs; ci++) {
      nat l1= max (T1::first_col (ri, ls), T2::first_row (ci, ls));
      nat l2= min (T1::last_col  (ri, ls), T2::last_row  (ci, ls));
      D tmp= dest [Mat::index (ri, ci, rr, cc)];
      Clr::set_op (tmp);
      S1* s1c= s1 + Mat::index (ri, 0, rr, ll);
      S1* s2c= s2 + Mat::index (0, ci, ll, cc);
      for (nat li=0; li<ls; li++, s1c += s1p, s2c += s2p)
        Acc::set_op (tmp, *s1c, *s2c);
      dest [Mat::index (ri, ci, rr, cc)]= tmp;
    }
  }
}

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

} // namespace mmx
#endif //__MMX__STRUCTURED_NAIVE__HPP