source: Sophya/trunk/SophyaLib/TArray/toeplitzMatrix.cc@ 2153

Last change on this file since 2153 was 1941, checked in by lemeur, 24 years ago

matrice de Toeplitz
File size: 21.6 KB
Line 
1
2#include "toeplitzMatrix.h"
3#include "timing.h"
4Toeplitz::Toeplitz() : hermitian_(false)
5 {
6 fftIntfPtr_=new FFTPackServer;
7 fftIntfPtr_->setNormalize(false);
8 }
9
10// matrice de Toplitz hermitienne. La premiere ligne est conjuguee de la
11// premeire colonne donnee (vecteur firstCol)
12Toeplitz::Toeplitz(const TVector<complex<double> >& firstCol) : hermitian_(false)
13 {
14 fftIntfPtr_=new FFTPackServer;
15 fftIntfPtr_->setNormalize(false);
16 setMatrix(firstCol);
17 }
18
19
20Toeplitz::~Toeplitz(){ if (fftIntfPtr_!=NULL) delete fftIntfPtr_;};
21
22
23// initialise matrice de Toeplitz hermitienne
24void Toeplitz::setMatrix(const TVector<complex<double> >& firstCol)
25 {
26 hermitian_ = true;
27 dimTop_ = firstCol.Size();
28 extensionACirculanteDyadique(firstCol);
29 // transformeeFourier(vecteurCirculant_, CirculanteFourier_);
30 }
31
32// intialise la matrice de Toeplitz generale complexe
33void Toeplitz::setMatrix(const TVector<complex<double> >& firstCol, const TVector<complex<double> >& firstRow)
34 {
35 hermitian_ = false;
36 dimTop_ = firstCol.Size();
37 if (dimTop_ != firstRow.Size())
38 {
39 cout << " Toeplitz::setMatrix : nb col different de nbLignes " << endl;
40 }
41 extensionACirculanteDyadique(firstCol, firstRow);
42 // transformeeFourier(vecteurCirculant_, CirculanteFourier_);
43 }
44
45
46// intialise la matrice de Toeplitz generale reelle
47void Toeplitz::setMatrix(const TVector<double>& firstCol, const TVector<double>& firstRow)
48 {
49 hermitian_ = false;
50 dimTop_ = firstCol.Size();
51 if (dimTop_ != firstRow.Size())
52 {
53 cout << " Toeplitz::setMatrix : nb col different de nbLignes " << endl;
54 }
55 extensionACirculanteDyadique(firstCol, firstRow);
56 // transformeeFourier(vecteurCirculantD_, CirculanteFourier_);
57 }
58
59// intialise la matrice de Toeplitz generale reelle symetrique
60void Toeplitz::setMatrix(const TVector<double>& firstCol)
61 {
62 hermitian_ = true;
63 dimTop_ = firstCol.Size();
64 extensionACirculanteDyadique(firstCol);
65 // transformeeFourier(vecteurCirculantD_, CirculanteFourier_);
66 }
67
68
69// resolution par gradient conjugue d'un systeme dont la matrice est
70// la matrice de Toeplitz
71// (le vecteur b est le second membre)
72// renvoie le nombre d'iterations (-1 si non convergence)
73int Toeplitz::gradientToeplitz(TVector<complex<double> >& b) const
74{
75 // PrtTim(" entree gradient " );
76 int k;
77 if (dimTop_ != b.Size())
78 {
79 throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n");
80 }
81 if (!hermitian_)
82 {
83 throw PException(" TOEPLITZMATRIX : conjugate gradient for hermitian matrix only\n");
84 }
85 int ndyad = vecteurCirculant_.Size();
86
87 TVector<complex<double> > a(dimTop_);
88 TVector<complex<double> > r(dimTop_);
89 TVector<complex<double> > q(ndyad);
90
91 a =complex<double>(0.,0.);
92
93 r = b;
94 q =complex<double>(0.,0.);
95 // int dimTopM1 = dimTop_-1;
96 for (k=0; k<dimTop_; k++) q(k) = r(k);
97 //q(Range(0,dimTopM1)) = r;
98 //q.Print(cout, q.Size());
99
100 double norm2Rk = prodScalaire(dimTop_,r,r).real();
101 double norm2R0 = norm2Rk;
102 int niter =0;
103 do
104 {
105 int k;
106 TVector<complex<double> > Aq;
107 produitParVecFourier(q, Aq);
108 //
109 double norme2Aq = prodScalaire(dimTop_,Aq,q).real();
110 complex<double> alpha = prodScalaire(dimTop_,r,r)/norme2Aq;
111 for (k=0; k<dimTop_; k++) a(k) += alpha*q(k);
112 for (k=0; k<dimTop_; k++) r(k) -= alpha*Aq(k);
113 double norm2Rkp1 = prodScalaire(dimTop_,r,r).real();
114 double beta = norm2Rkp1/norm2Rk;
115 for (k=0; k<dimTop_; k++) q(k) = r(k) + beta*q(k);
116
117 norm2Rk = norm2Rkp1;
118
119 // cout << " iteration " << niter+1 << " norme residu= " << norm2Rk << " critere arret " << norm2Rk/norm2R0 << endl;
120 niter++;
121 }
122
123 while (norm2Rk/norm2R0 > 1.e-12 && niter < 2*dimTop_);
124
125 // PrtTim(" fin du gradient ");
126
127 // verification a supprimer
128 //TMatrix<complex<double> > ttest(dimTop_, dimTop_);
129 //expliciteToeplitz(ttest);
130 //TVector<complex<double> > btest;
131 //ttest.Show();
132 //cout << " long vec circ " << dimTop_ << " long a " << a.Size() << endl;
133 //btest = ttest*a;
134 //for (int k=0; k< b.Size(); k++)
135 // {
136 // cout << " verif gradient donnee " << b(k) << " retrouve " << btest(k) << endl;
137 // }
138 //
139 b = a;
140
141 if (niter == 2*dimTop_-1) niter = -1;
142 return niter;
143}
144
145
146// resolution par gradient conjugue d'un systeme dont la matrice est
147// la matrice de Toeplitz
148// (le vecteur b est le second membre)
149// renvoie le nombre d'iterations (-1 si non converge
150int Toeplitz::gradientToeplitz(TVector<double>& b) const
151{
152 int k;
153 if (dimTop_ != b.Size())
154 {
155 throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n");
156 }
157 if (!hermitian_)
158 {
159 throw PException(" TOEPLITZMATRIX : conjugate gradient for hermitian matrix only\n");
160 }
161
162 int ndyad = vecteurCirculantD_.Size();
163
164 TVector<double> a(dimTop_);
165 TVector<double> r(dimTop_);
166 TVector<double> q(ndyad);
167
168 a =0.;
169
170 r = b;
171 q =0.;
172 for (k=0; k<dimTop_; k++) q(k) = r(k);
173
174
175 double norm2Rk = prodScalaire(dimTop_,r,r);
176 double norm2R0 = norm2Rk;
177 int niter =0;
178 do
179 {
180 int k;
181 TVector<double> Aq;
182 produitParVecFourier(q, Aq);
183
184 //
185 double norme2Aq = prodScalaire(dimTop_,Aq,q);
186 double alpha = prodScalaire(dimTop_,r,r)/norme2Aq;
187 for (k=0; k<dimTop_; k++)
188 {
189 a(k) += alpha*q(k);
190 r(k) -= alpha*Aq(k);
191 }
192 double norm2Rkp1 = prodScalaire(dimTop_,r,r);
193 double beta = norm2Rkp1/norm2Rk;
194 for (k=0; k<dimTop_; k++) q(k) = r(k) + beta*q(k);
195
196 norm2Rk = norm2Rkp1;
197
198 // cout << " iteration " << niter+1 << " norme residu= " << norm2Rk << " critere arret " << norm2Rk/norm2R0 << endl;
199 niter++;
200 }
201
202 while (norm2Rk/norm2R0 > 1.e-12 && niter < 2*dimTop_);
203 // while (norm2Rk/norm2R0 > 1.e-12 && niter < 1);
204
205
206
207 // verification a supprimer
208 //TMatrix<double> ttest(dimTop_, dimTop_);
209 //expliciteToeplitz(ttest);
210 // TVector<double> btest;
211 // ttest.Show();
212 //cout << " long vec circ " << dimTop_ << " long a " << a.Size() << endl;
213 // btest = ttest*a;
214 // for (int k=0; k< b.Size(); k++)
215 // {
216 // cout << " verif gradient donnee " << b(k) << " retrouve " << btest(k) << endl;
217 // }
218 //
219 b = a;
220
221 if (niter == 2*dimTop_-1) niter = -1;
222 return niter;
223}
224// resolution par Gradient Conjugue Generalise d'un systeme dont
225// la matrice est la matrice de Toeplitz
226// (le vecteur b est le second membre)
227// renvoie le nombre d'iterations (-1 si non converge
228int Toeplitz::CGSToeplitz(TVector<double>& b) const
229{
230 int k;
231 if (dimTop_ != b.Size())
232 {
233 throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n");
234 }
235
236 int ndyad = vecteurCirculantD_.Size();
237
238 TVector<double> x(dimTop_);
239 TVector<double> r(dimTop_);
240 TVector<double> r0(dimTop_);
241 TVector<double> q(ndyad);
242 TVector<double> p(ndyad);
243
244 // initialisations
245 x =0.;
246 r0 = b;
247 r = r0;
248 q =0.;
249 p =0.;
250 for (k=0; k<dimTop_; k++)
251 {
252 q(k) = r(k);
253 p(k) = r(k);
254 }
255 int niter =0;
256 double r0rk = prodScalaire(dimTop_,r0,r);
257 double r0r0 = r0rk;
258 double rkrk = r0r0;
259 double r0rkp1 = r0rk;
260 do
261 {
262 int k;
263 TVector<double> Aq;
264 produitParVecFourier(q, Aq);
265 TVector<double> Ap;
266 produitParVecFourier(p, Ap);
267
268 TVector<double> pMalphaAq(ndyad);
269
270
271 //
272
273 double r0Aqk = prodScalaire(dimTop_,r0,Aq);
274 double alpha = r0rk/r0Aqk;
275 pMalphaAq = p-alpha*Aq;
276 TVector<double> aux(ndyad);
277 aux = p + pMalphaAq;
278 TVector<double> Aaux;
279 produitParVecFourier(aux, Aaux);
280
281 for (k=0; k<dimTop_; k++)
282 {
283 x(k) += alpha*aux(k);
284 r(k) -= alpha*Aaux(k);
285 }
286 r0rkp1 = prodScalaire(dimTop_,r0,r);
287 double beta = r0rkp1/r0rk;
288 for (k=0; k<dimTop_; k++) p(k) = r(k) + beta*pMalphaAq(k);
289 for (k=0; k<dimTop_; k++)
290 {
291 q(k) = p(k) + beta*(beta*q(k) + pMalphaAq(k));
292 }
293 r0rk = r0rkp1;
294 rkrk = prodScalaire(dimTop_,r,r);
295 // cout << " iteration " << niter+1 << " rkrk= " << rkrk << " critere arret " << rkrk/r0r0 << endl;
296 niter++;
297 }
298
299 while (rkrk/r0r0 > 1.e-12 && niter < 2*dimTop_);
300 // while ( niter < 10*dimTop_);
301
302 // verification a supprimer
303 //TMatrix<double> ttest(dimTop_, dimTop_);
304 //expliciteToeplitz(ttest);
305 //TVector<double> btest;
306 // ttest.Show();
307 //cout << " long vec circ " << dimTop_ << " long a " << a.Size() << endl;
308 //btest = ttest*x;
309 // for (int k=0; k< b.Size(); k++)
310 // {
311 // cout << " verif gradient donnee " << b(k) << " retrouve " << btest(k) << endl;
312 // }
313 //
314 b = x;
315
316 if (niter == 2*dimTop_) niter = -1;
317 return niter;
318}
319// resolution par Double Gradient Conjugue d'un systeme dont
320// la matrice est la matrice de Toeplitz
321// (le vecteur b est le second membre)
322// renvoie le nombre d'iterations (-1 si non converge
323int Toeplitz::DCGToeplitz(TVector<double>& b)
324{
325 int k;
326 if (dimTop_ != b.Size())
327 {
328 throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n");
329 }
330
331 // on a besoin des produit (matrice trasposee)Xvectteur
332 // on initialise eventuellement la TF de Toeplitz transposee
333 if (CirculanteTransposeeFourier_.Size() == 0) initTFTransposeeToeplitzReelle();
334
335
336 int ndyad = vecteurCirculantD_.Size();
337
338 TVector<double> x(dimTop_);
339 TVector<double> r1(dimTop_);
340 TVector<double> r2(dimTop_);
341 TVector<double> d1(ndyad);
342 TVector<double> d2(ndyad);
343
344 // initialisations
345 x =0.;
346 r1 = b;
347 r2 = r1;
348 d1 =0.;
349 d2 =0.;
350 for (k=0; k<dimTop_; k++)
351 {
352 d1(k) = r1(k);
353 d2(k) = r2(k);
354 }
355 int niter =0;
356 double r1r2k = prodScalaire(dimTop_,r1,r2);
357 double r1r2ini = r1r2k;
358 do
359 {
360 int k;
361 TVector<double> Ad1;
362 produitParVecFourier(d1, Ad1);
363 TVector<double> tAd2;
364 produitTransposeeParVecFourier(d2, tAd2);
365
366 //
367 double d2Ad1 = prodScalaire(dimTop_,d2,Ad1);
368 double alpha = r1r2k/d2Ad1;
369
370 for (k=0; k<dimTop_; k++)
371 {
372 x(k) += alpha*d1(k);
373 r1(k) -= alpha*Ad1(k);
374 r2(k) -= alpha*tAd2(k);
375 }
376 prodScalaire(dimTop_,r1,r2);
377 double beta = prodScalaire(dimTop_,r1,r2)/r1r2k;
378 for (k=0; k<dimTop_; k++)
379 {
380 d1(k) = r1(k) + beta*d1(k);
381 d2(k) = r2(k) + beta*d2(k);
382 }
383 r1r2k = prodScalaire(dimTop_,r1,r2);
384 // cout << " iteration " << niter+1 << " r1r2= " << r1r2k << " critere arret " << r1r2k/r1r2ini << endl;
385 niter++;
386 }
387
388 while ( fabs(r1r2k/r1r2ini) > 1.e-12 && niter < 2*dimTop_);
389
390 // verification a supprimer
391 //TMatrix<double> ttest(dimTop_, dimTop_);
392 //expliciteToeplitz(ttest);
393 //TVector<double> btest;
394 // ttest.Show();
395 //cout << " long vec circ " << dimTop_ << " long a " << a.Size() << endl;
396 //btest = ttest*x;
397 //for (int k=0; k< b.Size(); k++)
398 // {
399 // cout << " verif gradient donnee " << b(k) << " retrouve " << btest(k) << endl;
400 // }
401 //
402 b = x;
403
404 if (niter == 2*dimTop_) niter = -1;
405 return niter;
406}
407
408
409void Toeplitz::expliciteCirculante(TMatrix<complex<double> >& m) const
410{
411 int k;
412 int n= vecteurCirculant_.Size();
413 m.ReSize(n,n);
414 int decal=0;
415 for (k=0; k<n; k++)
416 {
417 int j;
418 int courant = 0;
419 for (j=decal; j< n; j++)
420 {
421 m(j, k) = vecteurCirculant_(courant++);
422 }
423 for (j=0; j<decal; j++)
424 {
425 m(j,k)= vecteurCirculant_(courant++);
426 }
427 decal++;
428 }
429}
430void Toeplitz::expliciteToeplitz(TMatrix<complex<double> >& m) const
431{
432 int k;
433 int ndyad = vecteurCirculant_.Size();
434 m.ReSize(dimTop_, dimTop_);
435 for (k=0; k<dimTop_; k++)
436 {
437 int j;
438 for (j=k; j< dimTop_; j++)
439 {
440 m(j, k) = vecteurCirculant_(j-k);
441 }
442 for (j=0; j<k; j++)
443 {
444 m(k-j-1,k)= vecteurCirculant_(ndyad-j-1);
445 }
446 }
447}
448void Toeplitz::expliciteToeplitz(TMatrix<double>& m) const
449{
450 int k;
451 int ndyad = vecteurCirculantD_.Size();
452 m.ReSize(dimTop_, dimTop_);
453 for (k=0; k<dimTop_; k++)
454 {
455 int j;
456 for (j=k; j< dimTop_; j++)
457 {
458 m(j, k) = vecteurCirculantD_(j-k);
459 }
460 for (j=0; j < k; j++)
461 {
462 m(k-j-1,k)= vecteurCirculantD_(ndyad-j-1);
463 }
464 }
465}
466
467
468// resolution par gradient conjugue preconditionne d'un systeme
469// dont la matrice est la matrice de Toeplitz
470// (le vecteur b est le second membre)
471// renvoie le nombre d'iterations (-1 si non converge
472int Toeplitz::gradientToeplitzPreconTChang(TVector<complex<double> >& b) const
473{
474 int k;
475 if (dimTop_ != b.Size())
476 {
477 throw SzMismatchError("TOEPLITZMATRIX : gradientToeplitz : RHS dimensions mismatch \n");
478 }
479 if (!hermitian_)
480 {
481 throw PException(" TOEPLITZMATRIX : conjugate gradient for hermitian matrix only\n");
482 }
483 int ndyad = vecteurCirculant_.Size();
484 // preconditionneur (en Fourier)
485 TVector<complex<double> > TFourierC;
486 fabricationTChangPreconHerm(TFourierC);
487
488
489 TVector<complex<double> > a(dimTop_);
490 TVector<complex<double> > r(dimTop_);
491 TVector<complex<double> > z(dimTop_);
492 TVector<complex<double> > q(ndyad);
493
494 a =complex<double>(0.,0.);
495 r = b;
496 inverseSystemeCirculantFourier(TFourierC, r, z);
497 q =complex<double>(0.,0.);
498 for (k=0; k<dimTop_; k++) q(k) = z(k);
499
500 complex<double> norm2zRk = prodScalaire(dimTop_,z,r);
501 complex<double> norm2zR0 = norm2zRk;
502 int niter =0;
503 double critere = 1.e+8;
504 do
505 {
506 int k;
507 TVector<complex<double> > Aq;
508 produitParVecFourier(q, Aq);
509
510 // matrice hermitienne
511 double norme2Aq = prodScalaire(dimTop_,Aq,q).real();
512 complex<double> alpha = norm2zRk/norme2Aq;
513 for (k=0; k<dimTop_; k++) a(k) += alpha*q(k);
514 for (k=0; k<dimTop_; k++) r(k) -= alpha*Aq(k);
515
516 inverseSystemeCirculantFourier(TFourierC, r, z);
517
518 complex<double> norm2zRkp1 = prodScalaire(dimTop_,z,r);
519 complex<double> beta = norm2zRkp1/norm2zRk;
520 for (k=0; k<dimTop_; k++) q(k) = z(k) + beta*q(k);
521
522 norm2zRk = norm2zRkp1;
523 complex<double> rapport = norm2zRk/norm2zR0;
524 critere = fabs(rapport.real())+fabs(rapport.imag());
525 // cout << " norme residu= " << norm2zRk << " critere arret " << critere << endl;
526 niter++;
527 }
528
529 while ( critere > 1.e-12 && niter < 2*dimTop_);
530 b = a;
531 if (niter == 2*dimTop_-1) niter = -1;
532 return niter;
533}
534
535// extension de la matrice Toeplitz a une matrice circulante, avec zeros
536// additionnels : pour pouvoir utiliser les FFT
537void Toeplitz::extensionACirculanteDyadique(const TVector<complex<double> >& firstCol, const TVector<complex<double> >& firstRow)
538{
539 // la matrice de Toeplitz doit etre etendue a (2N-1)x(2N_1) pour
540 // la rendre circulante
541
542 int k;
543 int N = dimTop_;
544 int nEtendu = 2*N-1;
545 // on va completer avec des zeros, pour avoir une puissance de 2 (pour fft)
546 int ndyad=1;
547 while ( nEtendu > ndyad)
548 {
549 ndyad *=2;
550 }
551 // cout << " extensionToeplitz: taille etendue de la matrice " << nEtendu << endl;
552 // cout << " puissance de 2 trouvee " << ndyad << endl;
553
554
555 // extension de la matrice
556 vecteurCirculant_.ReSize(ndyad);
557 vecteurCirculant_ = complex<double>(0.,0.);
558 vecteurCirculant_(0) = firstCol(0);
559 for (k=1; k< N; k++)
560 {
561 vecteurCirculant_(k) = firstCol(k);
562 vecteurCirculant_(ndyad-k) = firstRow(k);
563 }
564 transformeeFourier(vecteurCirculant_, CirculanteFourier_);
565}
566
567
568// extension de la matrice Toeplitz a une matrice circulante, avec zeros
569// additionnels : pour pouvoir utiliser les FFT
570void Toeplitz::extensionACirculanteDyadique(const TVector<complex<double> >& firstCol)
571{
572 // la matrice de Toeplitz doit etre etendue a (2N-1)x(2N_1) pour
573 // la rendre circulante
574
575 int k;
576 int N = dimTop_;
577 int nEtendu = 2*N-1;
578 // on va completer avec des zeros, pour avoir une puissance de 2 (pour fft)
579 int ndyad=1;
580 while ( nEtendu > ndyad)
581 {
582 ndyad *=2;
583 }
584 // cout << " extensionToeplitz: taille etendue de la matrice " << nEtendu << endl;
585 // cout << " puissance de 2 trouvee " << ndyad << endl;
586
587
588 // extension de la matrice
589 vecteurCirculant_.ReSize(ndyad);
590 vecteurCirculant_ = complex<double>(0.,0.);
591 vecteurCirculant_(0) = firstCol(0);
592 for (k=1; k< N; k++)
593 {
594 complex<double> aux = firstCol(k);
595 vecteurCirculant_(k) = aux;
596 vecteurCirculant_(ndyad-k) = conj(aux);
597 }
598 transformeeFourier(vecteurCirculant_, CirculanteFourier_);
599}
600
601
602// extension de la matrice Toeplitz a une matrice circulante, avec zeros
603// additionnels : pour pouvoir utiliser les FFT
604void Toeplitz::extensionACirculanteDyadique(const TVector<double>& firstCol, const TVector<double>& firstRow)
605{
606 // la matrice de Toeplitz doit etre etendue a (2N-1)x(2N_1) pour
607 // la rendre circulante
608
609 int k;
610 int N = dimTop_;
611 int nEtendu = 2*N-1;
612 // on va completer avec des zeros, pour avoir une puissance de 2 (pour fft)
613 int ndyad=1;
614 while ( nEtendu > ndyad)
615 {
616 ndyad *=2;
617 }
618 // cout << " extensionToeplitz: taille etendue de la matrice " << nEtendu << endl;
619 // cout << " puissance de 2 trouvee " << ndyad << endl;
620
621
622 // extension de la matrice
623 vecteurCirculantD_.ReSize(ndyad);
624 vecteurCirculantD_ = 0.;
625
626 vecteurCirculantD_(0) = firstCol(0);
627 for (k=1; k< N; k++)
628 {
629 vecteurCirculantD_(k) = firstCol(k);
630 vecteurCirculantD_(ndyad-k) = firstRow(k);
631 }
632 transformeeFourier(vecteurCirculantD_, CirculanteFourier_);
633}
634
635
636
637// extension de la matrice Toeplitz a une matrice circulante, avec zeros
638// additionnels : pour pouvoir utiliser les FFT
639void Toeplitz::extensionACirculanteDyadique(const TVector<double>& firstCol)
640{
641 // la matrice de Toeplitz doit etre etendue a (2N-1)x(2N_1) pour
642 // la rendre circulante
643
644
645 int k;
646 int N = dimTop_;
647 int nEtendu = 2*N-1;
648 // on va completer avec des zeros, pour avoir une puissance de 2 (pour fft)
649 int ndyad=1;
650 while ( nEtendu > ndyad)
651 {
652 ndyad *=2;
653 }
654 // cout << " extensionToeplitz: taille etendue de la matrice " << nEtendu << endl;
655 // cout << " puissance de 2 trouvee " << ndyad << endl;
656
657
658 // extension de la matrice
659 vecteurCirculantD_.ReSize(ndyad);
660 vecteurCirculantD_ = 0.;
661 vecteurCirculantD_(0) = firstCol(0);
662 for (k=1; k< N; k++)
663 {
664 double aux = firstCol(k);
665 vecteurCirculantD_(k) = aux;
666 vecteurCirculantD_(ndyad-k) = aux;
667 }
668 transformeeFourier(vecteurCirculantD_, CirculanteFourier_);
669}
670
671
672// produit de la matrice de toeplitz par un vecteur (par transformee de
673// fourier) resultat dans Tq
674void Toeplitz::produitParVecFourier(const TVector<complex<double> >& q, TVector<complex<double> >& Tq) const
675{
676 int k;
677 int n = q.Size();
678 if (n != vecteurCirculant_.Size() )
679 {
680 cout << " produitToepVecFourier: dimensions matrice vecteur incompatibles " << endl;
681 return;
682 }
683 TVector<complex<double> > qFourier;
684 transformeeFourier(q, qFourier);
685
686
687 TVector<complex<double> > produitFourier(n);
688 for (k=0; k<qFourier.Size(); k++)
689 {
690 produitFourier(k) = CirculanteFourier_(k)*qFourier(k);
691 }
692 transformeeInverseFourier(produitFourier,Tq);
693}
694
695
696
697
698// produit de la matrice de toeplitz par un vecteur (par transformee de
699// fourier) resultat dans Tq
700void Toeplitz::produitParVecFourier(const TVector<double>& q, TVector<double>& Tq) const
701{
702 int k;
703 int n = q.Size();
704 if (n != vecteurCirculantD_.Size() )
705 {
706 cout << " produitParVecFourier: dimensions matrice vecteur incompatibles " << "vecteur " << n << " matrice "<< vecteurCirculantD_.Size() << endl;
707 return;
708 }
709 // cout << " le vecteur " << endl;
710 TVector<complex<double> > qFourier;
711 transformeeFourier(q, qFourier);
712 int M = qFourier.Size();
713 TVector<complex<double> > produitFourier(M);
714 for (k=0; k< M; k++)
715 {
716 produitFourier(k) = CirculanteFourier_(k)*qFourier(k);
717 }
718 transformeeInverseFourier(produitFourier,Tq);
719}
720// produit de la transposee de matrice de toeplitz par un vecteur
721// (par transformee de fourier) resultat dans Tq
722void Toeplitz::produitTransposeeParVecFourier(const TVector<double>& q, TVector<double>& Tq) const
723{
724 int k;
725 int n = q.Size();
726 if (n != vecteurCirculantD_.Size() )
727 {
728 cout << " produitTransposeeParVecFourier: dimensions matrice vecteur incompatibles " << "vecteur " << n << " matrice "<< vecteurCirculantD_.Size() << endl;
729 return;
730 }
731 // cout << " le vecteur " << endl;
732 TVector<complex<double> > qFourier;
733 transformeeFourier(q, qFourier);
734 int M = qFourier.Size();
735 TVector<complex<double> > produitFourier(M);
736 for (k=0; k< M; k++)
737 {
738 produitFourier(k) = CirculanteTransposeeFourier_(k)*qFourier(k);
739 }
740 transformeeInverseFourier(produitFourier,Tq);
741}
742
743
744void Toeplitz::fabricationTChangPreconHerm(TVector<complex<double> >& TFourierC) const
745{
746 int k;
747 TVector<complex<double> > C(dimTop_);
748 C(0) = vecteurCirculant_(0);
749 for (k=1; k<dimTop_; k++)
750 {
751 C(k) = (double(dimTop_-k)*vecteurCirculant_(k) + double(k)*conj(vecteurCirculant_(dimTop_-k)))/double(dimTop_);
752 }
753 transformeeFourier(C, TFourierC);
754 // verifier qu'il n'y a pas de terme nul (pour pouvoir inverser)
755 int index =0;
756 for (k=0; k<dimTop_; k++)
757 {
758 if ( fabs(TFourierC(k).real())== 0. && fabs(TFourierC(k).imag()) == 0.)
759 {
760 index = 1;
761 break;
762 }
763 }
764 if (index == 1)
765 {
766 cout << " fabricationTChangPreconHerm ; ATTENTION : un terme de la TF du preconditionneur est nul! " << endl;
767 }
768
769}
770void Toeplitz::inverseSystemeCirculantFourier(const TVector<complex<double> >& TFourierC, const TVector<complex<double> >& secondMembre, TVector<complex<double> >& resul) const
771{
772 int k;
773 int n = TFourierC.Size();
774 if (n != secondMembre.Size() )
775 {
776 cout << " inverseSystemeCirculantFourier: dimensions matrice vecteur incompatibles " << endl;
777 return;
778 }
779 TVector<complex<double> > SmFourier;
780 transformeeFourier(secondMembre, SmFourier);
781 TVector<complex<double> > quotientFourier(n);
782 for (k=0; k<n; k++)
783 {
784 quotientFourier(k) = SmFourier(k)/TFourierC(k);
785 }
786 transformeeInverseFourier(quotientFourier,resul);
787}
Note: See TracBrowser for help on using the repository browser.