#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 #include #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 >& firstCol); ~Toeplitz(); // toeplitz complexe generale void setMatrix(const TVector >& firstCol, const TVector >& firstRow); // toeplitz hermitienne void setMatrix(const TVector >& firstCol); // toeplitz reelle generale void setMatrix(const TVector& firstCol, const TVector& firstRow); // toeplitz reelle symetrique void setMatrix(const TVector& firstCol); int gradientToeplitz(TVector >& b) const; int gradientToeplitz(TVector& b) const; int gradientToeplitzPreconTChang(TVector >& b) const; int CGSToeplitz(TVector& b) const; int DCGToeplitz(TVector& b) ; private: inline complex prodScalaire(int n, TVector >& v1, TVector >& v2) const { complex produit = complex(0.,0.); for (int k =0; k& v1, TVector& v2) const { double produit = 0.; for (int k =0; k >& v, TVector< complex >& tv) const { fftIntfPtr_-> FFTForward(v, tv); } inline void transformeeFourier(const TVector& v, TVector< complex >& tv) const { fftIntfPtr_-> FFTForward(v, tv); } inline void transformeeInverseFourier(const TVector >& tv, TVector >& v) const { int n= tv.Size(); fftIntfPtr_-> FFTBackward(tv, v); v/=n; } inline void transformeeInverseFourier(const TVector >& tv, TVector& 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 > vecteurCirculantTranspose(ndyad); vecteurCirculantTranspose = complex(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 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 >& firstCol, const TVector >& firstRow); void extensionACirculanteDyadique(const TVector& firstCol, const TVector& firstRow); void extensionACirculanteDyadique(const TVector >& firstCol); void extensionACirculanteDyadique(const TVector& firstCol); // matrice circulante entree par sa T. de Fourier void produitParVecFourier(const TVector >& q, TVector >& Tq) const; void produitParVecFourier(const TVector& q, TVector& Tq) const; void produitTransposeeParVecFourier(const TVector& q, TVector& Tq) const; // preconditionneur pour une Toeplitz HERMITIENNE // LE RESUTAT EST UNE TRANSFORMEE DE FOURIER void fabricationTChangPreconHerm(TVector >& TFourierC) const; // matrice circulante entree par sa T. de Fourier void inverseSystemeCirculantFourier(const TVector >& TFourierC, const TVector >& secondMembre, TVector >& resul) const; void expliciteCirculante(TMatrix >& m) const; void expliciteToeplitz(TMatrix >& m) const; void expliciteToeplitz(TMatrix& m) const; bool hermitian_; int dimTop_; TVector > vecteurCirculant_; TVector vecteurCirculantD_; TVector > CirculanteFourier_; TVector > CirculanteTransposeeFourier_; FFTServerInterface* fftIntfPtr_; }; // fin classe Toeplitz #endif