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