source: Sophya/trunk/SophyaLib/NTools/toeplitzMatrix.h@ 3440

#ifndef TOEPLITZMATRIX_SEEN
#define TOEPLITZMATRIX_SEEN



// matrice de Toeplitz reelle ou complexe
//------------------------------------------------------------------
//   *****   Guy Le Meur -- LAL-Orsay mars 2002 ****************
//-------------------------------------------------------------------


#include "machdefs.h"    // Definitions specifiques SOPHYA

#include <math.h>
#include <iostream>

#include "nbmath.h"
#include "timing.h"

#include "array.h"

#include "fftservintf.h"
#include "fftpserver.h"

// classe pour decrire une matrice de Toeplitz reelle ou complexe, 
// symetrique (hermitienne) ou non. 
// Methodes de gradient conjugues pour resolution de systemes (uniquement
// pour symetriques ou hermitiennes, pour le moment)

namespace SOPHYA {

class Toeplitz
{

 private: 
  // verouiller le clonage
Toeplitz(const Toeplitz&)  {}
Toeplitz &operator = (const Toeplitz&) {return *this;}

public:

Toeplitz();
Toeplitz(const TVector<complex<double> >& firstCol);

~Toeplitz();

  // toeplitz complexe generale
  void  setMatrix(const TVector<complex<double> >& firstCol, const TVector<complex<double> >& firstRow); 
  // toeplitz hermitienne
  void  setMatrix(const TVector<complex<double> >& firstCol); 
  // toeplitz reelle generale
  void  setMatrix(const TVector<double>& firstCol, const TVector<double>& firstRow); 
  // toeplitz reelle symetrique
  void  setMatrix(const TVector<double>& firstCol); 


  int gradientToeplitz(TVector<complex<double> >& b) const;
  int gradientToeplitz(TVector<double>& b) const;

  int gradientToeplitzPreconTChang(TVector<complex<double> >& b) const;

  int CGSToeplitz(TVector<double>& b) const;
  int DCGToeplitz(TVector<double>& b) ;

private:

inline complex<double> prodScalaire(int n,  TVector<complex<double> >& v1, TVector<complex<double> >& v2) const
{
  complex<double> produit = complex<double>(0.,0.);
  for (int k =0; k<n; k++) produit += v1(k)*conj(v2(k));
  return produit;
}
inline double prodScalaire(int n,  TVector<double>& v1, TVector<double>& v2) const
{
  double produit = 0.;
  for (int k =0; k<n; k++) produit += v1(k)*v2(k);
  return produit;
}


// Reza+CMV du 3/7/2006 : on a instancie le FFTPackServer avec le flag preserveinput=true
// On peut donc faire en principe les const_cast sans danger 
inline void transformeeFourier(const TVector<complex<double> >& v, TVector< complex<double> >& tv) const
{
  fftIntfPtr_-> FFTForward(const_cast<TVector<complex<double> > & >(v), tv);
} 
inline void transformeeFourier(const TVector<double>& v, TVector< complex<double> >& tv) const
{
  fftIntfPtr_-> FFTForward(const_cast<TVector<double> & >(v), tv);
} 


inline void transformeeInverseFourier(const TVector<complex<double> >& tv, TVector<complex<double> >& v) const
{
  int n= tv.Size();
  fftIntfPtr_-> FFTBackward(const_cast<TVector<complex<double> > & >(tv), v);
  v/=n;
} 
inline void transformeeInverseFourier(const TVector<complex<double> >& tv, TVector<double>& v) const
{
  int n= tv.Size();
  fftIntfPtr_-> FFTBackward(const_cast< TVector<complex<double> > & >(tv), v);
  double fac = 1./(2*(n-1));
  v *= fac;
} 

inline void initTFTransposeeToeplitzComplexe() 
{
  int k;
  int ndyad = vecteurCirculant_.Size();
  TVector<complex<double> > vecteurCirculantTranspose(ndyad);
  vecteurCirculantTranspose =  complex<double>(0.,0.);
  vecteurCirculantTranspose(0) = vecteurCirculant_(0);
  for (k=1; k< dimTop_; k++)
    {
      vecteurCirculantTranspose(k) = vecteurCirculant_(ndyad-k);
      vecteurCirculantTranspose(ndyad-k) = vecteurCirculant_(k);
    }

  transformeeFourier(vecteurCirculantTranspose, CirculanteTransposeeFourier_);
}
inline void initTFTransposeeToeplitzReelle() 
{
  int k;
  int ndyad = vecteurCirculantD_.Size();
  TVector<double> vecteurCirculantTranspose(ndyad);
  vecteurCirculantTranspose =  0.;
  vecteurCirculantTranspose(0) = vecteurCirculantD_(0);
  for (k=1; k< dimTop_; k++)
    {
      vecteurCirculantTranspose(k) = vecteurCirculantD_(ndyad-k);
      vecteurCirculantTranspose(ndyad-k) = vecteurCirculantD_(k);
    }

  transformeeFourier(vecteurCirculantTranspose, CirculanteTransposeeFourier_);
}

void extensionACirculanteDyadique(const TVector<complex<double> >& firstCol, const TVector<complex<double> >& firstRow);

void extensionACirculanteDyadique(const TVector<double>& firstCol, const TVector<double>& firstRow);
void extensionACirculanteDyadique(const TVector<complex<double> >& firstCol);

void extensionACirculanteDyadique(const TVector<double>& firstCol);


// matrice circulante entree par sa T. de Fourier
void produitParVecFourier(const TVector<complex<double> >& q, TVector<complex<double> >& Tq) const;
void produitParVecFourier(const TVector<double>& q, TVector<double>& Tq) const;

void produitTransposeeParVecFourier(const TVector<double>& q, TVector<double>& Tq) const;


// preconditionneur pour une Toeplitz HERMITIENNE
// LE RESUTAT EST UNE TRANSFORMEE DE FOURIER
void fabricationTChangPreconHerm(TVector<complex<double> >& TFourierC) const;

// matrice circulante entree par sa T. de Fourier
void inverseSystemeCirculantFourier(const TVector<complex<double> >& TFourierC, const TVector<complex<double> >& secondMembre, TVector<complex<double> >& resul) const;


void expliciteCirculante(TMatrix<complex<double> >& m) const;
void expliciteToeplitz(TMatrix<complex<double> >& m) const;
void expliciteToeplitz(TMatrix<double>& m) const;


  bool hermitian_;
  int dimTop_;
  TVector<complex<double> > vecteurCirculant_;
  TVector<double> vecteurCirculantD_;
  TVector<complex<double> > CirculanteFourier_;
  TVector<complex<double> > CirculanteTransposeeFourier_;
  FFTServerInterface* fftIntfPtr_;
};


// fin classe Toeplitz


}  // Fin du namespace

#endif