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

Last change on this file since 2788 was 2615, checked in by cmv, 21 years ago

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