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 
24 char 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 //--------------------------------
98 namespace Lorene {
99 void _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 
109 void _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 
213 void _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 
321 void _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
Lorene prototypes.
Definition: app_hor.h:64