#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 >& 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 >& firstCol) { hermitian_ = true; dimTop_ = firstCol.Size(); extensionACirculanteDyadique(firstCol); // transformeeFourier(vecteurCirculant_, CirculanteFourier_); } // intialise la matrice de Toeplitz generale complexe void Toeplitz::setMatrix(const TVector >& firstCol, const TVector >& 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& firstCol, const TVector& 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& 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 >& 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 > a(dimTop_); TVector > r(dimTop_); TVector > q(ndyad); a =complex(0.,0.); r = b; q =complex(0.,0.); // int dimTopM1 = dimTop_-1; for (k=0; k > Aq; produitParVecFourier(q, Aq); // double norme2Aq = prodScalaire(dimTop_,Aq,q).real(); complex alpha = prodScalaire(dimTop_,r,r)/norme2Aq; for (k=0; k 1.e-12 && niter < 2*dimTop_); // PrtTim(" fin du gradient "); // verification a supprimer //TMatrix > ttest(dimTop_, dimTop_); //expliciteToeplitz(ttest); //TVector > 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& 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 a(dimTop_); TVector r(dimTop_); TVector q(ndyad); a =0.; r = b; q =0.; for (k=0; k Aq; produitParVecFourier(q, Aq); // double norme2Aq = prodScalaire(dimTop_,Aq,q); double alpha = prodScalaire(dimTop_,r,r)/norme2Aq; for (k=0; k 1.e-12 && niter < 2*dimTop_); // while (norm2Rk/norm2R0 > 1.e-12 && niter < 1); // verification a supprimer //TMatrix ttest(dimTop_, dimTop_); //expliciteToeplitz(ttest); // TVector 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& b) const { int k; if (dimTop_ != b.Size()) { throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n"); } int ndyad = vecteurCirculantD_.Size(); TVector x(dimTop_); TVector r(dimTop_); TVector r0(dimTop_); TVector q(ndyad); TVector p(ndyad); // initialisations x =0.; r0 = b; r = r0; q =0.; p =0.; for (k=0; k Aq; produitParVecFourier(q, Aq); TVector Ap; produitParVecFourier(p, Ap); TVector pMalphaAq(ndyad); // double r0Aqk = prodScalaire(dimTop_,r0,Aq); double alpha = r0rk/r0Aqk; pMalphaAq = p-alpha*Aq; TVector aux(ndyad); aux = p + pMalphaAq; TVector Aaux; produitParVecFourier(aux, Aaux); for (k=0; k 1.e-12 && niter < 2*dimTop_); // while ( niter < 10*dimTop_); // verification a supprimer //TMatrix ttest(dimTop_, dimTop_); //expliciteToeplitz(ttest); //TVector 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& 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 x(dimTop_); TVector r1(dimTop_); TVector r2(dimTop_); TVector d1(ndyad); TVector d2(ndyad); // initialisations x =0.; r1 = b; r2 = r1; d1 =0.; d2 =0.; for (k=0; k Ad1; produitParVecFourier(d1, Ad1); TVector tAd2; produitTransposeeParVecFourier(d2, tAd2); // double d2Ad1 = prodScalaire(dimTop_,d2,Ad1); double alpha = r1r2k/d2Ad1; for (k=0; k 1.e-12 && niter < 2*dimTop_); // verification a supprimer //TMatrix ttest(dimTop_, dimTop_); //expliciteToeplitz(ttest); //TVector 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 >& m) const { int k; int n= vecteurCirculant_.Size(); m.ReSize(n,n); int decal=0; for (k=0; k >& m) const { int k; int ndyad = vecteurCirculant_.Size(); m.ReSize(dimTop_, dimTop_); for (k=0; k& m) const { int k; int ndyad = vecteurCirculantD_.Size(); m.ReSize(dimTop_, dimTop_); for (k=0; k >& 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 > TFourierC; fabricationTChangPreconHerm(TFourierC); TVector > a(dimTop_); TVector > r(dimTop_); TVector > z(dimTop_); TVector > q(ndyad); a =complex(0.,0.); r = b; inverseSystemeCirculantFourier(TFourierC, r, z); q =complex(0.,0.); for (k=0; k norm2zRk = prodScalaire(dimTop_,z,r); complex norm2zR0 = norm2zRk; int niter =0; double critere = 1.e+8; do { int k; TVector > Aq; produitParVecFourier(q, Aq); // matrice hermitienne double norme2Aq = prodScalaire(dimTop_,Aq,q).real(); complex alpha = norm2zRk/norme2Aq; for (k=0; k norm2zRkp1 = prodScalaire(dimTop_,z,r); complex beta = norm2zRkp1/norm2zRk; for (k=0; k 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 >& firstCol, const TVector >& 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(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 >& 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(0.,0.); vecteurCirculant_(0) = firstCol(0); for (k=1; k< N; k++) { complex 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& firstCol, const TVector& 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& 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 >& q, TVector >& Tq) const { int k; int n = q.Size(); if (n != vecteurCirculant_.Size() ) { cout << " produitToepVecFourier: dimensions matrice vecteur incompatibles " << endl; return; } TVector > qFourier; transformeeFourier(q, qFourier); TVector > produitFourier(n); for (k=0; k& q, TVector& 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 > qFourier; transformeeFourier(q, qFourier); int M = qFourier.Size(); TVector > 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& q, TVector& 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 > qFourier; transformeeFourier(q, qFourier); int M = qFourier.Size(); TVector > produitFourier(M); for (k=0; k< M; k++) { produitFourier(k) = CirculanteTransposeeFourier_(k)*qFourier(k); } transformeeInverseFourier(produitFourier,Tq); } void Toeplitz::fabricationTChangPreconHerm(TVector >& TFourierC) const { int k; TVector > C(dimTop_); C(0) = vecteurCirculant_(0); for (k=1; k >& TFourierC, const TVector >& secondMembre, TVector >& resul) const { int k; int n = TFourierC.Size(); if (n != secondMembre.Size() ) { cout << " inverseSystemeCirculantFourier: dimensions matrice vecteur incompatibles " << endl; return; } TVector > SmFourier; transformeeFourier(secondMembre, SmFourier); TVector > quotientFourier(n); for (k=0; k