source: Sophya/trunk/SophyaLib/NTools/toeplitzMatrix.cc@ 4068

#include "sopnamsp.h"
#include "toeplitzMatrix.h"
#include "timing.h"   

/*!
  \class SOPHYA::Toeplitz
  \ingroup NTools
  Class for describing a real Toeplitz or complex matrix.
  Symmetic (/hermitian) or not.
  Conjugate gradient system solver (for symmetric or hermitian matrices
  only)
  Toeplitz matrices have  constant values along negative slope diagonals.
  A general toeplitz nxn matrix is defined by 2n-1 values and a symmetric
  (/hermitian) by n values. 
  General toeplitz matrix form :     
\verbatim
  a_0   b_1   b_2   ...     b_n-1
  a1   a0   b1   b2  ... 
  a2   a1   a0   b1   ...
   .....
   .....
  a_n-1  ....                a0
\endverbatim
*/


Toeplitz::Toeplitz() : hermitian_(false)
  {
    fftIntfPtr_=new FFTPackServer(true);  // true -> on ne surecrit pas les input
    fftIntfPtr_->setNormalize(false);
  }

// matrice de Toplitz hermitienne. La premiere ligne est conjuguee de la 
// premeire colonne donnee (vecteur firstCol)
Toeplitz::Toeplitz(const TVector<complex<double> >& firstCol) : hermitian_(false)
  {
    fftIntfPtr_=new FFTPackServer(true); // true -> on ne surecrit pas les input
    fftIntfPtr_->setNormalize(false);
    setMatrix(firstCol); 
  }


Toeplitz::~Toeplitz(){ if (fftIntfPtr_!=NULL) delete fftIntfPtr_;};


// initialise matrice de Toeplitz hermitienne
void  Toeplitz::setMatrix(const TVector<complex<double> >& firstCol) 
  {
    hermitian_ = true;
    dimTop_ = firstCol.Size();
    extensionACirculanteDyadique(firstCol);
    //    transformeeFourier(vecteurCirculant_, CirculanteFourier_);
  }

// intialise la matrice de Toeplitz generale complexe
void  Toeplitz::setMatrix(const TVector<complex<double> >& firstCol, const TVector<complex<double> >& firstRow) 
  {
    hermitian_ = false;
    dimTop_ = firstCol.Size();
  if (dimTop_ != firstRow.Size())
    {
      cout << " Toeplitz::setMatrix : nb col different de nbLignes " << endl; 
    }
    extensionACirculanteDyadique(firstCol, firstRow);
    //    transformeeFourier(vecteurCirculant_, CirculanteFourier_);
  }


// intialise la matrice de Toeplitz generale reelle
void  Toeplitz::setMatrix(const TVector<double>& firstCol, const TVector<double>& firstRow) 
  {
    hermitian_ = false;
    dimTop_ = firstCol.Size();
  if (dimTop_ != firstRow.Size())
    {
      cout << " Toeplitz::setMatrix : nb col different de nbLignes " << endl; 
    }
    extensionACirculanteDyadique(firstCol, firstRow);
    //    transformeeFourier(vecteurCirculantD_, CirculanteFourier_);
  }

// intialise la matrice de Toeplitz generale reelle symetrique
void  Toeplitz::setMatrix(const TVector<double>& firstCol) 
  {
    hermitian_ = true;
    dimTop_ = firstCol.Size();
    extensionACirculanteDyadique(firstCol);
    //    transformeeFourier(vecteurCirculantD_, CirculanteFourier_);
  }


// resolution par gradient conjugue d'un systeme dont la matrice est 
// la matrice de Toeplitz
// (le vecteur b est le second membre)
// renvoie le nombre d'iterations (-1 si non convergence)
int Toeplitz::gradientToeplitz(TVector<complex<double> >& b) const
{
  // PrtTim(" entree gradient " );
  int k;
  if (dimTop_ != b.Size())
    {
      throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n");
    }
  if (!hermitian_) 
    {
      throw PException(" TOEPLITZMATRIX : conjugate gradient for hermitian matrix only\n");
    }
  int ndyad = vecteurCirculant_.Size();

  TVector<complex<double> > a(dimTop_);
  TVector<complex<double> > r(dimTop_);
  TVector<complex<double> > q(ndyad);

  a =complex<double>(0.,0.);
  
  r = b;
  q =complex<double>(0.,0.);
  // int dimTopM1 = dimTop_-1;
   for (k=0; k<dimTop_; k++) q(k) = r(k);
  //q(Range(0,dimTopM1)) = r;
  //q.Print(cout, q.Size());  

  double norm2Rk = prodScalaire(dimTop_,r,r).real();  
  double norm2R0 = norm2Rk;  
  int niter =0;
    do 
    {
      int k;
      TVector<complex<double> > Aq;
      produitParVecFourier(q, Aq);
      //
      double  norme2Aq = prodScalaire(dimTop_,Aq,q).real();
      complex<double> alpha = prodScalaire(dimTop_,r,r)/norme2Aq;
      for (k=0; k<dimTop_; k++) a(k) += alpha*q(k);
      for (k=0; k<dimTop_; k++) r(k) -= alpha*Aq(k);
      double norm2Rkp1 = prodScalaire(dimTop_,r,r).real(); 
      double beta = norm2Rkp1/norm2Rk;
      for (k=0; k<dimTop_; k++) q(k) = r(k) + beta*q(k);

      norm2Rk = norm2Rkp1;
      
      //  cout << " iteration " << niter+1 << " norme residu= " << norm2Rk << " critere arret " << norm2Rk/norm2R0 << endl;
      niter++;
    }

    while (norm2Rk/norm2R0 > 1.e-12 && niter < 2*dimTop_);

    //   PrtTim(" fin du gradient ");

    // verification a supprimer 
    //TMatrix<complex<double> > ttest(dimTop_, dimTop_);
    //expliciteToeplitz(ttest);
    //TVector<complex<double> > btest;
    //ttest.Show();
    //cout << " long vec circ " << dimTop_ << " long a " << a.Size() << endl;
    //btest = ttest*a;
    //for (int k=0; k< b.Size(); k++)
    // {
    //  cout << " verif gradient donnee " << b(k) << " retrouve " << btest(k) << endl;
    //  }
    //
    b = a;
    
    if (niter == 2*dimTop_-1) niter = -1;
    return niter;
}


// resolution par gradient conjugue d'un systeme dont la matrice est 
// la matrice de Toeplitz
// (le vecteur b est le second membre)
// renvoie le nombre d'iterations (-1 si non converge
int Toeplitz::gradientToeplitz(TVector<double>& b) const
{
  int k;
  if (dimTop_ != b.Size())
    {
      throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n");
    }
  if (!hermitian_) 
    {
      throw PException(" TOEPLITZMATRIX : conjugate gradient for hermitian matrix only\n");
    }

  int ndyad = vecteurCirculantD_.Size();

  TVector<double> a(dimTop_);
  TVector<double> r(dimTop_);
  TVector<double> q(ndyad);

  a =0.;

  r = b;
  q =0.;
  for (k=0; k<dimTop_; k++) q(k) = r(k);


  double norm2Rk = prodScalaire(dimTop_,r,r);  
  double norm2R0 = norm2Rk;  
  int niter =0;
  do 
    {
      int k;
      TVector<double> Aq;
      produitParVecFourier(q, Aq);

      //
      double  norme2Aq = prodScalaire(dimTop_,Aq,q);
      double alpha = prodScalaire(dimTop_,r,r)/norme2Aq;
      for (k=0; k<dimTop_; k++)
        {
          a(k) += alpha*q(k);
          r(k) -= alpha*Aq(k);
        }
      double norm2Rkp1 = prodScalaire(dimTop_,r,r); 
      double beta = norm2Rkp1/norm2Rk;
      for (k=0; k<dimTop_; k++) q(k) = r(k) + beta*q(k);

      norm2Rk = norm2Rkp1;
      
      //  cout << " iteration " << niter+1 << " norme residu= " << norm2Rk << " critere arret " << norm2Rk/norm2R0 << endl;
      niter++;
        }

      while (norm2Rk/norm2R0 > 1.e-12 && niter < 2*dimTop_);
    //   while (norm2Rk/norm2R0 > 1.e-12 && niter < 1);



    // verification a supprimer 
  //TMatrix<double> ttest(dimTop_, dimTop_);
  //expliciteToeplitz(ttest);
  // TVector<double> btest;
  // ttest.Show();
    //cout << " long vec circ " << dimTop_ << " long a " << a.Size() << endl;
  // btest = ttest*a;
  // for (int k=0; k< b.Size(); k++)
  //  {
  //    cout << " verif gradient donnee " << b(k) << " retrouve " << btest(k) << endl;
  //   }
    //
    b = a;
    
    if (niter == 2*dimTop_-1) niter = -1;
    return niter;
}
// resolution par Gradient Conjugue Generalise d'un systeme dont 
// la matrice est la matrice de Toeplitz
// (le vecteur b est le second membre)
// renvoie le nombre d'iterations (-1 si non converge
int Toeplitz::CGSToeplitz(TVector<double>& b) const
{
  int k;
  if (dimTop_ != b.Size())
    {
      throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n");
    }

  int ndyad = vecteurCirculantD_.Size();

  TVector<double> x(dimTop_);
  TVector<double> r(dimTop_);
  TVector<double> r0(dimTop_);
  TVector<double> q(ndyad);
  TVector<double> p(ndyad);

  // initialisations 
  x =0.;
  r0 = b;
  r = r0;
  q =0.;
  p =0.;
  for (k=0; k<dimTop_; k++)
    {
      q(k) = r(k);
      p(k) = r(k);
    }
  int niter =0;
  double r0rk = prodScalaire(dimTop_,r0,r);
  double r0r0 = r0rk;
  double rkrk = r0r0;
  double r0rkp1 = r0rk;
  do 
    {
      int k;
      TVector<double> Aq;
      produitParVecFourier(q, Aq);
      TVector<double> Ap;
      produitParVecFourier(p, Ap);

      TVector<double> pMalphaAq(ndyad);


      //

      double  r0Aqk = prodScalaire(dimTop_,r0,Aq);
      double alpha = r0rk/r0Aqk;
      pMalphaAq = p-alpha*Aq;
      TVector<double> aux(ndyad);
      aux = p + pMalphaAq;
      TVector<double> Aaux;
      produitParVecFourier(aux, Aaux);

      for (k=0; k<dimTop_; k++)
        {
          x(k) += alpha*aux(k);
          r(k) -= alpha*Aaux(k);
        }
      r0rkp1 = prodScalaire(dimTop_,r0,r); 
      double beta = r0rkp1/r0rk;
      for (k=0; k<dimTop_; k++) p(k) = r(k) + beta*pMalphaAq(k);
      for (k=0; k<dimTop_; k++)
        {
          q(k)  = p(k) + beta*(beta*q(k) + pMalphaAq(k));
        }
      r0rk  = r0rkp1;
      rkrk = prodScalaire(dimTop_,r,r); 
      //    cout << " iteration " << niter+1 << " rkrk= " << rkrk << " critere arret " << rkrk/r0r0 << endl;
      niter++;
        }

     while (rkrk/r0r0 > 1.e-12 && niter < 2*dimTop_);
  // while ( niter < 10*dimTop_);

    // verification a supprimer 
  //TMatrix<double> ttest(dimTop_, dimTop_);
  //expliciteToeplitz(ttest);
  //TVector<double> btest;
  // ttest.Show();
    //cout << " long vec circ " << dimTop_ << " long a " << a.Size() << endl;
  //btest = ttest*x;
  // for (int k=0; k< b.Size(); k++)
  // {
  //    cout << " verif gradient donnee " << b(k) << " retrouve " << btest(k) << endl;
  //  }
    //
    b = x;
    
    if (niter == 2*dimTop_) niter = -1;
    return niter;
}
// resolution par Double Gradient Conjugue d'un systeme dont 
// la matrice est la matrice de Toeplitz
// (le vecteur b est le second membre)
// renvoie le nombre d'iterations (-1 si non converge
int Toeplitz::DCGToeplitz(TVector<double>& b) 
{
  int k;
  if (dimTop_ != b.Size())
    {
      throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n");
    }

  // on a besoin des produit (matrice trasposee)Xvectteur
  // on initialise eventuellement la TF de Toeplitz transposee
  if (CirculanteTransposeeFourier_.Size() == 0) initTFTransposeeToeplitzReelle();


  int ndyad = vecteurCirculantD_.Size();

  TVector<double> x(dimTop_);
  TVector<double> r1(dimTop_);
  TVector<double> r2(dimTop_);
  TVector<double> d1(ndyad);
  TVector<double> d2(ndyad);

  // initialisations 
  x =0.;
  r1 = b;
  r2 = r1;
  d1 =0.;
  d2 =0.;
  for (k=0; k<dimTop_; k++)
    {
      d1(k) = r1(k);
      d2(k) = r2(k);
    }
  int niter =0;
  double r1r2k = prodScalaire(dimTop_,r1,r2);
  double r1r2ini = r1r2k;
  do 
    {
      int k;
      TVector<double> Ad1;
      produitParVecFourier(d1, Ad1);
      TVector<double> tAd2;
      produitTransposeeParVecFourier(d2, tAd2);

      //
      double  d2Ad1 = prodScalaire(dimTop_,d2,Ad1);
      double alpha = r1r2k/d2Ad1;

      for (k=0; k<dimTop_; k++)
        {
          x(k) += alpha*d1(k);
          r1(k) -= alpha*Ad1(k);
          r2(k) -= alpha*tAd2(k);
        }
      prodScalaire(dimTop_,r1,r2); 
      double beta = prodScalaire(dimTop_,r1,r2)/r1r2k;
      for (k=0; k<dimTop_; k++)
        {
          d1(k)  = r1(k) + beta*d1(k);
          d2(k)  = r2(k) + beta*d2(k);
        }
      r1r2k = prodScalaire(dimTop_,r1,r2); 
      //      cout << " iteration " << niter+1 << " r1r2= " << r1r2k << " critere arret " << r1r2k/r1r2ini << endl;
      niter++;
        }

     while ( fabs(r1r2k/r1r2ini) > 1.e-12 && niter < 2*dimTop_);

    // verification a supprimer 
  //TMatrix<double> ttest(dimTop_, dimTop_);
  //expliciteToeplitz(ttest);
  //TVector<double> btest;
  // ttest.Show();
    //cout << " long vec circ " << dimTop_ << " long a " << a.Size() << endl;
  //btest = ttest*x;
  //for (int k=0; k< b.Size(); k++)
  // {
  //    cout << " verif gradient donnee " << b(k) << " retrouve " << btest(k) << endl;
  //   }
    //
    b = x;
    
    if (niter == 2*dimTop_) niter = -1;
    return niter;
}


void  Toeplitz::expliciteCirculante(TMatrix<complex<double> >& m) const
{
  int k;
  int n= vecteurCirculant_.Size();
  m.ReSize(n,n);
  int decal=0;
  for (k=0; k<n; k++)
    {
      int j;
      int courant = 0;
      for (j=decal; j< n; j++)
        {
          m(j, k) = vecteurCirculant_(courant++);
        }
      for (j=0; j<decal; j++)
        {
          m(j,k)= vecteurCirculant_(courant++);
        }
      decal++;
    }
}
void  Toeplitz::expliciteToeplitz(TMatrix<complex<double> >& m) const
{
  int k;
  int ndyad = vecteurCirculant_.Size();
  m.ReSize(dimTop_, dimTop_);
  for (k=0; k<dimTop_; k++)
    {
      int j;
      for (j=k; j< dimTop_; j++)
        {
          m(j, k) = vecteurCirculant_(j-k);
        }
      for (j=0; j<k; j++)
        {
          m(k-j-1,k)= vecteurCirculant_(ndyad-j-1);
        }
    }
}
void  Toeplitz::expliciteToeplitz(TMatrix<double>& m) const
{
  int k;  
  int ndyad = vecteurCirculantD_.Size();
  m.ReSize(dimTop_, dimTop_);
  for (k=0; k<dimTop_; k++)
    {
      int j;
      for (j=k; j< dimTop_; j++)
        {
          m(j, k) = vecteurCirculantD_(j-k);
        }
      for (j=0; j < k; j++)
        {
          m(k-j-1,k)= vecteurCirculantD_(ndyad-j-1);
        }
    }
}


// resolution par gradient conjugue preconditionne d'un systeme 
// dont la matrice est  la matrice de Toeplitz
// (le vecteur b est le second membre)
// renvoie le nombre d'iterations (-1 si non converge
int Toeplitz::gradientToeplitzPreconTChang(TVector<complex<double> >& b) const
{
  int k;
  if (dimTop_ != b.Size())
    {
      throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n");
    }
  if (!hermitian_) 
    {
      throw PException(" TOEPLITZMATRIX : conjugate gradient for hermitian matrix only\n");
    }
  int ndyad = vecteurCirculant_.Size();
  // preconditionneur (en Fourier)
  TVector<complex<double> > TFourierC;
  fabricationTChangPreconHerm(TFourierC);


  TVector<complex<double> > a(dimTop_);
  TVector<complex<double> > r(dimTop_);
  TVector<complex<double> > z(dimTop_);
  TVector<complex<double> > q(ndyad);

  a =complex<double>(0.,0.);
  r = b;
  inverseSystemeCirculantFourier(TFourierC, r, z);
  q =complex<double>(0.,0.);
  for (k=0; k<dimTop_; k++) q(k) = z(k);

  complex<double> norm2zRk = prodScalaire(dimTop_,z,r);  
  complex<double> norm2zR0 = norm2zRk;  
  int niter =0;
  double critere = 1.e+8;
    do 
    {
      int k;
      TVector<complex<double> > Aq;
      produitParVecFourier(q, Aq);

      // matrice hermitienne
      double  norme2Aq = prodScalaire(dimTop_,Aq,q).real();
      complex<double> alpha = norm2zRk/norme2Aq;
      for (k=0; k<dimTop_; k++) a(k) += alpha*q(k);
      for (k=0; k<dimTop_; k++) r(k) -= alpha*Aq(k);

      inverseSystemeCirculantFourier(TFourierC, r, z);

      complex<double> norm2zRkp1 = prodScalaire(dimTop_,z,r); 
      complex<double> beta = norm2zRkp1/norm2zRk;
      for (k=0; k<dimTop_; k++) q(k) = z(k) + beta*q(k);

      norm2zRk = norm2zRkp1;
      complex<double> rapport = norm2zRk/norm2zR0;
      critere = fabs(rapport.real())+fabs(rapport.imag());
      //   cout << " norme residu= " << norm2zRk << " critere arret " << critere << endl;
      niter++;
        }

    while ( critere > 1.e-12 && niter < 2*dimTop_);
    b = a;
    if (niter == 2*dimTop_-1) niter = -1;
    return niter;
}

// extension de la matrice Toeplitz a une matrice circulante, avec zeros
// additionnels : pour pouvoir utiliser les FFT
void Toeplitz::extensionACirculanteDyadique(const TVector<complex<double> >& firstCol, const TVector<complex<double> >& firstRow)
{
  // la matrice de Toeplitz doit etre etendue a (2N-1)x(2N_1) pour 
  // la rendre circulante

  int k;
  int N = dimTop_;
  int nEtendu = 2*N-1;
  // on va completer avec des zeros, pour avoir une puissance de 2 (pour fft)
  int ndyad=1;
  while ( nEtendu >  ndyad)  
    {
      ndyad *=2;
    }
  //  cout << " extensionToeplitz: taille etendue de la matrice " << nEtendu << endl;
  //  cout << " puissance de 2 trouvee " << ndyad << endl;


  // extension de la matrice
  vecteurCirculant_.ReSize(ndyad);
  vecteurCirculant_ = complex<double>(0.,0.);
  vecteurCirculant_(0) = firstCol(0);
  for (k=1; k< N; k++)
    {
      vecteurCirculant_(k) = firstCol(k);
      vecteurCirculant_(ndyad-k) = firstRow(k);
    }
  transformeeFourier(vecteurCirculant_, CirculanteFourier_);
}


// extension de la matrice Toeplitz a une matrice circulante, avec zeros
// additionnels : pour pouvoir utiliser les FFT
void Toeplitz::extensionACirculanteDyadique(const TVector<complex<double> >& firstCol)
{
  // la matrice de Toeplitz doit etre etendue a (2N-1)x(2N_1) pour 
  // la rendre circulante

  int k;
  int N = dimTop_;
  int nEtendu = 2*N-1;
  // on va completer avec des zeros, pour avoir une puissance de 2 (pour fft)
  int ndyad=1;
  while ( nEtendu >  ndyad)  
    {
      ndyad *=2;
    }
  //  cout << " extensionToeplitz: taille etendue de la matrice " << nEtendu << endl;
  //  cout << " puissance de 2 trouvee " << ndyad << endl;


  // extension de la matrice
  vecteurCirculant_.ReSize(ndyad);
  vecteurCirculant_ = complex<double>(0.,0.);
  vecteurCirculant_(0) = firstCol(0);
  for (k=1; k< N; k++)
    {
      complex<double> aux = firstCol(k);
      vecteurCirculant_(k) = aux;
      vecteurCirculant_(ndyad-k) = conj(aux);
    }
  transformeeFourier(vecteurCirculant_, CirculanteFourier_);
}


// extension de la matrice Toeplitz a une matrice circulante, avec zeros
// additionnels : pour pouvoir utiliser les FFT
void Toeplitz::extensionACirculanteDyadique(const TVector<double>& firstCol, const TVector<double>& firstRow)
{
  // la matrice de Toeplitz doit etre etendue a (2N-1)x(2N_1) pour 
  // la rendre circulante

  int k;
  int N = dimTop_;
  int nEtendu = 2*N-1;
  // on va completer avec des zeros, pour avoir une puissance de 2 (pour fft)
  int ndyad=1;
  while ( nEtendu >  ndyad)  
    {
      ndyad *=2;
    }
  //  cout << " extensionToeplitz: taille etendue de la matrice " << nEtendu << endl;
  //  cout << " puissance de 2 trouvee " << ndyad << endl;


  // extension de la matrice
  vecteurCirculantD_.ReSize(ndyad);
  vecteurCirculantD_ = 0.;

  vecteurCirculantD_(0) = firstCol(0);
  for (k=1; k< N; k++)
    {
      vecteurCirculantD_(k) = firstCol(k);
      vecteurCirculantD_(ndyad-k) = firstRow(k);
    }
  transformeeFourier(vecteurCirculantD_, CirculanteFourier_);
}



// extension de la matrice Toeplitz a une matrice circulante, avec zeros
// additionnels : pour pouvoir utiliser les FFT
void Toeplitz::extensionACirculanteDyadique(const TVector<double>& firstCol)
{
  // la matrice de Toeplitz doit etre etendue a (2N-1)x(2N_1) pour 
  // la rendre circulante


  int k;
  int N = dimTop_;
  int nEtendu = 2*N-1;
  // on va completer avec des zeros, pour avoir une puissance de 2 (pour fft)
  int ndyad=1;
  while ( nEtendu >  ndyad)  
    {
      ndyad *=2;
    }
  //  cout << " extensionToeplitz: taille etendue de la matrice " << nEtendu << endl;
  //  cout << " puissance de 2 trouvee " << ndyad << endl;


  // extension de la matrice
  vecteurCirculantD_.ReSize(ndyad);
  vecteurCirculantD_ = 0.;
  vecteurCirculantD_(0) = firstCol(0);
  for (k=1; k< N; k++)
    {
      double aux = firstCol(k);
      vecteurCirculantD_(k) = aux;
      vecteurCirculantD_(ndyad-k) = aux;
    }
  transformeeFourier(vecteurCirculantD_, CirculanteFourier_);
}


// produit de la matrice de toeplitz par un vecteur (par transformee de 
// fourier) resultat dans Tq
void Toeplitz::produitParVecFourier(const TVector<complex<double> >& q, TVector<complex<double> >& Tq) const
{
  int k;
  int n = q.Size();
  if (n != vecteurCirculant_.Size() )
    {
      cout << " produitToepVecFourier: dimensions matrice vecteur incompatibles " << endl;
      return;
    }
  TVector<complex<double> > qFourier;
  transformeeFourier(q, qFourier);


    TVector<complex<double> > produitFourier(n);
    for (k=0; k<qFourier.Size(); k++)
     {
      produitFourier(k) = CirculanteFourier_(k)*qFourier(k);
     }
     transformeeInverseFourier(produitFourier,Tq); 
}




// produit de la matrice de toeplitz par un vecteur (par transformee de 
// fourier) resultat dans Tq
void Toeplitz::produitParVecFourier(const TVector<double>& q, TVector<double>& Tq) const
{
  int k;
  int n = q.Size();
  if (n != vecteurCirculantD_.Size() )
    {
      cout << " produitParVecFourier: dimensions matrice vecteur incompatibles " << "vecteur  " << n << " matrice "<< vecteurCirculantD_.Size()  << endl;
      return;
    }
  // cout << " le vecteur " << endl;
  TVector<complex<double> > qFourier;
  transformeeFourier(q, qFourier);
  int  M = qFourier.Size();  
  TVector<complex<double> > produitFourier(M);
  for (k=0; k< M; k++)
    { 
      produitFourier(k) = CirculanteFourier_(k)*qFourier(k);
    }
  transformeeInverseFourier(produitFourier,Tq); 
}
// produit de la transposee de matrice de toeplitz par un vecteur 
// (par transformee de fourier) resultat dans Tq
void Toeplitz::produitTransposeeParVecFourier(const TVector<double>& q, TVector<double>& Tq) const
{
  int k;
  int n = q.Size();
  if (n != vecteurCirculantD_.Size() )
    {
      cout << " produitTransposeeParVecFourier: dimensions matrice vecteur incompatibles " << "vecteur  " << n << " matrice "<< vecteurCirculantD_.Size()  << endl;
      return;
    }
  // cout << " le vecteur " << endl;
  TVector<complex<double> > qFourier;
  transformeeFourier(q, qFourier);
  int  M = qFourier.Size();  
  TVector<complex<double> > produitFourier(M);
  for (k=0; k< M; k++)
    { 
      produitFourier(k) = CirculanteTransposeeFourier_(k)*qFourier(k);
    }
  transformeeInverseFourier(produitFourier,Tq); 
}


void Toeplitz::fabricationTChangPreconHerm(TVector<complex<double> >& TFourierC) const
{
  int k;
  TVector<complex<double> > C(dimTop_);
  C(0) = vecteurCirculant_(0);
  for (k=1; k<dimTop_; k++)
    {
      C(k) = (double(dimTop_-k)*vecteurCirculant_(k) + double(k)*conj(vecteurCirculant_(dimTop_-k)))/double(dimTop_);
    }
  transformeeFourier(C, TFourierC);
  // verifier qu'il n'y a pas de terme nul (pour pouvoir inverser)
  int index =0;
  for (k=0; k<dimTop_; k++)
    {
      if ( fabs(TFourierC(k).real())== 0. &&  fabs(TFourierC(k).imag()) == 0.)
        {
          index = 1;
          break;
        }
    }
  if (index == 1) 
    {
      cout << " fabricationTChangPreconHerm  ; ATTENTION : un terme de la TF du preconditionneur est nul! " << endl;
    }

}
void Toeplitz::inverseSystemeCirculantFourier(const TVector<complex<double> >& TFourierC, const TVector<complex<double> >& secondMembre, TVector<complex<double> >& resul) const
{
  int k;
  int n = TFourierC.Size();
  if (n != secondMembre.Size() )
    {
      cout << " inverseSystemeCirculantFourier: dimensions matrice vecteur incompatibles " << endl;
      return;
    }
  TVector<complex<double> > SmFourier;
  transformeeFourier(secondMembre, SmFourier);
  TVector<complex<double> > quotientFourier(n);
  for (k=0; k<n; k++)
    {
      quotientFourier(k) = SmFourier(k)/TFourierC(k);
    }
  transformeeInverseFourier(quotientFourier,resul); 
}