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