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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_complex.hpp
* DESCRIPTION: Complex multiplication of matrices
* COPYRIGHT  : (C) 2008 Raouf Dridi and 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_COMPLEX__HPP
#define __MMX__MATRIX_COMPLEX__HPP
#include <algebramix/matrix.hpp>
#include <numerix/complex.hpp>
#include <numerix/integer.hpp>

namespace mmx {
#define TMPL template<typename C>
#define Complex complex<C>

/******************************************************************************
* Complex matrix variant
******************************************************************************/

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

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

TMPL
struct matrix_variant_helper<Complex > {
  typedef typename matrix_variant_helper<C>::MV CMV;
  typedef matrix_complex<CMV> MV;
};

STMPL
struct matrix_variant_helper<complex<generic> > {
  typedef matrix_naive MV;
};

/******************************************************************************
* Complex matrix multiplication
******************************************************************************/

template<typename V, typename CV>
struct implementation<matrix_multiply,V,matrix_complex<CV> >:
  public implementation<matrix_linear,V>
{
  typedef implementation<vector_linear,CV> Vec_C;
  typedef implementation<matrix_multiply,CV> Mat_C;

TMPL static inline void
mul (Complex* dest, const Complex* s1, const Complex* s2,
     nat r, nat l, nat c)
{
  nat n1 = r*l, spc1= aligned_size<C,CV> (n1);
  nat n2 = l*c, spc2= aligned_size<C,CV> (n2);
  nat nd = r*c, spcd= aligned_size<C,CV> (nd);
  nat spc= 2*spc1 + 2*spc2 + 3*spcd;
  C* real_x1= mmx_new<C> (spc);
  C* imag_x1= real_x1 + spc1;
  C* real_x2= imag_x1 + spc1;
  C* imag_x2= real_x2 + spc2;
  C* real_d = imag_x2 + spc2;
  C* imag_d = real_d  + spcd;
  C* temp_d = imag_d  + spcd;

  Vec_C::half_copy (real_x1, (C*) ((void*) s1), n1);
  Vec_C::half_copy (imag_x1, ((C*) ((void*) s1)) + 1, n1);
  Vec_C::half_copy (real_x2, (C*) ((void*) s2), n2);
  Vec_C::half_copy (imag_x2, ((C*) ((void*) s2)) + 1, n2);
  //mmout << "real_x1= "; Vec_C::print (mmout, real_x1, n1); mmout << "\n";
  //mmout << "real_x2= "; Vec_C::print (mmout, real_x2, n2); mmout << "\n";
  Mat_C::mul (real_d, real_x1, real_x2, r, l, c);
  //mmout << "mul1_d = "; Vec_C::print (mmout, real_d , nd); mmout << "\n";
  //mmout << "imag_x1= "; Vec_C::print (mmout, imag_x1, n1); mmout << "\n";
  //mmout << "imag_x2= "; Vec_C::print (mmout, imag_x2, n2); mmout << "\n";
  Mat_C::mul (temp_d, imag_x1, imag_x2, r, l, c);
  //mmout << "mul2_d = "; Vec_C::print (mmout, temp_d , nd); mmout << "\n";
  Vec_C::add (real_x1, imag_x1, n1);
  Vec_C::add (real_x2, imag_x2, n2);
  //mmout << "sum_x1 = "; Vec_C::print (mmout, real_x1, n1); mmout << "\n";
  //mmout << "sum_x2 = "; Vec_C::print (mmout, real_x2, n2); mmout << "\n";
  Mat_C::mul (imag_d, real_x1, real_x2, r, l, c);
  //mmout << "mul3_d = "; Vec_C::print (mmout, imag_d , nd); mmout << "\n";
  Vec_C::sub (imag_d, real_d, nd);
  Vec_C::sub (imag_d, temp_d, nd);
  Vec_C::sub (real_d, temp_d, nd);
  Vec_C::double_copy ((C*) ((void*) dest), real_d, nd);
  Vec_C::double_copy (((C*) ((void*) dest)) + 1, imag_d, nd);

  mmx_delete<C> (real_x1, spc);
}

}; // implementation<matrix_multiply,V,matrix_complex<CV> >

#undef TMPL
#undef Complex
} // namespace mmx
#endif //__MMX__MATRIX_COMPLEX__HPP