include/algebramix/matrix_tft.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_tft.hpp
* DESCRIPTION: Tft based multiplication of matrices over polynomials
* COPYRIGHT  : (C) 2012  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_TFT_HPP
#define __MMX_MATRIX_TFT_HPP
#include <algebramix/matrix.hpp>
#include <algebramix/fft_truncated.hpp>
namespace mmx {

/******************************************************************************
* Variant for matrices with TFT based operations
******************************************************************************/

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

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

/******************************************************************************
* Size helper (for polynomials only)
******************************************************************************/

template<typename C>
struct matrix_tft_multiply_helper {

  typedef Scalar_type (C) F; // Fft field
  typedef fft_truncated_transformer<F> tft_transformer;
  static nat size (const C* s1, nat s1_rs, nat s1_cs,
                   const C* s2, nat s2_rs, nat s2_cs,
                   nat r, nat l, nat c) {
    nat sz= 0;
    for (nat k= 0; k < l; k++) {
      nat sz1= 0, sz2= 0;
      const C* ss1= s1 + k * s1_cs;
      const C* ss2= s2 + k * s2_rs;
      for (nat i= 0; i < r; i++, ss1 += s1_rs) sz1= max (sz1, N (*ss1));
      for (nat j= 0; j < c; j++, ss2 += s2_cs) sz2= max (sz2, N (*ss2));
      sz= max (sz, (sz1 == 0 || sz2 == 0) ? 0 : (sz1 + sz2 - 1));
    }
    return sz; }
};

/******************************************************************************
* Matrix product variant
******************************************************************************/

template<typename V, typename W>
struct implementation<matrix_multiply_base,V,matrix_tft<W> >:
  public implementation<matrix_linear,V>
{
  typedef implementation<matrix_multiply,W> Mat;

template<typename C, typename F, typename MV, typename Tfter>
static void
mat_direct_tft (matrix<F,MV>* dest, const C* s,
                nat s_rs, nat s_cs, nat r, nat c, Tfter& tfter) {
  nat n= tfter.len;
  nat m= next_power_of_two (n);
  for (nat k= 0; k < n; k++)
    dest[k]= matrix<F,MV> (F (), r, c);
  F* aux= mmx_new<F> (m);
  for (nat i= 0; i < r; i++)
    for (nat j= 0; j < c; j++) {
      for (nat k= 0; k < n ; k++)
        aux[k]= s[i * s_rs + j * s_cs] [k];
      tfter.direct_transform (aux);
      for (nat k= 0; k < n ; k++)
        dest[k](i,j)= aux[k];
    }
  mmx_delete<F> (aux, m); }

template<typename Op,typename C, typename F, typename MV, typename Tfter>
static void
mat_inverse_tft (C* d, nat d_rs, nat d_cs, nat r, nat c,
                 const matrix<F,MV>* s, Tfter& tfter) {
  typedef typename Op::nomul_op Acc;
  typedef typename Vector_variant(F) VV;
  nat n= tfter.len;
  nat m= next_power_of_two (n);
  nat l= aligned_size<F,VV> (m);
  for (nat i= 0; i < r; i++)
    for (nat j= 0; j < c; j++) {
      F* aux= mmx_new<F> (l);
      for (nat k= 0; k < n ; k++)
        aux[k]= s[k](i,j);
      tfter.inverse_transform (aux);
      vector<F,VV> tmp (new vector_rep<F,VV> (aux, n, l, false, format<F> ()));
      Acc::set_op (d[i * d_rs + j * d_cs], as<C> (tmp)); } }

template<typename Op, typename C, typename Tfter> static void
mul (C* d, const C* s1, const C* s2,
     nat r, nat rr, nat l, nat ll, nat c, nat cc, Tfter& tfter) {
  typedef Scalar_type (C) F;
  typedef matrix<F> Matrix_F;
  nat n= tfter.len;
  Matrix_F* mm1= mmx_new<Matrix_F> (n);
  Matrix_F* mm2= mmx_new<Matrix_F> (n);
  Matrix_F* mmd= mmx_new<Matrix_F> (n);
  mat_direct_tft (mm1, s1, Mat::index (1, 0, rr, ll), Mat::index (0, 1, rr, ll),
                  r, l, tfter);
  mat_direct_tft (mm2, s2, Mat::index (1, 0, ll, cc), Mat::index (0, 1, ll, cc),
                  l, c, tfter);
  for (nat k= 0; k < n; k++)
    mmd[k]= mm1[k] * mm2[k];
  mat_inverse_tft<Op> (d, Mat::index (1, 0, rr, cc), Mat::index (0, 1, rr, cc),
                       r, c, mmd, tfter);
  mmx_delete<Matrix_F> (mm1, n);
  mmx_delete<Matrix_F> (mm2, n);
  mmx_delete<Matrix_F> (mmd, n); }

template<typename Op, typename C> static void
mul (C* d, const C* s1, const C* s2,
     nat r, nat rr, nat l, nat ll, nat c, nat cc) {
  typedef matrix_tft_multiply_helper<C> Matrix_tft;
  typedef typename Matrix_tft::tft_transformer Tfter;
  typedef typename Matrix_tft::F F;
  typedef typename Matrix_variant(F) MVF;
  typedef implementation<matrix_multiply,MVF> MatF;
  nat sz= matrix_tft_multiply_helper<C>
    ::size (s1, Mat::index (1, 0, rr, ll), Mat::index (0, 1, rr, ll),
            s2, Mat::index (1, 0, ll, cc), Mat::index (0, 1, ll, cc),
            r, l, c);
  if (sz == 0)
    Mat::template clr<Op> (d, r, rr, c, cc);
  else if (sz == 1) {
    nat l1= aligned_size<F,MVF> (r * l);
    nat l2= aligned_size<F,MVF> (l * c);
    nat ld= aligned_size<F,MVF> (r * c);
    F* m1= mmx_new<F> (l1);
    F* m2= mmx_new<F> (l2);
    F* md= mmx_new<F> (ld);
    Mat::template mat_binary_scalar_stride<access_op>
      (m1, 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), 0, r, l);
    Mat::template mat_binary_scalar_stride<access_op>
      (m2, 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), 0, l, c);
    Mat::template mat_binary_scalar_stride<access_op>
      (md, Mat::index (1, 0, r, c), Mat::index (0, 1, r, c),
       d, Mat::index (1, 0, rr, cc), Mat::index (0, 1, rr, cc), 0, r, c);
    MatF::template mul<Op> (md, m1, m2, r, r, l, l, c, c);
    Mat::template mat_unary_stride<id_op>
      (d, Mat::index (1, 0, rr, cc), Mat::index (0, 1, rr, cc),
       md, Mat::index (1, 0, r, c), Mat::index (0, 1, r, c), r, c);
    mmx_delete<F> (m1, l1);
    mmx_delete<F> (m2, l2);
    mmx_delete<F> (md, ld);
  }
  else {
    Tfter tfter (sz, format<F> ());
    mul<Op> (d, s1, s2, r, rr, l, ll, c, cc, tfter); } }

template<typename C> static void
mul (C* d, const C* s1, const C* s2, nat r, nat l, nat c) {
  mul<mul_op> (d, s1, s2, r, r, l, l, c, c); }

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::template mul<Op> (d, s1, s2, r, rr, l, ll, c, cc); }

template<typename D, typename S1, typename S2>
static inline void
mul (D* d, const S1* s1, const S2* s2,
     nat r, nat l, nat c) {
  Mat::template mul<mul_op> (d, s1, s2, r, r, l, l, c, c); }

}; // implementation<matrix_multiply_base,V,matrix_tft>

} // namespace mmx
#endif // __MMX_MATRIX_TFT_HPP

---