Index: /trunk/SophyaLib/TArray/toeplitzMatrix.cc
===================================================================
--- /trunk/SophyaLib/TArray/toeplitzMatrix.cc	(revision 1941)
+++ /trunk/SophyaLib/TArray/toeplitzMatrix.cc	(revision 1941)
@@ -0,0 +1,787 @@
+
+#include "toeplitzMatrix.h"
+#include "timing.h"   
+Toeplitz::Toeplitz() : hermitian_(false)
+  {
+    fftIntfPtr_=new FFTPackServer;
+    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;
+    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); 
+}
Index: /trunk/SophyaLib/TArray/toeplitzMatrix.h
===================================================================
--- /trunk/SophyaLib/TArray/toeplitzMatrix.h	(revision 1941)
+++ /trunk/SophyaLib/TArray/toeplitzMatrix.h	(revision 1941)
@@ -0,0 +1,175 @@
+#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.h>
+
+#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)
+
+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;
+}
+
+
+inline void transformeeFourier(const TVector<complex<double> >& v, TVector< complex<double> >& tv) const
+{
+  fftIntfPtr_-> FFTForward(v, tv);
+} 
+inline void transformeeFourier(const TVector<double>& v, TVector< complex<double> >& tv) const
+{
+  fftIntfPtr_-> FFTForward(v, tv);
+} 
+
+
+inline void transformeeInverseFourier(const TVector<complex<double> >& tv, TVector<complex<double> >& v) const
+{
+  int n= tv.Size();
+  fftIntfPtr_-> FFTBackward(tv, v);
+  v/=n;
+} 
+inline void transformeeInverseFourier(const TVector<complex<double> >& tv, TVector<double>& v) const
+{
+  int n= tv.Size();
+  fftIntfPtr_-> FFTBackward(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
+
+
+#endif
