synaps/linalg/hankel.h

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.        *
*   Author(s): B. Mourrain, GALAAD, INRIA                            *
**********************************************************************
History: 
$Id: hankel.h,v 1.1 2005/07/11 08:33:35 mourrain Exp $
**********************************************************************/
#ifndef Hankel_H
#define Hankel_H

#include <synaps/init.h>
#include <synaps/linalg/rep1d.h>
#include <synaps/base/Range.h>
#include <synaps/linalg/FFT.h>

__BEGIN_NAMESPACE_SYNAPS


template<char O,class A, class B> struct OP;
//OPEN
namespace linalg {
//----------------------------------------------------------------------
template<class C> 
struct hankel : public linalg::rep1d<C> {

  unsigned int nbrow_, nbcol_;

  typedef C           coeff_t;
  typedef typename linalg::rep1d<C>::size_type   size_type;
  typedef rep1d<C>  row_t;
  typedef rep1d<C>  col_t;

  hankel(): rep1d<C>(){}
  hankel(const hankel<C> & v): rep1d<C>(v), 
    nbrow_(v.nbrow_), nbcol_(v.nbcol_){}
  hankel(unsigned int i, unsigned int j): 
    rep1d<C>(i+j-1),nbrow_(i),nbcol_(j) {}
  hankel(unsigned int i, unsigned int j, C* t): 
    rep1d<C>(i+j-1,t),nbrow_(i),nbcol_(j){}

  void resize(unsigned int i, unsigned int j)
    {
      this->rep1d<C>::resize(i+j-1); nbrow_ = i; nbcol_=j;
    }

  //  C &   operator()(unsigned int i, unsigned int j) {return tab_[i+j];} 
  C   operator()(unsigned int i, unsigned int j) const 
        {return this->tab_[i+j];} 

  size_type nbrow() const {return nbrow_;}
  size_type nbcol() const {return nbcol_;}

  void transpose() { std::swap(nbcol_,nbrow_);};
};
//----------------------------------------------------------------------
// template<class C>
// void assign(hankel<C> & r, const hankel<C>* v) 
// {
//   assign((rep1d<C> &)r, (rep1d<C>*) v);
//   r.nbrow_=v->nbrow_;
//   r.nbcol_=v->nbcol_;
// }
//----------------------------------------------------------------------

template<class C, class M>
void submatrix( hankel<C> & r, const M & m, const Range2d  & rg )
{
  typedef typename hankel<C>::size_type size_type;
  for(size_type i=rg.j1; i< rg.j2+1; i++)
    r[i-rg.j1]=m(rg.i1,i);
  for(size_type i=rg.i1+1; i< rg.i2+1; i++)
    r[i+rg.j2-rg.j1-rg.i1]=m(i,rg.j2);
}

  //-----------------------------------------------------------------
  // Assignment of a matrix to a container {\tt hankel}
  template<class S, class C>
  inline void assign(hankel<C> & r, S M)
  {
    for(unsigned int i=0;i<M.nbcol();i++) r(0,i)=M(0,i);  
    unsigned int c= M.nbcol()-1;
    for(unsigned int i=1;i<M.nbrow();i++) r(i,c)=M(i,c);  
  }
  //----------------------------------------------------------------------
  template<class C,class V> 
  void mult_vect(V & r, const hankel<C> & t, const V & v);
  
} //CLOSE linalg

__END_NAMESPACE_SYNAPS

//-------------------------------------------------------------------
#ifndef SEPARATE_COMPILATION
#include "synaps/linalg/hankel.C"
#endif
//-------------------------------------------------------------------

#endif //hankel_H