LORENE
FFT991/cfrcheb_interp.C
1 /*
2  * Copyright (c) 2024 Jerome Novak
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  * Chebyshev transforms (cases R_CHEB, R_CHEBP & R_CHEBI)
27  * using the fft991 routine in the "intermediate" case
28  * where there is no coefficicent transform in theta nor in phi.
29  *
30  * Entree:
31  * -------
32  * int* deg : tableau du nombre effectif de degres de liberte dans chacune
33  * des 3 dimensions: le nombre de points de collocation
34  * en r est nr = deg[2] et doit etre de la forme
35  * nr = 2*p + 1
36  * int* dimf : tableau du nombre d'elements de ff dans chacune des trois
37  * dimensions.
38  * On doit avoir dimf[2] >= deg[2] = nr.
39  *
40  * double* ff : tableau des valeurs de la fonction aux nr points de
41  * de collocation
42  *
43  * x_i = - cos( pi i/(nr-1) ) 0 <= i <= nr-1
44  *
45  * Les valeurs de la fonction doivent etre stokees dans le
46  * tableau ff comme suit
47  * f( x_i ) = ff[ dimf[1]*dimf[2] * j + dimf[2] * k + i ]
48  * ou j et k sont les indices correspondant a phi et theta
49  * respectivement.
50  * L'espace memoire correspondant a ce pointeur doit etre
51  * dimf[0]*dimf[1]*dimf[2] et doit avoir ete alloue avant
52  * l'appel a la routine.
53  *
54  * int* dimc : tableau du nombre d'elements de cf dans chacune des trois
55  * dimensions.
56  * On doit avoir dimc[2] >= deg[2] = nr.
57  *
58  * Sortie:
59  * -------
60  * double* cf : tableau des coefficients c_i de la fonction definis
61  * comme suit (a theta et phi fixes)
62  *
63  * f(x) = som_{i=0}^{nr-1} c_i T_i(x) ,
64  *
65  * ou T_i(x) designe le polynome de Tchebyshev de degre i.
66  * Les coefficients c_i (0 <= i <= nr-1) sont stokes dans
67  * le tableau cf comme suit
68  * c_i = cf[ dimc[1]*dimc[2] * j + dimc[2] * k + i ]
69  * ou j et k sont les indices correspondant a phi et theta
70  * respectivement.
71  * L'espace memoire correspondant au pointeur cf doit etre
72  * dimc[0]*dimc[1]*dimc[2] et doit avoir ete alloue avant
73  * l'appel a la routine.
74  *
75  * NB: Si le pointeur ff est egal a cf, la routine ne travaille que sur un
76  * seul tableau, qui constitue une entree/sortie.
77  */
78 
79 /*
80  * $Id: cfrcheb_interp.C,v 1.1 2024/10/08 17:45:18 j_novak Exp $
81  *
82  *
83  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFT991/cfrcheb_interp.C,v 1.1 2024/10/08 17:45:18 j_novak Exp $
84  *
85  */
86 
87 
88 // headers du C
89 #include <cstdlib>
90 #include <cassert>
91 
92 // Prototypes of F77 subroutines
93 #include "headcpp.h"
94 #include "proto_f77.h"
95 
96 // Prototypage des sous-routines utilisees:
97 namespace Lorene {
98 int* facto_ini(int ) ;
99 double* trigo_ini(int ) ;
100 double* cheb_ini(const int) ;
101 double* chebimp_ini(const int ) ;
102 
103  //*****************************************************************************
104 
105  void cfrcheb_interp(const int* deg, const int* dimf, double* ff, const int* dimc,
106  double* cf)
107 
108  {
109  int i, j, k ;
110 
111 // Dimensions des tableaux ff et cf :
112  int n1f = dimf[0] ;
113  int n2f = dimf[1] ;
114  int n3f = dimf[2] ;
115  int n1c = dimc[0] ;
116  int n2c = dimc[1] ;
117  int n3c = dimc[2] ;
118 
119 // Nombres de degres de liberte en r :
120  int nr = deg[2] ;
121 
122 // Tests de dimension:
123  if (nr > n3f) {
124  cout << "cfrcheb: nr > n3f : nr = " << nr << " , n3f = "
125  << n3f << endl ;
126  abort () ;
127  exit(-1) ;
128  }
129  if (nr > n3c) {
130  cout << "cfrcheb: nr > n3c : nr = " << nr << " , n3c = "
131  << n3c << endl ;
132  abort () ;
133  exit(-1) ;
134  }
135  if (n1f > n1c) {
136  cout << "cfrcheb: n1f > n1c : n1f = " << n1f << " , n1c = "
137  << n1c << endl ;
138  abort () ;
139  exit(-1) ;
140  }
141  if (n2f > n2c) {
142  cout << "cfrcheb: n2f > n2c : n2f = " << n2f << " , n2c = "
143  << n2c << endl ;
144  abort () ;
145  exit(-1) ;
146  }
147 
148 
149 // Nombre de points pour la FFT:
150  int nm1 = nr - 1;
151  int nm1s2 = nm1 / 2;
152 
153 // Recherche des tables pour la FFT:
154  int* facto = facto_ini(nm1) ;
155  double* trigo = trigo_ini(nm1) ;
156 
157 // Recherche de la table des sin(psi) :
158  double* sinp = cheb_ini(nr);
159 
160 // Reallocation eventuelle de l'espace memoire decrivant G et le
161 // tableau de travail t1
162 // (la dimension nm1+2 = nr+1 est exigee par la routine fft991)
163  double* g = (double*)( malloc( (nm1+2)*sizeof(double) ) );
164  double* t1 = (double*)( malloc( (nm1+2)*sizeof(double) ) ) ;
165 
166 // Parametres pour la routine FFT991
167  int jump = 1 ;
168  int inc = 1 ;
169  int lot = 1 ;
170  int isign = - 1 ;
171 
172 // boucle sur phi et theta
173 
174  int n2n3f = n2f * n3f ;
175  int n2n3c = n2c * n3c ;
176 
177  for (j=0; j< n1f; j++) {
178 
179  for (k=0; k<n2f; k++) {
180 
181  int i0 = n2n3f * j + n3f * k ; // indice de depart
182  double* ff0 = ff + i0 ; // tableau des donnees a transformer
183 
184  i0 = n2n3c * j + n3c * k ; // indice de depart
185  double* cf0 = cf + i0 ; // tableau resultat
186 
187 /*
188  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
189  * reliee a x par x = - cos(psi) et F(psi) = f(x(psi)).
190  */
191 
192 // Valeur en psi=0 de la partie antisymetrique de F, F_ :
193  double fmoins0 = 0.5 * ( ff0[0] - ff0[nm1] );
194 
195 // Fonction G(psi) = F+(psi) + F_(psi) sin(psi)
196 //---------------------------------------------
197  for ( i = 1; i < nm1s2 ; i++ ) {
198 // ... indice du pt symetrique de psi par rapport a pi/2:
199  int isym = nm1 - i ;
200 // ... F+(psi)
201  double fp = 0.5 * ( ff0[i] + ff0[isym] ) ;
202 // ... F_(psi) sin(psi)
203  double fms = 0.5 * ( ff0[i] - ff0[isym] ) * sinp[i] ;
204  g[i] = fp + fms ;
205  g[isym] = fp - fms ;
206  }
207 //... cas particuliers:
208  g[0] = 0.5 * ( ff0[0] + ff0[nm1] );
209  g[nm1s2] = ff0[nm1s2];
210 
211 // Developpement de G en series de Fourier par une FFT
212 //----------------------------------------------------
213 
214  F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
215 
216 // Coefficients pairs du developmt. de Tchebyshev de f
217 //----------------------------------------------------
218 // Ces coefficients sont egaux aux coefficients en cosinus du developpement
219 // de G en series de Fourier (le facteur 2 vient de la normalisation
220 // de fft991) :
221 
222  cf0[0] = g[0] ;
223  for (i=2; i<nm1; i += 2) cf0[i] = 2*g[i] ;
224  cf0[nm1] = g[nm1] ;
225 
226 // Coefficients impairs du developmt. de Tchebyshev de f
227 //------------------------------------------------------
228 // 1. Coef. c'_k (recurrence amorcee a partir de zero):
229 // NB: Le -4 en facteur de g[i] est du a la normalisation de fft991
230 // (si fft991 donnait reellement les coef. en sinus, il faudrait le
231 // remplacer par un +2.)
232  cf0[1] = 0 ;
233  double som = 0;
234  for ( i = 3; i < nr; i += 2 ) {
235  cf0[i] = cf0[i-2] -4 * g[i] ;
236  som += cf0[i] ;
237  }
238 
239 // 2. Calcul de c_1 :
240  double c1 = - ( fmoins0 + som ) / nm1s2 ;
241 
242 // 3. Coef. c_k avec k impair:
243  cf0[1] = c1 ;
244  for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
245 
246  } // fin de la boucle sur theta
247  } // fin de la boucle sur phi
248 
249  // Menage
250  free(t1) ;
251  free(g) ;
252  } // end of cfrcheb_interp
253 
254 //*****************************************************************************
255 
256  void cfrchebp_interp(const int* deg, const int* dimf, double* ff, const int* dimc,
257  double* cf)
258 
259  {
260 
261  int i, j, k ;
262 
263 // Dimensions des tableaux ff et cf :
264  int n1f = dimf[0] ;
265  int n2f = dimf[1] ;
266  int n3f = dimf[2] ;
267  int n1c = dimc[0] ;
268  int n2c = dimc[1] ;
269  int n3c = dimc[2] ;
270 
271 // Nombres de degres de liberte en r :
272  int nr = deg[2] ;
273 
274 // Tests de dimension:
275  if (nr > n3f) {
276  cout << "cfrchebp: nr > n3f : nr = " << nr << " , n3f = "
277  << n3f << endl ;
278  abort () ;
279  exit(-1) ;
280  }
281  if (nr > n3c) {
282  cout << "cfrchebp: nr > n3c : nr = " << nr << " , n3c = "
283  << n3c << endl ;
284  abort () ;
285  exit(-1) ;
286  }
287  if (n1f > n1c) {
288  cout << "cfrchebp: n1f > n1c : n1f = " << n1f << " , n1c = "
289  << n1c << endl ;
290  abort () ;
291  exit(-1) ;
292  }
293  if (n2f > n2c) {
294  cout << "cfrchebp: n2f > n2c : n2f = " << n2f << " , n2c = "
295  << n2c << endl ;
296  abort () ;
297  exit(-1) ;
298  }
299 
300 // Nombre de points pour la FFT:
301  int nm1 = nr - 1;
302  int nm1s2 = nm1 / 2;
303 
304 // Recherche des tables pour la FFT:
305  int* facto = facto_ini(nm1) ;
306  double* trigo = trigo_ini(nm1) ;
307 
308 // Recherche de la table des sin(psi) :
309  double* sinp = cheb_ini(nr);
310 
311 // tableau de travail G et t1
312 // (la dimension nm1+2 = nr+1 est exigee par la routine fft991)
313  double* g = (double*)( malloc( (nm1+2)*sizeof(double) ) );
314  double* t1 = (double*)( malloc( (nm1+2)*sizeof(double) ) ) ;
315 
316 // Parametres pour la routine FFT991
317  int jump = 1 ;
318  int inc = 1 ;
319  int lot = 1 ;
320  int isign = - 1 ;
321 
322 // boucle sur phi et theta
323 
324  int n2n3f = n2f * n3f ;
325  int n2n3c = n2c * n3c ;
326 
327  for (j=0; j< n1f; j++) {
328 
329  for (k=0; k<n2f; k++) {
330 
331  int i0 = n2n3f * j + n3f * k ; // indice de depart
332  double* ff0 = ff + i0 ; // tableau des donnees a transformer
333 
334  i0 = n2n3c * j + n3c * k ; // indice de depart
335  double* cf0 = cf + i0 ; // tableau resultat
336 
337 /*
338  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
339  * reliee a x par x = cos(psi/2) et F(psi) = f(x(psi)).
340  */
341 
342 // Valeur en psi=0 de la partie antisymetrique de F, F_ :
343  double fmoins0 = 0.5 * ( ff0[nm1] - ff0[0] );
344 
345 // Fonction G(psi) = F+(psi) + F_(psi) sin(psi)
346 //---------------------------------------------
347  for ( i = 1; i < nm1s2 ; i++ ) {
348 // ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
349  int isym = nm1 - i ;
350 // ... indice (dans le tableau ff0) du point x correspondant a psi
351  int ix = nm1 - i ;
352 // ... indice (dans le tableau ff0) du point x correspondant a sym(psi)
353  int ixsym = nm1 - isym ;
354 
355 // ... F+(psi)
356  double fp = 0.5 * ( ff0[ix] + ff0[ixsym] ) ;
357 // ... F_(psi) sin(psi)
358  double fms = 0.5 * ( ff0[ix] - ff0[ixsym] ) * sinp[i] ;
359  g[i] = fp + fms ;
360  g[isym] = fp - fms ;
361  }
362 //... cas particuliers:
363  g[0] = 0.5 * ( ff0[nm1] + ff0[0] );
364  g[nm1s2] = ff0[nm1s2];
365 
366 // Developpement de G en series de Fourier par une FFT
367 //----------------------------------------------------
368 
369  F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
370 
371 // Coefficients pairs du developmt. de Tchebyshev de f
372 //----------------------------------------------------
373 // Ces coefficients sont egaux aux coefficients en cosinus du developpement
374 // de G en series de Fourier (le facteur 2 vient de la normalisation
375 // de fft991) :
376 
377  cf0[0] = g[0] ;
378  for (i=2; i<nm1; i += 2 ) cf0[i] = 2. * g[i] ;
379  cf0[nm1] = g[nm1] ;
380 
381 // Coefficients impairs du developmt. de Tchebyshev de f
382 //------------------------------------------------------
383 // 1. Coef. c'_k (recurrence amorcee a partir de zero)
384 // Le +4. en facteur de g[i] est du a la normalisation de fft991
385 // (si fft991 donnait reellement les coef. en sinus, il faudrait le
386 // remplacer par un -2.)
387  cf0[1] = 0 ;
388  double som = 0;
389  for ( i = 3; i < nr; i += 2 ) {
390  cf0[i] = cf0[i-2] + 4. * g[i] ;
391  som += cf0[i] ;
392  }
393 
394 // 2. Calcul de c_1 :
395  double c1 = ( fmoins0 - som ) / nm1s2 ;
396 
397 // 3. Coef. c_k avec k impair:
398  cf0[1] = c1 ;
399  for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
400 
401 
402  } // fin de la boucle sur theta
403  } // fin de la boucle sur phi
404 
405 // Menage
406  free (t1) ;
407  free (g) ;
408  }
409 //*****************************************************************************
410 
411 void cfrchebi_interp(const int* deg, const int* dimf, double* ff, const int* dimc,
412  double* cf)
413 
414 {
415 
416 int i, j, k ;
417 
418 // Dimensions des tableaux ff et cf :
419  int n1f = dimf[0] ;
420  int n2f = dimf[1] ;
421  int n3f = dimf[2] ;
422  int n1c = dimc[0] ;
423  int n2c = dimc[1] ;
424  int n3c = dimc[2] ;
425 
426 // Nombres de degres de liberte en theta et r :
427  int nr = deg[2] ;
428 
429 // Tests de dimension:
430  if (nr > n3f) {
431  cout << "cfrchebi: nr > n3f : nr = " << nr << " , n3f = "
432  << n3f << endl ;
433  abort () ;
434  exit(-1) ;
435  }
436  if (nr > n3c) {
437  cout << "cfrchebi: nr > n3c : nr = " << nr << " , n3c = "
438  << n3c << endl ;
439  abort () ;
440  exit(-1) ;
441  }
442  if (n1f > n1c) {
443  cout << "cfrchebi: n1f > n1c : n1f = " << n1f << " , n1c = "
444  << n1c << endl ;
445  abort () ;
446  exit(-1) ;
447  }
448  if (n2f > n2c) {
449  cout << "cfrchebi: n2f > n2c : n2f = " << n2f << " , n2c = "
450  << n2c << endl ;
451  abort () ;
452  exit(-1) ;
453  }
454 
455 // Nombre de points pour la FFT:
456  int nm1 = nr - 1;
457  int nm1s2 = nm1 / 2;
458 
459 // Recherche des tables pour la FFT:
460  int* facto = facto_ini(nm1) ;
461  double* trigo = trigo_ini(nm1) ;
462 
463 // Recherche de la table des sin(psi) :
464  double* sinp = cheb_ini(nr);
465 
466 // Recherche de la table des points de collocations x_k :
467  double* x = chebimp_ini(nr);
468 
469 // tableau de travail G et t1
470 // (la dimension nm1+2 = nr+1 est exigee par la routine fft991)
471  double* g = (double*)( malloc( (nm1+2)*sizeof(double) ) );
472  double* t1 = (double*)( malloc( (nm1+2)*sizeof(double) ) ) ;
473 
474 // Parametres pour la routine FFT991
475  int jump = 1 ;
476  int inc = 1 ;
477  int lot = 1 ;
478  int isign = - 1 ;
479 
480 // boucle sur phi et theta
481 
482  int n2n3f = n2f * n3f ;
483  int n2n3c = n2c * n3c ;
484 
485  for (j=0; j< n1f; j++) {
486 
487  for (k=0; k<n2f; k++) {
488 
489  int i0 = n2n3f * j + n3f * k ; // indice de depart
490  double* ff0 = ff + i0 ; // tableau des donnees a transformer
491 
492  i0 = n2n3c * j + n3c * k ; // indice de depart
493  double* cf0 = cf + i0 ; // tableau resultat
494 
495 // Multiplication de la fonction par x (pour la rendre paire)
496 // NB: dans les commentaires qui suivent, on note h(x) = x f(x).
497 // (Les valeurs de h dans l'espace des configurations sont stokees dans le
498 // tableau cf0).
499  cf0[0] = 0 ;
500  for (i=1; i<nr; i++) cf0[i] = x[i] * ff0[i] ;
501 
502 /*
503  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
504  * reliee a x par x = cos(psi/2) et F(psi) = h(x(psi)).
505  */
506 
507 // Valeur en psi=0 de la partie antisymetrique de F, F_ :
508  double fmoins0 = 0.5 * ( cf0[nm1] - cf0[0] );
509 
510 // Fonction G(psi) = F+(psi) + F_(psi) sin(psi)
511 //---------------------------------------------
512  for ( i = 1; i < nm1s2 ; i++ ) {
513 // ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
514  int isym = nm1 - i ;
515 // ... indice (dans le tableau cf0) du point x correspondant a psi
516  int ix = nm1 - i ;
517 // ... indice (dans le tableau cf0) du point x correspondant a sym(psi)
518  int ixsym = nm1 - isym ;
519 
520 // ... F+(psi)
521  double fp = 0.5 * ( cf0[ix] + cf0[ixsym] ) ;
522 // ... F_(psi) sin(psi)
523  double fms = 0.5 * ( cf0[ix] - cf0[ixsym] ) * sinp[i] ;
524  g[i] = fp + fms ;
525  g[isym] = fp - fms ;
526  }
527 //... cas particuliers:
528  g[0] = 0.5 * ( cf0[nm1] + cf0[0] );
529  g[nm1s2] = cf0[nm1s2];
530 
531 // Developpement de G en series de Fourier par une FFT
532 //----------------------------------------------------
533 
534  F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
535 
536 // Coefficients pairs du developmt. de Tchebyshev de h
537 //----------------------------------------------------
538 // Ces coefficients sont egaux aux coefficients en cosinus du developpement
539 // de G en series de Fourier (le facteur 2. vient de la normalisation
540 // de fft991; si fft991 donnait reellement les coef. en cosinus, il faudrait le
541 // remplacer par un +1.) :
542 
543  cf0[0] = g[0] ;
544  for (i=2; i<nm1; i += 2 ) cf0[i] = 2. * g[i] ;
545  cf0[nm1] = g[nm1] ;
546 
547 // Coefficients impairs du developmt. de Tchebyshev de h
548 //------------------------------------------------------
549 // 1. Coef. c'_k (recurrence amorcee a partir de zero)
550 // Le +4. en facteur de g[i] est du a la normalisation de fft991
551 // (si fft991 donnait reellement les coef. en sinus, il faudrait le
552 // remplacer par un -2.)
553  cf0[1] = 0 ;
554  double som = 0;
555  for ( i = 3; i < nr; i += 2 ) {
556  cf0[i] = cf0[i-2] + 4. * g[i] ;
557  som += cf0[i] ;
558  }
559 
560 // 2. Calcul de c_1 :
561  double c1 = ( fmoins0 - som ) / nm1s2 ;
562 
563 // 3. Coef. c_k avec k impair:
564  cf0[1] = c1 ;
565  for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
566 
567 // Coefficients de f en fonction de ceux de h
568 //-------------------------------------------
569 
570  cf0[0] = 2* cf0[0] ;
571  for (i=1; i<nm1; i++) {
572  cf0[i] = 2 * cf0[i] - cf0[i-1] ;
573  }
574  cf0[nm1] = 0 ;
575 
576 
577  } // fin de la boucle sur theta
578  } // fin de la boucle sur phi
579 
580  // Menage
581  free (t1) ;
582  free (g) ;
583 
584 }
585 
586 } // end of namespace Lorene
587 
Lorene
Lorene prototypes.
Definition: app_hor.h:67