synaps/linalg/MatrStr.h

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.        *
*   Author(s): B. Mourrain, GALAAD, INRIA                            *
**********************************************************************
History: 
$Id: MatrStr.h,v 1.1 2005/07/11 08:33:34 mourrain Exp $
**********************************************************************/
#ifndef SYNAPS_LINALG_MATRSTR_H
#define SYNAPS_LINALG_MATRSTR_H
//----------------------------------------------------------------------
#include <sstream>
#include <synaps/init.h>
#include <synaps/base/Range.h>
#include <synaps/arithm/traits.h>
#include <synaps/arithm/texp.h>
#include <synaps/linalg/VECTOR.m>
#include <synaps/linalg/MATRIX.m>
#include <synaps/linalg/VectDse.h>
#include <synaps/linalg/texp.h>

__BEGIN_NAMESPACE_SYNAPS

//----------------------------------------------------------------------
template<class R> struct MatrStr  {
  
  typedef VectDse<typename R::value_type> vector_type;
  typedef vector_type col_type;
  typedef vector_type row_type;
  shared_object<R>  data;
  R &       rep()       {return (*data);}
  R  rep() const {return (*data);}

  typedef typename R::size_type  size_type;
  typedef typename R::value_type value_type;
  typedef MatrStr<R>             self_t;

  MatrStr(): data() {}
  MatrStr(const R  & r): data(r) {}
  MatrStr(int m, int n) {rep().resize(m,n);}
  MatrStr(int m, int n, value_type* t): data(R(m,n,t)) {}
  MatrStr(unsigned m, unsigned n, char *t);

  //Copy constructor
  MatrStr(const self_t & M): data(M.data) {} 
  //Assignement
  self_t & operator=(const self_t & M)  {rep()=M.rep(); return *this;} 
  template<class X>
  self_t & operator=(const X & M);

  size_type nbrow() const {return rep().nbrow();}
  size_type nbcol() const {return rep().nbcol();}

  value_type     operator()(size_type i, size_type j) const {return rep()(i,j);} 
  
  typedef typename binary_operator_prototype< linalg::_submatrix_, self_t, Range2d >::F submatrix_t;
  submatrix_t  operator()(const Range & I, const Range & J) 
  { return submatrix(*this,Range2d(I,J)); };
  
  template<class X>
  self_t & operator+=(const X & x) { return *this = *this + x; };
  template<class X>
  self_t & operator-=(const X & x) { return *this = *this - x; };
  template<class X>
  self_t & operator*=(const X & x) { return *this = *this * x; };
  
  self_t & transpose() {rep().transpose(); return *this;}
  
};

template<class M, class S> inline
void submatrix( M & r, const MatrStr<S> & a, const Range2d & b )
{
  std::cout << a << " " << b << std::endl;
  using namespace MATRIX;
  submatrix(r.rep(),a.rep(),b);
};

template<class V, class R> inline
void row( V & v, const MatrStr<R> & m, int i )
{
  using MATRIX::getrow;
  v.resize( m.nbcol() );
  getrow( v.rep(), m.rep(), i );
};
template<class V, class R> inline
void col( V & v, const MatrStr<R> & m, int i )
{
  using MATRIX::getcol;
  v.resize( m.nbrow() );
  getcol( v.rep(), m.rep(), i );
};
 
template<class S> inline
void col(  MatrStr<S>& sm, const MatrStr<S> & m, const Range & rg )
{
  submatrix(sm,m,Range2d(0,m.nbrow()-1,rg.i1,rg.i2));
};

template<class S> inline
void row(  MatrStr<S>& sm, const MatrStr<S> & m, const Range & rg )
{
  submatrix(sm,m,Range2d(rg.i1,rg.i2,0,m.nbcol()-1));
};



template<class R,class S> inline
void add( MatrStr<R> & r, const MatrStr<S> & a, const MatrStr<S>  &  b )
{
  r.rep().resize(a.nbrow(),b.nbcol());
  using namespace VECTOR; 
  add(r.rep(),a.rep(),b.rep());
};

template<class R> inline
void sub(MatrStr<R> & r, const MatrStr<R> & a ,const MatrStr<R> & b )
{
  r.rep().resize(a.nbrow(),b.nbcol());
  using namespace VECTOR; 
  sub(r.rep(),a.rep(),b.rep());
}

template<class S> inline
void mul(MatrStr<S> & r, const MatrStr<S> & a, const typename MatrStr<S>::value_type & b )
{
  r.rep().resize(a.nbrow(),b.nbcol());
  using namespace VECTOR; mul_ext(r.rep(),a.rep(),b);
};
template<class R,class S,class C> inline
void mul(MatrStr<R> & r, const MatrStr<S> & a, const C & b )
{
  r.rep().resize(a.nbrow(),b.nbcol());
  using namespace VECTOR; mul_ext(r.rep(),a.rep(),b);
}

template<class R,class S,class C> inline
void mul(MatrStr<R> & r, const C & b, const MatrStr<S> & a  )
{
  mul(r,a,b);
}

template<class Ca, class Ra,
         class S,
         class Cb, class Rb>
void mul(VectDse<Ca,Ra> & r, const MatrStr<S> & m, const VectDse<Cb,Rb> & v )
{
  assert(m.nbcol()==v.size());
  r.resize(m.nbrow());
  using namespace MATRIX;
  mul_vect(r.rep(),m.rep(),v.rep());
};
//======================================================================
// ASSIGNEMENTS
//======================================================================
template<class C,class S> struct MatrDse;
namespace let {
  
  template<class R,class C, class S> inline
  void assign(MatrStr<R> & r, const MatrDse<C,S> & M) 
  { 
    using MATRIX::assign;
    assign(r.rep(), M.rep()); 
  }
  
  //----------------------------------------------------------------------
  template<class R> inline
  void assign(MatrStr<R> & r, const MatrStr<R>* M)
  {
    assign(r.rep(), &M->rep());
  }
}

template<class R>
std::ostream & operator<<(std::ostream & os, const  MatrStr<R> & p)
{
  using namespace MATRIX; return print(os,p.rep());
}
//----------------------------------------------------------------------
template<class R>  inline
std::istream & operator>>(std::istream & is, MatrStr<R> & M)
{
  return read(is,M.rep());
}
//======================================================================

namespace type
{
  template<class C, class S>
    struct instanceof< MatrStr<S>, C > { typedef MatrStr< typename instanceof<S,C>::T > T; };
};
namespace arithm 
{
  template<class S> struct structureof< MatrStr<S> > { typedef structure::matrix T; };

#define head  template<class R>
#define parm0 MatrStr<R>
#define parm1 MatrStr<R>
  declare_binary_operator(head,parm0,int,linalg::_col_,col);
  declare_binary_operator(head,parm0,int,linalg::_row_,row);
  declare_binary_operator(head,parm0,Range,linalg::_col_,col);
  declare_binary_operator(head,parm0,Range,linalg::_row_,row);
  declare_binary_operator(head,parm0,Range2d,linalg::_submatrix_,submatrix);
  declare_binary_operator(head,parm0,parm1,_add_,operator+);
  declare_binary_operator(head,parm0,parm1,_sub_,operator-);
  declare_binary_operator(head,parm0,parm1,_mul_,operator*);
  declare_unary_operator(head,parm0,_neg_,operator-);
#undef head
#undef parm0
#undef parm1
#define head  template<class R, class C, class Rc>
#define parm0 MatrStr<R>
#define parm1 VectDse<C,Rc>
  declare_binary_operator(head,parm0,parm1,_mul_,operator*);
#undef head
#undef parm0
#undef parm1
#define head  template<class R>
#define parm0 MatrStr<R>
#define parm1 typename MatrStr<R>::value_type
  declare_binary_operator(head,parm0,parm1,_mul_,operator*);
  declare_binary_operator(head,parm1,parm0,_mul_,operator*);
  declare_binary_operator(head,parm0,parm1,_div_,operator/);
#undef head
#undef parm0
#undef parm1
  

}

template<class S>
template<class X>
MatrStr<S> & MatrStr<S>::operator=( const X & x )
{
  using let::assign; assign(*this,x); 
  return *this;
};




__END_NAMESPACE_SYNAPS

#ifndef SEPARATE_COMPILATION
#include "synaps/linalg/MatrStr.C"
#endif
//======================================================================

#endif // SYNAPS_LINALG_MATRSTR_H