source: Sophya/trunk/SophyaLib/NTools/toeplitzMatrix.cc@ 2908

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

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