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

Last change on this file since 2806 was 2806, checked in by ansari, 20 years ago

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