LORENE
FFT991/cipcossin.C
1 /*
2  * Copyright (c) 1999-2002 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 cipcossin_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFT991/cipcossin.C,v 1.4 2014/10/15 12:48:21 j_novak Exp $" ;
24 
25 
26 /*
27  * Transformation de Fourier inverse sur le premier indice d'un tableau 3-D
28  * par le biais de la routine FFT Fortran FFT991
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 phi est np = deg[0] et doit etre de la forme
35  * np = 2^p 3^q 5^r
36  * int* dimc : tableau du nombre d'elements dans chacune des trois
37  * dimensions du tableau cf
38  * On doit avoir dimc[0] >= deg[0] + 2 = np + 2.
39  *
40  * int* dimf : tableau du nombre d'elements dans chacune des trois
41  * dimensions du tableau ff
42  * On doit avoir dimf[0] >= deg[0] = np .
43  *
44  *
45  * Entree / sortie :
46  * -----------------
47  * double* cf : entree: tableau des coefficients de la fonction f;
48  * L'espace memoire correspondant a ce
49  * pointeur doit etre dimc[0]*dimc[1]*dimc[2] et doit
50  * avoir ete alloue avant l'appel a la routine.
51  * La convention de stokage est la suivante:
52  * cf[ dimc[2]*dimc[1]*k + dimc[2]*j + i ] = c_k 0<= k <= np,
53  * ou les indices j et i correspondent respectivement
54  * a theta et a r et ou les c_k sont les coefficients
55  * du developpement de f en series de Fourier:
56  *
57  * f(phi) = som_{l=0}^{np/2} c_{2l} cos( 2 pi/T l phi )
58  * + c_{2l+1} sin( 2 pi/T l phi ),
59  *
60  * ou T est la periode de f.
61  * En particulier cf[1] = 0.
62  * De plus, cf[np+1] n'est pas egal a c_{np+1}
63  * mais a zero.
64  * !!!! CE TABLEAU EST DETRUIT EN SORTIE !!!!!
65  *
66  * Sortie:
67  * -------
68  * double* ff : tableau des valeurs de la fonction aux points de
69  * de collocation
70  *
71  * phi_k = T k/np 0 <= k <= np-1
72  *
73  * suivant la convention de stokage:
74  * ff[ dimf[2]*dimf[1]*k + dimf[2]*j + i ] = f(phi_k) 0 <= k <= np-1,
75  * les indices j et i correspondant respectivement
76  * a theta et a r.
77  * L'espace memoire correspondant a ce
78  * pointeur doit etre dimf[0]*dimf[1]*dimf[2] et doit
79  * avoir ete alloue avant l'appel a la routine.
80  */
81 
82 /*
83  * $Id: cipcossin.C,v 1.4 2014/10/15 12:48:21 j_novak Exp $
84  * $Log: cipcossin.C,v $
85  * Revision 1.4 2014/10/15 12:48:21 j_novak
86  * Corrected namespace declaration.
87  *
88  * Revision 1.3 2014/10/13 08:53:16 j_novak
89  * Lorene classes and functions now belong to the namespace Lorene.
90  *
91  * Revision 1.2 2014/10/06 15:18:45 j_novak
92  * Modified #include directives to use c++ syntax.
93  *
94  * Revision 1.1 2004/12/21 17:06:01 j_novak
95  * Added all files for using fftw3.
96  *
97  * Revision 1.4 2003/01/31 10:31:23 e_gourgoulhon
98  * Suppressed the directive #include <malloc.h> for malloc is defined
99  * in <stdlib.h>
100  *
101  * Revision 1.3 2002/10/16 14:36:53 j_novak
102  * Reorganization of #include instructions of standard C++, in order to
103  * use experimental version 3 of gcc.
104  *
105  * Revision 1.2 2002/09/09 13:00:40 e_gourgoulhon
106  * Modification of declaration of Fortran 77 prototypes for
107  * a better portability (in particular on IBM AIX systems):
108  * All Fortran subroutine names are now written F77_* and are
109  * defined in the new file C++/Include/proto_f77.h.
110  *
111  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
112  * LORENE
113  *
114  * Revision 2.0 1999/02/22 15:43:58 hyc
115  * *** empty log message ***
116  *
117  *
118  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFT991/cipcossin.C,v 1.4 2014/10/15 12:48:21 j_novak Exp $
119  *
120  */
121 
122 
123 // headers du C
124 #include <cassert>
125 #include <cstdlib>
126 
127 // Prototypes of F77 subroutines
128 #include "headcpp.h"
129 #include "proto_f77.h"
130 
131 // Prototypage des sous-routines utilisees:
132 namespace Lorene {
133 int* facto_ini(int ) ;
134 double* trigo_ini(int ) ;
135 //*****************************************************************************
136 
137 void cipcossin(const int* deg, const int* dimc, const int* dimf,
138  double* cf, double* ff)
139 {
140 
141 int i, j, k, index ;
142 
143 // Dimensions des tableaux ff et cf :
144  int n1f = dimf[0] ;
145  int n2f = dimf[1] ;
146  int n3f = dimf[2] ;
147  int n1c = dimc[0] ;
148  int n2c = dimc[1] ;
149  int n3c = dimc[2] ;
150 
151 // Nombres de degres de liberte en phi:
152  int np = deg[0] ;
153 
154 // Tests de dimension:
155  if (np+2 > n1c) {
156  cout << "cipcossin: np+2 > n1c : np = " << np << " , n1c = "
157  << n1c << endl ;
158  abort () ;
159  exit(-1) ;
160  }
161  if (np > n1f) {
162  cout << "cipcossin: np > n1f : np = " << np << " , n1f = "
163  << n1f << endl ;
164  abort () ;
165  exit(-1) ;
166  }
167  if (n3f > n3c) {
168  cout << "cipcossin: n3f > n3c : n3f = " << n3f << " , n3c = "
169  << n3c << endl ;
170  abort () ;
171  exit(-1) ;
172  }
173  if (n2f > n2c) {
174  cout << "cipcossin: n2f > n2c : n2f = " << n2f << " , n2c = "
175  << n2c << endl ;
176  abort () ;
177  exit(-1) ;
178  }
179 
180  // Recherche des tables
181  int* facto = facto_ini(np) ;
182  double* trigo = trigo_ini(np) ;
183 
184  // Tableau de travail
185  double* t1 = (double*)( malloc( (np+2)*sizeof(double) ) ) ;
186 
187 // Denormalisation des cosinus
188  int n2n3c = n2c * n3c ;
189  for (i=2; i<np; i += 2 ) {
190  for (j=0; j<n2c; j++) {
191  for (k=0; k<n3c; k++) {
192  index = n2n3c * i + n3c * j + k ;
193  cf[index] *= .5 ;
194  }
195  }
196  }
197 
198 // Normalisation des sinus (les termes sin(0*phi) et sin(np/2 *phi) ne sont pas
199 // traites)
200  for (i=3; i<np+1; i += 2 ) {
201  for (j=0; j<n2c; j++) {
202  for (k=0; k<n3c; k++) {
203  index = n2n3c * i + n3c * j + k ;
204  cf[index] *= -.5 ;
205  }
206  }
207  }
208 
209 // Parametres pour la routine FFT991
210  int jump = 1 ;
211  int inc = n2n3c ;
212  int lot = 1 ;
213  int isign = 1 ;
214 
215 // boucle sur r et theta
216 
217  for (j=0; j<n2c; j++) {
218  for (k=0; k<n3c; k++) {
219 
220  index = n3c * j + k ;
221 
222 // FFT inverse
223  double* debut = cf + index ;
224 
225  F77_fft991( debut, t1, trigo, facto, &inc, &jump, &np,
226  &lot, &isign) ;
227  } // fin de la boucle sur r
228  } // fin de la boucle sur theta
229 
230 // On recopie le resultat dans ff si besoin est (c'est-a-dire si le
231 // pointeur ff est different de cf) :
232 
233  if ( ff != cf ) {
234  int n2n3f = n2f * n3f ;
235  for (i=0; i<np; i++) {
236  for (j=0; j<n2f; j++) {
237  for (k=0; k<n3f; k++) {
238  int indexc = n2n3c * i + n3c * j + k ;
239  int indexf = n2n3f * i + n3f * j + k ;
240  ff[indexf] = cf[indexc] ;
241  }
242  }
243  }
244 
245  }
246 
247  // Menage
248  free (t1) ;
249 
250 }
251 }
Lorene prototypes.
Definition: app_hor.h:64