LORENE
op_dsdphi.C
1/*
2 * Copyright (c) 1999-2000 Jean-Alain Marck
3 * Copyright (c) 1999-2001 Eric Gourgoulhon
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24char op_dsdphi_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_dsdphi.C,v 1.6 2014/10/13 08:53:25 j_novak Exp $" ;
25
26/*
27 * Ensemble des routines de base pour la derivation par rapport a phi
28 * (Utilisation interne)
29 *
30 * void _dsdphi_XXXX(Tbl * t, int & b)
31 * t pointeur sur le Tbl d'entree-sortie
32 * b base spectrale
33 *
34 */
35
36/*
37 * $Id: op_dsdphi.C,v 1.6 2014/10/13 08:53:25 j_novak Exp $
38 * $Log: op_dsdphi.C,v $
39 * Revision 1.6 2014/10/13 08:53:25 j_novak
40 * Lorene classes and functions now belong to the namespace Lorene.
41 *
42 * Revision 1.5 2013/04/25 15:46:06 j_novak
43 * Added special treatment in the case np = 1, for type_p = NONSYM.
44 *
45 * Revision 1.4 2006/03/10 12:45:38 j_novak
46 * Use of C++-style cast.
47 *
48 * Revision 1.3 2003/12/19 16:21:48 j_novak
49 * Shadow hunt
50 *
51 * Revision 1.2 2002/09/09 13:00:40 e_gourgoulhon
52 * Modification of declaration of Fortran 77 prototypes for
53 * a better portability (in particular on IBM AIX systems):
54 * All Fortran subroutine names are now written F77_* and are
55 * defined in the new file C++/Include/proto_f77.h.
56 *
57 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
58 * LORENE
59 *
60 * Revision 2.8 2000/11/14 15:08:16 eric
61 * Traitement du cas np=1 dans P_COSSIN_I.
62 *
63 * Revision 2.7 2000/10/04 15:54:09 eric
64 * *** empty log message ***
65 *
66 * Revision 2.6 2000/10/04 12:50:38 eric
67 * Ajout de la base P_COSSIN_I.
68 *
69 * Revision 2.5 2000/10/04 11:50:32 eric
70 * Ajout d' abort() dans le cas non prevu.
71 *
72 * Revision 2.4 2000/08/18 13:20:01 eric
73 * Traitement du cas np = 1.
74 *
75 * Revision 2.3 2000/02/24 16:41:03 eric
76 * Initialisation a zero du tableau xo avant le calcul.
77 *
78 * Revision 2.2 1999/11/15 16:38:03 eric
79 * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
80 *
81 * Revision 2.1 1999/02/23 11:36:34 hyc
82 * *** empty log message ***
83 *
84 * Revision 2.0 1999/02/22 17:24:59 hyc
85 * *** empty log message ***
86 *
87 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_dsdphi.C,v 1.6 2014/10/13 08:53:25 j_novak Exp $
88 *
89 */
90
91// Headers Lorene
92#include "tbl.h"
93#include "proto_f77.h"
94
95
96// Routine pour les cas non prevus
97//--------------------------------
98namespace Lorene {
99void _dsdphi_pas_prevu(Tbl* , int & b) {
100 cout << "Unknown phi basis in Mtbl_cf::dsdp() !" << endl ;
101 cout << " basis: " << hex << b << endl ;
102 abort() ;
103}
104
105 //------------------//
106 // cas P_COSSIN //
107 //------------------//
108
109void _dsdphi_p_cossin(Tbl* tb, int & )
110{
111
112 // Peut-etre rien a faire ?
113 if (tb->get_etat() == ETATZERO) {
114 return ;
115 }
116
117 // Protection
118 assert(tb->get_etat() == ETATQCQ) ;
119
120 // Pour le confort
121 int nr = (tb->dim).dim[0] ; // Nombre
122 int nt = (tb->dim).dim[1] ; // de points
123 int np = (tb->dim).dim[2] ; // physiques REELS
124 np = np - 2 ; // Nombre de points physiques
125
126 // Cas particulier de la symetrie axiale :
127 if (np == 1) {
128 tb->set_etat_zero() ;
129 return ;
130 }
131
132 // Variables statiques
133 static double* cx = 0 ;
134 static int np_pre =0 ;
135
136 // Test sur np pour initialisation eventuelle
137 //#pragma critical (loch_dsdphi_p_cossin)
138 {
139 if (np > np_pre) {
140 np_pre = np ;
141 cx = reinterpret_cast<double*>(realloc(cx, (np+2) * sizeof(double))) ;
142 for (int i=0 ; i<np+2 ; i += 2) {
143 cx[i] = - (i/2) ;
144 cx[i+1] = (i/2) ;
145 }
146 }
147 } // Fin de region critique
148
149 // pt. sur le tableau de double resultat
150 double* xo = new double[(tb->dim).taille] ;
151
152 // Initialisation a zero :
153 for (int i=0; i<(tb->dim).taille; i++) {
154 xo[i] = 0 ;
155 }
156
157 // On y va...
158 double* xi = tb->t ;
159 double* xci = xi ; // Pointeurs
160 double* xco = xo ; // courants
161 int k, j ;
162
163 // k = 0: inutile, deviendra k=1
164 xci += nr*nt ; // Positionnement
165 xco += nr*nt ;
166
167 // k = 1: 0
168 for (j=0 ; j<nr*nt ; j++) {
169 xco[j] = 0 ;
170 } // Fin de la boucle sur r * theta
171 xci += nr*nt ; // Positionnement
172 xco += nr*nt ;
173
174 // 2 =< k <= np-1
175 for (k=2 ; k<np ; k++) {
176 for (j=0 ; j<nr*nt ; j++) {
177 xco[j] = cx[k] * xci[j] ;
178 } // Fin de la boucle sur r * theta
179 xci += nr*nt ; // Positionnement
180 xco += nr*nt ;
181 } // Fin de la boucle sur phi
182
183 // k = np: inutile, deviendra k=np+1
184 xci += nr*nt ; // Positionnement
185 xco += nr*nt ;
186
187 // k = np+1
188 for (j=0 ; j<nr*nt ; j++) {
189 xco[j] = 0 ;
190 } // Fin de la boucle sur r * theta
191
192 // Inversion des sinus et cosinus (appel a dswap de BLAS)
193 int nbr = nr*nt ;
194 int inc = 1 ;
195 double* p1 = xo ;
196 double* p2 = p1 + nr*nt ;
197 for (k=0 ; k<np+1 ; k +=2, p1 += 2*nr*nt, p2 += 2*nr*nt) {
198 F77_dswap(&nbr, p1, &inc, p2, &inc) ;
199 }
200
201 // On remet les choses la ou il faut
202 delete [] tb->t ;
203 tb->t = xo ;
204
205 // base de developpement
206 // inchangee
207}
208
209 //------------------//
210 // cas P_COSSIN_P //
211 //------------------//
212
213void _dsdphi_p_cossin_p(Tbl* tb, int & )
214{
215
216 // Peut-etre rien a faire ?
217 if (tb->get_etat() == ETATZERO) {
218 return ;
219 }
220
221 // Protection
222 assert(tb->get_etat() == ETATQCQ) ;
223
224 // Pour le confort
225 int nr = (tb->dim).dim[0] ; // Nombre
226 int nt = (tb->dim).dim[1] ; // de points
227 int np = (tb->dim).dim[2] ; // physiques REELS
228 np = np - 2 ; // Nombre de points physiques
229
230 // Cas particulier de la symetrie axiale :
231 if (np == 1) {
232 tb->set_etat_zero() ;
233 return ;
234 }
235
236 // Variables statiques
237 static double* cx = 0 ;
238 static int np_pre =0 ;
239
240 // Test sur np pour initialisation eventuelle
241 //#pragma critical (loch_dsdphi_p_cossin_p)
242 {
243 if (np > np_pre) {
244 np_pre = np ;
245 cx = reinterpret_cast<double*>(realloc(cx, (np+2) * sizeof(double))) ;
246 int i ;
247 for (i=0 ; i<np+2 ; i += 2) {
248 cx[i] = - (i/2) ;
249 cx[i+1] = (i/2) ;
250 }
251 for (i=0 ; i<np+2 ; i++) {
252 cx[i] = 2 * cx[i] ;
253 }
254 }
255 } // Fin de region critique
256
257 // pt. sur le tableau de double resultat
258 double* xo = new double[(tb->dim).taille] ;
259
260 // Initialisation a zero :
261 for (int i=0; i<(tb->dim).taille; i++) {
262 xo[i] = 0 ;
263 }
264
265 // On y va...
266 double* xi = tb->t ;
267 double* xci = xi ; // Pointeurs
268 double* xco = xo ; // courants
269 int k, j ;
270
271 // k = 0: inutile, deviendra k=1
272 xci += nr*nt ; // Positionnement
273 xco += nr*nt ;
274
275 // k = 1: 0
276 for (j=0 ; j<nr*nt ; j++) {
277 xco[j] = 0 ;
278 } // Fin de la boucle sur r * theta
279 xci += nr*nt ; // Positionnement
280 xco += nr*nt ;
281
282 // 2 =< k <= np-1
283 for (k=2 ; k<np ; k++) {
284 for (j=0 ; j<nr*nt ; j++) {
285 xco[j] = cx[k] * xci[j] ;
286 } // Fin de la boucle sur r * theta
287 xci += nr*nt ; // Positionnement
288 xco += nr*nt ;
289 } // Fin de la boucle sur phi
290
291 // k = np: inutile, deviendra k=np+1
292 xci += nr*nt ; // Positionnement
293 xco += nr*nt ;
294
295 // k = np+1
296 for (j=0 ; j<nr*nt ; j++) {
297 xco[j] = 0 ;
298 } // Fin de la boucle sur r * theta
299
300 // Inversion des sinus et cosinus (appel a dswap de BLAS)
301 int nbr = nr*nt ;
302 int inc = 1 ;
303 double* p1 = xo ;
304 double* p2 = p1 + nr*nt ;
305 for (k=0 ; k<np+1 ; k +=2, p1 += 2*nr*nt, p2 += 2*nr*nt) {
306 F77_dswap(&nbr, p1, &inc, p2, &inc) ;
307 }
308
309 // On remet les choses la ou il faut
310 delete [] tb->t ;
311 tb->t = xo ;
312
313 // base de developpement
314 // inchangee
315}
316
317 //------------------//
318 // cas P_COSSIN_I //
319 //------------------//
320
321void _dsdphi_p_cossin_i(Tbl* tb, int & )
322{
323
324 // Peut-etre rien a faire ?
325 if (tb->get_etat() == ETATZERO) {
326 return ;
327 }
328
329 // Protection
330 assert(tb->get_etat() == ETATQCQ) ;
331
332 // Pour le confort
333 int nr = (tb->dim).dim[0] ; // Nombre
334 int nt = (tb->dim).dim[1] ; // de points
335 int np = (tb->dim).dim[2] ; // physiques REELS
336 np = np - 2 ; // Nombre de points physiques
337
338 int ntnr = nt*nr ; // increment en phi
339
340 // Cas particulier de la symetrie axiale :
341 if (np == 1) {
342 tb->set_etat_zero() ;
343 return ;
344 }
345
346 // pt. sur le tableau de double resultat
347 double* xo = new double[(tb->dim).taille] ;
348
349 // pt. sur le tableau de double entree
350 const double* xi = tb->t ;
351
352 // k = 0 : resultat sur cos(phi)
353 // ------------------------------
354
355 double* xco = xo ; // coef de cos(phi)
356 const double* xci = xi + 2*ntnr ; // coef de sin(phi)
357 for (int i=0; i<ntnr; i++) {
358 xco[i] = xci[i] ;
359 }
360
361 // k = 1 : mise a zero
362 // --------------------
363
364 xco = xo + ntnr ;
365 for (int i=0; i<ntnr; i++) {
366 xco[i] = 0 ;
367 }
368
369 // k = 2 : resultat sur sin(phi)
370 // ------------------------------
371
372 xco = xo + 2*ntnr ; // coef de sin(phi)
373 xci = xi ; // coef de cos(phi)
374 for (int i=0; i<ntnr; i++) {
375 xco[i] = - xci[i] ;
376 }
377
378 // k >= 3
379 // ------
380
381 for (int k=3; k<np; k+=2) {
382
383 // resultat sur cos(k phi)
384 // -----------------------
385
386 xco = xo + k*ntnr ; // coef de cos(k phi)
387 xci = xi + (k+1)*ntnr ; // coef de sin(k phi)
388
389 for (int i=0; i<ntnr; i++) {
390 xco[i] = k * xci[i] ;
391 }
392
393 // resultat sur sin(k phi)
394 // -----------------------
395
396 xco = xo + (k+1)*ntnr ; // coef de sin(k phi)
397 xci = xi + k*ntnr ; // coef de cos(k phi)
398
399 for (int i=0; i<ntnr; i++) {
400 xco[i] = - k * xci[i] ;
401 }
402
403 }
404
405 // k = np+1 : mise a zero
406 // -----------------------
407
408 xco = xo + (np+1)*ntnr ;
409 for (int i=0; i<ntnr; i++) {
410 xco[i] = 0 ;
411 }
412
413 // On remet les choses la ou il faut
414 // ---------------------------------
415 delete [] tb->t ;
416 tb->t = xo ;
417
418 // base de developpement
419 // inchangee
420}
421}
Lorene prototypes.
Definition app_hor.h:64