LORENE
FFTW3/citcossincp.C
1 /*
2  * Copyright (c) 1999-2001 Eric Gourgoulhon
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 char citcossincp_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/citcossincp.C,v 1.4 2014/10/13 08:53:20 j_novak Exp $" ;
24 
25 
26 /*
27  * Transformation inverse cos(2*l*theta) ou sin((2*l+1)*theta) (suivant la
28  * parite de l'indice m en phi) sur le deuxieme indice (theta)
29  * d'un tableau 3-D representant une fonction symetrique par rapport
30  * au plan z=0.
31  * Utilise la bibliotheque fftw
32  *
33  * Entree:
34  * -------
35  * int* deg : tableau du nombre effectif de degres de liberte dans chacune
36  * des 3 dimensions: le nombre de points de collocation
37  * en theta est nt = deg[1] et doit etre de la forme
38  * nt = 2*p + 1
39  * int* dimc : tableau du nombre d'elements de cc dans chacune des trois
40  * dimensions.
41  * On doit avoir dimc[1] >= deg[1] = nt.
42  *
43  * double* cf : tableau des coefficients c_l de la fonction definis
44  * comme suit (a r et phi fixes)
45  *
46  * pour m pair:
47  * f(theta) = som_{l=0}^{nt-1} c_l cos( 2 l theta ) .
48  * pour m impair:
49  * f(theta) = som_{l=0}^{nt-2} c_l sin( (2 l+1) theta ) .
50  *
51  * L'espace memoire correspondant a ce
52  * pointeur doit etre dimc[0]*dimc[1]*dimc[2] et doit
53  * avoir ete alloue avant l'appel a la routine.
54  * Le coefficient c_l (0 <= l <= nt-1) doit etre stoke dans
55  * le tableau cf comme suit
56  * c_l = cf[ dimc[1]*dimc[2] * j + k + dimc[2] * l ]
57  * ou j et k sont les indices correspondant a
58  * phi et r respectivement.
59  *
60  * int* dimf : tableau du nombre d'elements de ff dans chacune des trois
61  * dimensions.
62  * On doit avoir dimf[1] >= deg[1] = nt.
63  *
64  * Sortie:
65  * -------
66  * double* ff : tableau des valeurs de la fonction aux nt points de
67  * de collocation
68  *
69  * theta_l = pi/2 l/(nt-1) 0 <= l <= nt-1
70  *
71  * L'espace memoire correspondant a ce
72  * pointeur doit etre dimf[0]*dimf[1]*dimf[2] et doit
73  * avoir ete alloue avant l'appel a la routine.
74  * Les valeurs de la fonction sont stokees
75  * dans le tableau ff comme suit
76  * f( theta_l ) = ff[ dimf[1]*dimf[2] * j + k + dimf[2] * l ]
77  * ou j et k sont les indices correspondant a
78  * phi et r respectivement.
79  *
80  * NB: Si le pointeur cf est egal a ff, la routine ne travaille que sur un
81  * seul tableau, qui constitue une entree/sortie.
82  *
83  */
84 
85 /*
86  * $Id: citcossincp.C,v 1.4 2014/10/13 08:53:20 j_novak Exp $
87  * $Log: citcossincp.C,v $
88  * Revision 1.4 2014/10/13 08:53:20 j_novak
89  * Lorene classes and functions now belong to the namespace Lorene.
90  *
91  * Revision 1.3 2014/10/06 15:18:50 j_novak
92  * Modified #include directives to use c++ syntax.
93  *
94  * Revision 1.2 2013/04/25 15:46:06 j_novak
95  * Added special treatment in the case np = 1, for type_p = NONSYM.
96  *
97  * Revision 1.1 2004/12/21 17:06:03 j_novak
98  * Added all files for using fftw3.
99  *
100  * Revision 1.4 2003/01/31 10:31:23 e_gourgoulhon
101  * Suppressed the directive #include <malloc.h> for malloc is defined
102  * in <stdlib.h>
103  *
104  * Revision 1.3 2002/10/16 14:36:53 j_novak
105  * Reorganization of #include instructions of standard C++, in order to
106  * use experimental version 3 of gcc.
107  *
108  * Revision 1.2 2002/09/09 13:00:40 e_gourgoulhon
109  * Modification of declaration of Fortran 77 prototypes for
110  * a better portability (in particular on IBM AIX systems):
111  * All Fortran subroutine names are now written F77_* and are
112  * defined in the new file C++/Include/proto_f77.h.
113  *
114  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
115  * LORENE
116  *
117  * Revision 2.0 1999/02/22 15:42:27 hyc
118  * *** empty log message ***
119  *
120  *
121  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/citcossincp.C,v 1.4 2014/10/13 08:53:20 j_novak Exp $
122  *
123  */
124 // headers du C
125 #include <cstdlib>
126 #include <fftw3.h>
127 
128 //Lorene prototypes
129 #include "tbl.h"
130 
131 // Prototypage des sous-routines utilisees:
132 namespace Lorene {
133 fftw_plan back_fft(int, Tbl*&) ;
134 double* cheb_ini(const int) ;
135 double* chebimp_ini(const int ) ;
136 //*****************************************************************************
137 
138 void citcossincp(const int* deg, const int* dimc, double* cf, const int* dimf,
139  double* ff)
140 {
141 
142  int i, j, k ;
143 
144  // Dimensions des tableaux ff et cf :
145  int n1f = dimf[0] ;
146  int n2f = dimf[1] ;
147  int n3f = dimf[2] ;
148  int n1c = dimc[0] ;
149  int n2c = dimc[1] ;
150  int n3c = dimc[2] ;
151 
152  // Nombres de degres de liberte en theta :
153  int nt = deg[1] ;
154 
155  // Tests de dimension:
156  if (nt > n2f) {
157  cout << "citcossincp: nt > n2f : nt = " << nt << " , n2f = "
158  << n2f << endl ;
159  abort () ;
160  exit(-1) ;
161  }
162  if (nt > n2c) {
163  cout << "citcossincp: nt > n2c : nt = " << nt << " , n2c = "
164  << n2c << endl ;
165  abort () ;
166  exit(-1) ;
167  }
168  if ( (n1f > 1) && (n1c > n1f) ) {
169  cout << "citcossincp: n1c > n1f : n1c = " << n1c << " , n1f = "
170  << n1f << endl ;
171  abort () ;
172  exit(-1) ;
173  }
174  if (n3c > n3f) {
175  cout << "citcossincp: n3c > n3f : n3c = " << n3c << " , n3f = "
176  << n3f << endl ;
177  abort () ;
178  exit(-1) ;
179  }
180 
181  // Nombre de points pour la FFT:
182  int nm1 = nt - 1;
183  int nm1s2 = nm1 / 2;
184 
185  // Recherche des tables pour la FFT:
186  Tbl* pg = 0x0 ;
187  fftw_plan p = back_fft(nm1, pg) ;
188  Tbl& g = *pg ;
189  double* t1 = new double[nt] ;
190 
191  // Recherche de la table des sin(psi) :
192  double* sinp = cheb_ini(nt);
193 
194  // Recherche de la table des sin( theta_l ) :
195  double* sinth = chebimp_ini(nt);
196 
197  // boucle sur r de 0 a dimf[2])
198 
199  int n2n3f = n2f * n3f ;
200  int n2n3c = n2c * n3c ;
201 
202 /*
203  * Borne de la boucle sur phi:
204  * si n1f = 1, on effectue la boucle une fois seulement.
205  * si n1f > 1, on va jusqu'a j = n1f-2 en sautant j = 1 (les coefficients
206  * j=n1f-1 et j=0 ne sont pas consideres car nuls).
207  */
208  int borne_phi = n1f-1 ;
209  if (n1f == 1) borne_phi = 1 ;
210 
211  //=======================================================================
212  // Cas m pair
213  //=======================================================================
214 
215  j = 0 ;
216 
217  while (j < borne_phi) { //le dernier coef en phi n'est pas considere
218  // (car nul)
219 
220  //-----------------------------------------------------------------------
221  // partie cos(m phi) avec m pair : transformation cos(2 l theta) inverse
222  //-----------------------------------------------------------------------
223 
224  for (k=0; k<n3c; k++) {
225 
226  int i0 = n2n3c * j + k ; // indice de depart
227  double* cf0 = cf + i0 ; // tableau des donnees a transformer
228 
229  i0 = n2n3f * j + k ; // indice de depart
230  double* ff0 = ff + i0 ; // tableau resultat
231 
232  /*
233  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
234  * reliee a x par x = cos(psi/2) et F(psi) = f(x(psi)).
235  */
236 
237  // Calcul des coefficients de Fourier de la fonction
238  // G(psi) = F+(psi) + F_(psi) sin(psi)
239  // en fonction des coefficients en cos(2l theta) de f:
240 
241  // Coefficients impairs de G
242  //--------------------------
243 
244  double c1 = cf0[n3c] ;
245 
246  double som = 0;
247  ff0[n3f] = 0 ;
248  for ( i = 3; i < nt; i += 2 ) {
249  ff0[ n3f*i ] = cf0[ n3c*i ] - c1 ;
250  som += ff0[ n3f*i ] ;
251  }
252 
253  // Valeur en psi=0 de la partie antisymetrique de F, F_ :
254  double fmoins0 = nm1s2 * c1 + som ;
255 
256  // Coef. impairs de G
257  // NB: le facteur 0.25 est du a la normalisation de fftw; si fftw
258  // donnait exactement les coef. des sinus, ce facteur serait -0.5.
259  for ( i = 3; i < nt; i += 2 ) {
260  g.set(nm1-i/2) = 0.25 * ( ff0[ n3f*i ] - ff0[ n3f*(i-2) ] ) ;
261  }
262 
263 
264  // Coefficients pairs de G
265  //------------------------
266  // Ces coefficients sont egaux aux coefficients pairs du developpement de
267  // f.
268  // NB: le facteur 0.5 est du a la normalisation de fftw; si fftw
269  // donnait exactement les coef. des cosinus, ce facteur serait 1.
270 
271  g.set(0) = cf0[0] ;
272  for (i=1; i<nm1s2; i++ ) g.set(i) = 0.5 * cf0[ n3c*2*i ] ;
273  g.set(nm1s2) = cf0[ n3c*nm1 ] ;
274 
275  // Transformation de Fourier inverse de G
276  //---------------------------------------
277 
278  // FFT inverse
279  fftw_execute(p) ;
280 
281  // Valeurs de f deduites de celles de G
282  //-------------------------------------
283 
284  for ( i = 1; i < nm1s2 ; i++ ) {
285  // ... indice du pt symetrique de psi par rapport a pi/2:
286  int isym = nm1 - i ;
287 
288  double fp = 0.5 * ( g(i) + g(isym) ) ;
289  double fm = 0.5 * ( g(i) - g(isym) ) / sinp[i] ;
290  ff0[ n3f*i ] = fp + fm ;
291  ff0[ n3f*isym ] = fp - fm ;
292  }
293 
294  //... cas particuliers:
295  ff0[0] = g(0) + fmoins0 ;
296  ff0[ n3f*nm1 ] = g(0) - fmoins0 ;
297  ff0[ n3f*nm1s2 ] = g(nm1s2) ;
298 
299  } // fin de la boucle sur r
300 
301  //-----------------------------------------------------------------------
302  // partie sin(m phi) avec m pair : transformation cos(2 l theta) inverse
303  //-----------------------------------------------------------------------
304 
305  j++ ;
306 
307  if ( (j != 1) && (j != borne_phi ) ) {
308  // on effectue le calcul seulement dans les cas ou les coef en phi ne sont
309  // pas nuls
310 
311  for (k=0; k<n3c; k++) {
312 
313  int i0 = n2n3c * j + k ; // indice de depart
314  double* cf0 = cf + i0 ; // tableau des donnees a transformer
315 
316  i0 = n2n3f * j + k ; // indice de depart
317  double* ff0 = ff + i0 ; // tableau resultat
318 
319  /*
320  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
321  * reliee a x par x = cos(psi/2) et F(psi) = f(x(psi)).
322  */
323 
324  // Calcul des coefficients de Fourier de la fonction
325  // G(psi) = F+(psi) + F_(psi) sin(psi)
326  // en fonction des coefficients en cos(2l theta) de f:
327 
328  // Coefficients impairs de G
329  //--------------------------
330 
331  double c1 = cf0[n3c] ;
332 
333  double som = 0;
334  ff0[n3f] = 0 ;
335  for ( i = 3; i < nt; i += 2 ) {
336  ff0[ n3f*i ] = cf0[ n3c*i ] - c1 ;
337  som += ff0[ n3f*i ] ;
338  }
339 
340  // Valeur en psi=0 de la partie antisymetrique de F, F_ :
341  double fmoins0 = nm1s2 * c1 + som ;
342 
343  // Coef. impairs de G
344  // NB: le facteur 0.25 est du a la normalisation de fftw; si fftw
345  // donnait exactement les coef. des sinus, ce facteur serait -0.5.
346  for ( i = 3; i < nt; i += 2 ) {
347  g.set(nm1-i/2) = 0.25 * ( ff0[ n3f*i ] - ff0[ n3f*(i-2) ] ) ;
348  }
349 
350 
351  // Coefficients pairs de G
352  //------------------------
353  // Ces coefficients sont egaux aux coefficients pairs du developpement de
354  // f.
355  // NB: le facteur 0.5 est du a la normalisation de fftw; si fftw
356  // donnait exactement les coef. des cosinus, ce facteur serait 1.
357 
358  g.set(0) = cf0[0] ;
359  for (i=1; i<nm1s2; i++ ) g.set(i) = 0.5 * cf0[ n3c*2*i ] ;
360  g.set(nm1s2) = cf0[ n3c*nm1 ] ;
361 
362  // Transformation de Fourier inverse de G
363  //---------------------------------------
364 
365  // FFT inverse
366  fftw_execute(p) ;
367 
368  // Valeurs de f deduites de celles de G
369  //-------------------------------------
370 
371  for ( i = 1; i < nm1s2 ; i++ ) {
372  // ... indice du pt symetrique de psi par rapport a pi/2:
373  int isym = nm1 - i ;
374 
375  double fp = 0.5 * ( g(i) + g(isym) ) ;
376  double fm = 0.5 * ( g(i) - g(isym) ) / sinp[i] ;
377  ff0[ n3f*i ] = fp + fm ;
378  ff0[ n3f*isym ] = fp - fm ;
379  }
380 
381  //... cas particuliers:
382  ff0[0] = g(0) + fmoins0 ;
383  ff0[ n3f*nm1 ] = g(0) - fmoins0 ;
384  ff0[ n3f*nm1s2 ] = g(nm1s2) ;
385 
386  } // fin de la boucle sur r
387 
388  } // fin du cas ou le calcul etait necessaire (i.e. ou les
389  // coef en phi n'etaient pas nuls)
390 
391  // On passe au cas m pair suivant:
392  j+=3 ;
393 
394  } // fin de la boucle sur les cas m pair
395 
396  //##
397  if (n1f<=3) {
398  delete [] t1 ;
399  return ;
400  }
401 
402  //=======================================================================
403  // Cas m impair
404  //=======================================================================
405 
406  j = 2 ;
407 
408  while (j < borne_phi) { //le dernier coef en phi n'est pas considere
409  // (car nul)
410 
411 //--------------------------------------------------------------------------
412 // partie cos(m phi) avec m impair : transformation sin((2 l+1) theta) inv.
413 //--------------------------------------------------------------------------
414 
415  for (k=0; k<n3c; k++) {
416 
417  int i0 = n2n3c * j + k ; // indice de depart
418  double* cf0 = cf + i0 ; // tableau des donnees a transformer
419 
420  i0 = n2n3f * j + k ; // indice de depart
421  double* ff0 = ff + i0 ; // tableau resultat
422 
423  // Calcul des coefficients du developpement en cos(2 l theta)
424  // de la fonction h(theta) := f(theta) sin(theta)
425  // en fonction de ceux de f (le resultat est stoke dans le tableau t1) :
426  t1[0] = .5 * cf0[0] ;
427  for (i=1; i<nm1; i++) {
428  t1[i] = .5 * ( cf0[ n3c*i ] - cf0[ n3c*(i-1) ] ) ;
429  }
430  t1[nm1] = -.5 * cf0[ n3c*(nt-2) ] ;
431 
432  /*
433  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
434  * reliee a theta par psi = 2 theta et F(psi) = h(theta(psi)).
435  */
436 
437  // Calcul des coefficients de Fourier de la fonction
438  // G(psi) = F+(psi) + F_(psi) sin(psi)
439  // en fonction des coefficients en cos(2l theta) de h:
440 
441  // Coefficients impairs de G
442  //--------------------------
443 
444  double c1 = t1[1] ;
445 
446  double som = 0;
447  ff0[n3f] = 0 ;
448  for ( i = 3; i < nt; i += 2 ) {
449  ff0[ n3f*i ] = t1[i] - c1 ;
450  som += ff0[ n3f*i ] ;
451  }
452 
453  // Valeur en psi=0 de la partie antisymetrique de F, F_ :
454  double fmoins0 = nm1s2 * c1 + som ;
455 
456  // Coef. impairs de G
457  // NB: le facteur 0.25 est du a la normalisation de fftw; si fftw
458  // donnait exactement les coef. des sinus, ce facteur serait -0.5.
459  for ( i = 3; i < nt; i += 2 ) {
460  g.set(nm1-i/2) = 0.25 * ( ff0[ n3f*i ] - ff0[ n3f*(i-2) ] ) ;
461  }
462 
463  // Coefficients pairs de G
464  //------------------------
465  // Ces coefficients sont egaux aux coefficients pairs du developpement de
466  // h.
467  // NB: le facteur 0.5 est du a la normalisation de fftw; si fftw
468  // donnait exactement les coef. des cosinus, ce facteur serait 1.
469 
470  g.set(0) = t1[0] ;
471  for (i=1; i<nm1s2; i ++ ) g.set(i) = 0.5 * t1[2*i] ;
472  g.set(nm1s2) = t1[nm1] ;
473 
474  // Transformation de Fourier inverse de G
475  //---------------------------------------
476 
477  // FFT inverse
478  fftw_execute(p) ;
479 
480  // Valeurs de f deduites de celles de G
481  //-------------------------------------
482 
483  for ( i = 1; i < nm1s2 ; i++ ) {
484  // ... indice du pt symetrique de psi par rapport a pi/2:
485  int isym = nm1 - i ;
486 
487  double fp = 0.5 * ( g(i) + g(isym) ) ;
488  double fm = 0.5 * ( g(i) - g(isym) ) / sinp[i] ;
489  ff0[ n3f*i ] = ( fp + fm ) / sinth[i] ;
490  ff0[ n3f*isym ] = ( fp - fm ) / sinth[isym] ;
491  }
492 
493  //... cas particuliers:
494  ff0[0] = 0 ;
495  ff0[ n3f*nm1 ] = g(0) - fmoins0 ;
496  ff0[ n3f*nm1s2 ] = g(nm1s2) / sinth[nm1s2];
497 
498  } // fin de la boucle sur r
499 
500  //--------------------------------------------------------------------------
501  // partie sin(m phi) avec m impair : transformation sin((2 l+1) theta) inv.
502  //--------------------------------------------------------------------------
503 
504  j++ ;
505 
506  if ( j != borne_phi ) {
507  // on effectue le calcul seulement dans les cas ou les coef en phi ne sont
508  // pas nuls
509 
510  for (k=0; k<n3c; k++) {
511 
512  int i0 = n2n3c * j + k ; // indice de depart
513  double* cf0 = cf + i0 ; // tableau des donnees a transformer
514 
515  i0 = n2n3f * j + k ; // indice de depart
516  double* ff0 = ff + i0 ; // tableau resultat
517 
518  // Calcul des coefficients du developpement en cos(2 l theta)
519  // de la fonction h(theta) := f(theta) sin(theta)
520  // en fonction de ceux de f (le resultat est stoke dans le tableau t1) :
521  t1[0] = .5 * cf0[0] ;
522  for (i=1; i<nm1; i++) {
523  t1[i] = .5 * ( cf0[ n3c*i ] - cf0[ n3c*(i-1) ] ) ;
524  }
525  t1[nm1] = -.5 * cf0[ n3c*(nt-2) ] ;
526 
527  /*
528  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
529  * reliee a theta par psi = 2 theta et F(psi) = h(theta(psi)).
530  */
531 
532  // Calcul des coefficients de Fourier de la fonction
533  // G(psi) = F+(psi) + F_(psi) sin(psi)
534  // en fonction des coefficients en cos(2l theta) de h:
535 
536  // Coefficients impairs de G
537  //--------------------------
538 
539  double c1 = t1[1] ;
540 
541  double som = 0;
542  ff0[n3f] = 0 ;
543  for ( i = 3; i < nt; i += 2 ) {
544  ff0[ n3f*i ] = t1[i] - c1 ;
545  som += ff0[ n3f*i ] ;
546  }
547 
548  // Valeur en psi=0 de la partie antisymetrique de F, F_ :
549  double fmoins0 = nm1s2 * c1 + som ;
550 
551  // Coef. impairs de G
552  // NB: le facteur 0.25 est du a la normalisation de fftw; si fftw
553  // donnait exactement les coef. des sinus, ce facteur serait -0.5.
554  for ( i = 3; i < nt; i += 2 ) {
555  g.set(nm1-i/2) = 0.25 * ( ff0[ n3f*i ] - ff0[ n3f*(i-2) ] ) ;
556  }
557 
558  // Coefficients pairs de G
559  //------------------------
560  // Ces coefficients sont egaux aux coefficients pairs du developpement de
561  // h.
562  // NB: le facteur 0.5 est du a la normalisation de fftw; si fftw
563  // donnait exactement les coef. des cosinus, ce facteur serait 1.
564 
565  g.set(0) = t1[0] ;
566  for (i=1; i<nm1s2; i ++ ) g.set(i) = 0.5 * t1[2*i] ;
567  g.set(nm1s2) = t1[nm1] ;
568 
569  // Transformation de Fourier inverse de G
570  //---------------------------------------
571 
572  // FFT inverse
573  fftw_execute(p) ;
574 
575  // Valeurs de f deduites de celles de G
576  //-------------------------------------
577 
578  for ( i = 1; i < nm1s2 ; i++ ) {
579  // ... indice du pt symetrique de psi par rapport a pi/2:
580  int isym = nm1 - i ;
581 
582  double fp = 0.5 * ( g(i) + g(isym) ) ;
583  double fm = 0.5 * ( g(i) - g(isym) ) / sinp[i] ;
584  ff0[ n3f*i ] = ( fp + fm ) / sinth[i] ;
585  ff0[ n3f*isym ] = ( fp - fm ) / sinth[isym] ;
586  }
587 
588  //... cas particuliers:
589  ff0[0] = 0 ;
590  ff0[ n3f*nm1 ] = g(0) - fmoins0 ;
591  ff0[ n3f*nm1s2 ] = g(nm1s2) / sinth[nm1s2];
592 
593 
594  } // fin de la boucle sur r
595 
596  } // fin du cas ou le calcul etait necessaire (i.e. ou les
597  // coef en phi n'etaient pas nuls)
598 
599  // On passe au cas m impair suivant:
600  j+=3 ;
601 
602  } // fin de la boucle sur les cas m impair
603  delete [] t1 ;
604 
605 }
606 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64