LORENE
comb_lin.C
1 /*
2  * Copyright (c) 1999-2001 Philippe Grandclement
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 comb_lin_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/comb_lin.C,v 1.10 2014/10/13 08:53:28 j_novak Exp $" ;
24 
25 /*
26  * $Id: comb_lin.C,v 1.10 2014/10/13 08:53:28 j_novak Exp $
27  * $Log: comb_lin.C,v $
28  * Revision 1.10 2014/10/13 08:53:28 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.9 2014/10/06 15:16:08 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.8 2008/02/18 13:53:42 j_novak
35  * Removal of special indentation instructions.
36  *
37  * Revision 1.7 2007/12/12 12:30:48 jl_cornou
38  * *** empty log message ***
39  *
40  * Revision 1.6 2007/06/06 07:43:28 p_grandclement
41  * nmax increased to 200
42  *
43  * Revision 1.5 2004/02/09 09:33:56 j_novak
44  * Minor modif.
45  *
46  * Revision 1.4 2004/02/06 10:53:54 j_novak
47  * New dzpuis = 5 -> dzpuis = 3 case (not ready yet).
48  *
49  * Revision 1.3 2002/10/16 14:37:11 j_novak
50  * Reorganization of #include instructions of standard C++, in order to
51  * use experimental version 3 of gcc.
52  *
53  * Revision 1.2 2002/10/09 12:47:31 j_novak
54  * Execution speed improved
55  *
56  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
57  * LORENE
58  *
59  * Revision 2.15 2000/09/07 12:52:04 phil
60  * *** empty log message ***
61  *
62  * Revision 2.14 2000/05/22 13:39:01 phil
63  * ajout du cas dzpuis == 3
64  *
65  * Revision 2.13 2000/01/04 18:59:58 phil
66  * Double nmax
67  *
68  * Revision 2.12 1999/10/12 09:38:52 phil
69  * passage en const Matrice&
70  *
71  * Revision 2.11 1999/10/11 14:26:07 phil
72  * & _> &&
73  *
74  * Revision 2.10 1999/09/30 09:25:36 phil
75  * ajout condition sur dirac=0
76  *
77  * Revision 2.9 1999/09/30 09:21:51 phil
78  * remplacement des && en &
79  *
80  * Revision 2.8 1999/09/17 15:22:47 phil
81  * correction definition NMAX
82  *
83  * Revision 2.7 1999/06/23 12:34:29 phil
84  * ajout de dzpuis = 2
85  *
86  * Revision 2.6 1999/04/28 10:48:19 phil
87  * augmentation de NMAX a 50
88  *
89  * Revision 2.5 1999/04/19 14:01:45 phil
90  * *** empty log message ***
91  *
92  * Revision 2.4 1999/04/16 13:15:45 phil
93  * *** empty log message ***
94  *
95  * Revision 2.3 1999/04/14 13:56:17 phil
96  * Sauvegarde des Matrices deja calculees
97  *
98  * Revision 2.2 1999/04/07 14:37:34 phil
99  * *** empty log message ***
100  *
101  * Revision 2.1 1999/04/07 14:29:51 phil
102  * passage par reference
103  *
104  * Revision 2.0 1999/04/07 14:09:11 phil
105  * *** empty log message ***
106  *
107  *
108  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/comb_lin.C,v 1.10 2014/10/13 08:53:28 j_novak Exp $
109  *
110  */
111 
112 //fichiers includes
113 #include <cstdio>
114 #include <cstdlib>
115 #include <cmath>
116 
117 #include "matrice.h"
118 #include "type_parite.h"
119 #include "proto.h"
120 
121 /* FONCTIONS FAISANT DES COMBINAISONS LINEAIRES DES LIGNES.
122  * Existe en deux versions Tbl ou Matrice.
123  *
124  * Entree :
125  * La Matrice ou le Tbl ( a une dimension)
126  * l : nbre quantique
127  * puis = puissance dans la ZEC
128  * La base de developpement en R
129  *
130  * Sortie :
131  * Renvoie la matrice ou le tableau apres Combinaison lineaire
132  *
133  */
134 
135 // Version Matrice --> Matrice
136 namespace Lorene {
137 Matrice _cl_pas_prevu (const Matrice &source, int l, double echelle, int puis) {
138  cout << "Combinaison lineaire pas prevu..." << endl ;
139  cout << "Source : " << source << endl ;
140  cout << "l : " << l << endl ;
141  cout << "dzpuis : " << puis << endl ;
142  cout << "Echelle : " << echelle << endl ;
143  abort() ;
144  exit(-1) ;
145  return source;
146 }
147 
148 
149  //-------------------
150  //-- R_CHEB ------
151  //-------------------
152 
153 Matrice _cl_r_cheb (const Matrice &source, int l, double echelle, int) {
154  int n = source.get_dim(0) ;assert (n == source.get_dim(1)) ;
155 
156 
157  const int nmax = 200 ; // Nombre de Matrices stockees
158  static Matrice* tab[nmax] ; // les matrices calculees
159  static int nb_dejafait = 0 ; // nbre de matrices calculees
160  static int l_dejafait[nmax] ;
161  static int nr_dejafait[nmax] ;
162  static double vieux_echelle = 0 ;
163 
164  // Si on a change l'echelle : on detruit tout :
165  if (vieux_echelle != echelle) {
166  for (int i=0 ; i<nb_dejafait ; i++) {
167  l_dejafait[i] = -1 ;
168  nr_dejafait[i] = -1 ;
169  delete tab[i] ;
170  }
171  nb_dejafait = 0 ;
172  vieux_echelle = echelle ;
173  }
174 
175  int indice = -1 ;
176 
177  // On determine si la matrice a deja ete calculee :
178  for (int conte=0 ; conte<nb_dejafait ; conte ++)
179  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
180  indice = conte ;
181 
182  // Calcul a faire :
183  if (indice == -1) {
184  if (nb_dejafait >= nmax) {
185  cout << "_cl_r_cheb : trop de matrices" << endl ;
186  abort() ;
187  exit (-1) ;
188  }
189 
190  l_dejafait[nb_dejafait] = l ;
191  nr_dejafait[nb_dejafait] = n ;
192 
193  Matrice barre(source) ;
194  int dirac = 1 ;
195  for (int i=0 ; i<n-2 ; i++) {
196  for (int j=i ; j<(n>(i+7)? i+7 : n) ; j++)
197  barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j))
198  /(i+1) ;
199  if (i==0) dirac = 0 ;
200  }
201 
202  Matrice res(barre) ;
203  for (int i=0 ; i<n-4 ; i++)
204  for (int j=i ; j<(n>(i+5)? i+5 : n) ; j++)
205  res.set(i, j) = barre(i, j)-barre(i+2, j) ;
206  tab[nb_dejafait] = new Matrice(res) ;
207  nb_dejafait ++ ;
208  return res ;
209  }
210 
211  // Cas ou le calcul a deja ete effectue :
212  else
213  return *tab[indice] ;
214 }
215 
216 
217  //-------------------
218  //-- R_JACO02 ------
219  //-------------------
220 
221 Matrice _cl_r_jaco02 (const Matrice &source, int l, double echelle, int) {
222  int n = source.get_dim(0) ;assert (n == source.get_dim(1)) ;
223 
224 
225  const int nmax = 200 ; // Nombre de Matrices stockees
226  static Matrice* tab[nmax] ; // les matrices calculees
227  static int nb_dejafait = 0 ; // nbre de matrices calculees
228  static int l_dejafait[nmax] ;
229  static int nr_dejafait[nmax] ;
230  static double vieux_echelle = 0 ;
231 
232  // Si on a change l'echelle : on detruit tout :
233  if (vieux_echelle != echelle) {
234  for (int i=0 ; i<nb_dejafait ; i++) {
235  l_dejafait[i] = -1 ;
236  nr_dejafait[i] = -1 ;
237  delete tab[i] ;
238  }
239  nb_dejafait = 0 ;
240  vieux_echelle = echelle ;
241  }
242 
243  int indice = -1 ;
244 
245  // On determine si la matrice a deja ete calculee :
246  for (int conte=0 ; conte<nb_dejafait ; conte ++)
247  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
248  indice = conte ;
249 
250  // Calcul a faire :
251  if (indice == -1) {
252  if (nb_dejafait >= nmax) {
253  cout << "_cl_r_jaco02 : trop de matrices" << endl ;
254  abort() ;
255  exit (-1) ;
256  }
257 
258  l_dejafait[nb_dejafait] = l ;
259  nr_dejafait[nb_dejafait] = n ;
260 
261  Matrice barre(source) ;
262  for (int i=0 ; i<n ; i++) {
263  for (int j=i ; j<n ; j++)
264  barre.set(i, j) = source(i, j) ;
265  }
266 
267  Matrice res(barre) ;
268  for (int i=0 ; i<n ; i++)
269  for (int j=i ; j<n ; j++)
270  res.set(i, j) = barre(i, j);
271  tab[nb_dejafait] = new Matrice(res) ;
272  nb_dejafait ++ ;
273  return res ;
274  }
275 
276  // Cas ou le calcul a deja ete effectue :
277  else
278  return *tab[indice] ;
279 }
280 
281 
282  //-------------------
283  //-- R_CHEBP -----
284  //-------------------
285 
286 
287 Matrice _cl_r_chebp (const Matrice &source, int l, double, int) {
288 
289  int n = source.get_dim(0) ;
290  assert (n == source.get_dim(1)) ;
291 
292  const int nmax = 200 ; // Nombre de Matrices stockees
293  static Matrice* tab[nmax] ; // les matrices calculees
294  static int nb_dejafait = 0 ; // nbre de matrices calculees
295  static int l_dejafait[nmax] ;
296  static int nr_dejafait[nmax] ;
297 
298  int indice = -1 ;
299 
300  // On determine si la matrice a deja ete calculee :
301  for (int conte=0 ; conte<nb_dejafait ; conte ++)
302  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
303  indice = conte ;
304 
305  // Calcul a faire :
306  if (indice == -1) {
307  if (nb_dejafait >= nmax) {
308  cout << "_cl_r_chebp : trop de matrices" << endl ;
309  abort() ;
310  exit (-1) ;
311  }
312 
313  l_dejafait[nb_dejafait] = l ;
314  nr_dejafait[nb_dejafait] = n ;
315 
316  Matrice barre(source) ;
317 
318  int dirac = 1 ;
319  for (int i=0 ; i<n-2 ; i++) {
320  for (int j=0 ; j<n ; j++)
321  barre.set(i, j) = (1+dirac)*source(i, j)-source(i+2, j) ;
322  if (i==0) dirac = 0 ;
323  }
324 
325  Matrice tilde(barre) ;
326  for (int i=0 ; i<n-4 ; i++)
327  for (int j=0 ; j<n ; j++)
328  tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
329 
330  Matrice res(tilde) ;
331  for (int i=0 ; i<n-4 ; i++)
332  for (int j=0 ; j<n ; j++)
333  res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
334  tab[nb_dejafait] = new Matrice(res) ;
335  nb_dejafait ++ ;
336  return res ;
337  }
338 
339  // Cas ou le calcul a deja ete effectue :
340  else
341  return *tab[indice] ;
342 }
343 
344  //-------------------
345  //-- R_CHEBI -----
346  //-------------------
347 
348 
349 Matrice _cl_r_chebi (const Matrice &source, int l, double, int) {
350  int n = source.get_dim(0) ;
351  assert (n == source.get_dim(1)) ;
352 
353 
354  const int nmax = 200 ; // Nombre de Matrices stockees
355  static Matrice* tab[nmax] ; // les matrices calculees
356  static int nb_dejafait = 0 ; // nbre de matrices calculees
357  static int l_dejafait[nmax] ;
358  static int nr_dejafait[nmax] ;
359 
360  int indice = -1 ;
361 
362  // On determine si la matrice a deja ete calculee :
363  for (int conte=0 ; conte<nb_dejafait ; conte ++)
364  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
365  indice = conte ;
366 
367  // Calcul a faire :
368  if (indice == -1) {
369  if (nb_dejafait >= nmax) {
370  cout << "_cl_r_chebi : trop de matrices" << endl ;
371  abort() ;
372  exit (-1) ;
373  }
374 
375  l_dejafait[nb_dejafait] = l ;
376  nr_dejafait[nb_dejafait] = n ;
377 
378  Matrice barre(source) ;
379 
380  for (int i=0 ; i<n-2 ; i++)
381  for (int j=0 ; j<n ; j++)
382  barre.set(i, j) = source(i, j)-source(i+2, j) ;
383 
384  Matrice tilde(barre) ;
385  for (int i=0 ; i<n-4 ; i++)
386  for (int j=0 ; j<n ; j++)
387  tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
388 
389  Matrice res(tilde) ;
390  for (int i=0 ; i<n-4 ; i++)
391  for (int j=0 ; j<n ; j++)
392  res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
393  tab[nb_dejafait] = new Matrice(res) ;
394  nb_dejafait ++ ;
395  return res ;
396  }
397 
398  // Cas ou le calcul a deja ete effectue :
399  else
400  return *tab[indice] ;
401 }
402  //-------------------
403  //-- R_CHEBU -----
404  //-------------------
405 
406 Matrice _cl_r_chebu (const Matrice &source, int l, double, int puis) {
407  int n = source.get_dim(0) ;
408  assert (n == source.get_dim(1)) ;
409 
410  Matrice res(n, n) ;
411  res.set_etat_qcq() ;
412 
413  switch (puis) {
414  case 5 :
415  res = _cl_r_chebu_cinq(source, l) ;
416  break ;
417  case 4 :
418  res = _cl_r_chebu_quatre(source, l) ;
419  break ;
420  case 3 :
421  res = _cl_r_chebu_trois (source, l) ;
422  break ;
423  case 2 :
424  res = _cl_r_chebu_deux(source, l) ;
425  break ;
426  default :
427  abort() ;
428  exit(-1) ;
429  }
430 
431  return res ;
432 }
433 
434 
435 // Cas dzpuis = 4
436 Matrice _cl_r_chebu_quatre (const Matrice &source, int l) {
437  int n = source.get_dim(0) ;
438  assert (n == source.get_dim(1)) ;
439 
440 
441  const int nmax = 200 ; // Nombre de Matrices stockees
442  static Matrice* tab[nmax] ; // les matrices calculees
443  static int nb_dejafait = 0 ; // nbre de matrices calculees
444  static int l_dejafait[nmax] ;
445  static int nr_dejafait[nmax] ;
446 
447  int indice = -1 ;
448 
449  // On determine si la matrice a deja ete calculee :
450  for (int conte=0 ; conte<nb_dejafait ; conte ++)
451  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
452  indice = conte ;
453 
454  // Calcul a faire :
455  if (indice == -1) {
456  if (nb_dejafait >= nmax) {
457  cout << "_cl_r_chebu_quatre : trop de matrices" << endl ;
458  abort() ;
459  exit (-1) ;
460  }
461 
462  l_dejafait[nb_dejafait] = l ;
463  nr_dejafait[nb_dejafait] = n ;
464 
465  Matrice barre(source) ;
466 
467  int dirac = 1 ;
468  for (int i=0 ; i<n-2 ; i++) {
469  for (int j=0 ; j<n ; j++)
470  barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
471  if (i==0) dirac = 0 ;
472  }
473 
474  Matrice tilde(barre) ;
475  for (int i=0 ; i<n-4 ; i++)
476  for (int j=0 ; j<n ; j++)
477  tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
478 
479  Matrice prime(tilde) ;
480  for (int i=0 ; i<n-4 ; i++)
481  for (int j=0 ; j<n ; j++)
482  prime.set(i, j) = (tilde(i, j)-tilde(i+1, j)) ;
483 
484  Matrice res(prime) ;
485  for (int i=0 ; i<n-4 ; i++)
486  for (int j=0 ; j<n ; j++)
487  res.set(i, j) = (prime(i, j)-prime(i+2, j)) ;
488  tab[nb_dejafait] = new Matrice(res) ;
489  nb_dejafait ++ ;
490  return res ;
491  }
492 
493  // Cas ou le calcul a deja ete effectue :
494  else
495  return *tab[indice] ;
496 }
497 
498 // Cas dzpuis == 3
499 Matrice _cl_r_chebu_trois (const Matrice &source, int l) {
500  int n = source.get_dim(0) ;
501  assert (n == source.get_dim(1)) ;
502 
503 
504  const int nmax = 200 ; // Nombre de Matrices stockees
505  static Matrice* tab[nmax] ; // les matrices calculees
506  static int nb_dejafait = 0 ; // nbre de matrices calculees
507  static int l_dejafait[nmax] ;
508  static int nr_dejafait[nmax] ;
509 
510  int indice = -1 ;
511 
512  // On determine si la matrice a deja ete calculee :
513  for (int conte=0 ; conte<nb_dejafait ; conte ++)
514  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
515  indice = conte ;
516 
517  // Calcul a faire :
518  if (indice == -1) {
519  if (nb_dejafait >= nmax) {
520  cout << "_cl_r_chebu_trois : trop de matrices" << endl ;
521  abort() ;
522  exit (-1) ;
523  }
524 
525  l_dejafait[nb_dejafait] = l ;
526  nr_dejafait[nb_dejafait] = n ;
527 
528  Matrice barre(source) ;
529 
530  int dirac = 1 ;
531  for (int i=0 ; i<n-2 ; i++) {
532  for (int j=0 ; j<n ; j++)
533  barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
534  if (i==0) dirac = 0 ;
535  }
536 
537  Matrice tilde(barre) ;
538  for (int i=0 ; i<n-4 ; i++)
539  for (int j=0 ; j<n ; j++)
540  tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
541 
542  Matrice res(tilde) ;
543  for (int i=0 ; i<n-4 ; i++)
544  for (int j=0 ; j<n ; j++)
545  res.set(i, j) = (tilde(i, j)+tilde(i+1, j)) ;
546 
547  tab[nb_dejafait] = new Matrice(res) ;
548  nb_dejafait ++ ;
549  return res ;
550  }
551 
552  // Cas ou le calcul a deja ete effectue :
553  else
554  return *tab[indice] ;
555 }
556 
557 
558 //Cas dzpuis == 2
559 Matrice _cl_r_chebu_deux (const Matrice &source, int l) {
560  int n = source.get_dim(0) ;
561  assert (n == source.get_dim(1)) ;
562 
563 
564  const int nmax = 200 ; // Nombre de Matrices stockees
565  static Matrice* tab[nmax] ; // les matrices calculees
566  static int nb_dejafait = 0 ; // nbre de matrices calculees
567  static int l_dejafait[nmax] ;
568  static int nr_dejafait[nmax] ;
569 
570  int indice = -1 ;
571 
572  // On determine si la matrice a deja ete calculee :
573  for (int conte=0 ; conte<nb_dejafait ; conte ++)
574  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
575  indice = conte ;
576 
577  // Calcul a faire :
578  if (indice == -1) {
579  if (nb_dejafait >= nmax) {
580  cout << "_cl_r_chebu_deux : trop de matrices" << endl ;
581  abort() ;
582  exit (-1) ;
583  }
584 
585  l_dejafait[nb_dejafait] = l ;
586  nr_dejafait[nb_dejafait] = n ;
587 
588  Matrice barre(source) ;
589 
590  int dirac = 1 ;
591  for (int i=0 ; i<n-2 ; i++) {
592  for (int j=0 ; j<n ; j++)
593  barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
594  if (i==0) dirac = 0 ;
595  }
596 
597  Matrice tilde(barre) ;
598  for (int i=0 ; i<n-4 ; i++)
599  for (int j=0 ; j<n ; j++)
600  tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
601 
602  Matrice res(tilde) ;
603  for (int i=0 ; i<n-4 ; i++)
604  for (int j=0 ; j<n ; j++)
605  res.set(i, j) = (tilde(i, j)+tilde(i+1, j)) ;
606 
607  return res ;
608  }
609 
610  // Cas ou le calcul a deja ete effectue :
611  else
612  return *tab[indice] ;
613 }
614 
615 
616 //Cas dzpuis == 5
617 Matrice _cl_r_chebu_cinq (const Matrice &source, int l) {
618  int n = source.get_dim(0) ;
619  assert (n == source.get_dim(1)) ;
620 
621 
622  const int nmax = 200 ; // Nombre de Matrices stockees
623  static Matrice* tab[nmax] ; // les matrices calculees
624  static int nb_dejafait = 0 ; // nbre de matrices calculees
625  static int l_dejafait[nmax] ;
626  static int nr_dejafait[nmax] ;
627 
628  int indice = -1 ;
629 
630  // On determine si la matrice a deja ete calculee :
631  for (int conte=0 ; conte<nb_dejafait ; conte ++)
632  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
633  indice = conte ;
634 
635  // Calcul a faire :
636  if (indice == -1) {
637  if (nb_dejafait >= nmax) {
638  cout << "_cl_r_chebu_cinq : trop de matrices" << endl ;
639  abort() ;
640  exit (-1) ;
641  }
642 
643  l_dejafait[nb_dejafait] = l ;
644  nr_dejafait[nb_dejafait] = n ;
645 
646  Matrice barre(source) ;
647 
648  int dirac = 1 ;
649  for (int i=0 ; i<n-2 ; i++) {
650  for (int j=0 ; j<n ; j++)
651  barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
652  if (i==0) dirac = 0 ;
653  }
654 
655  Matrice tilde(barre) ;
656  for (int i=0 ; i<n-4 ; i++)
657  for (int j=0 ; j<n ; j++)
658  tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
659 
660  Matrice res(tilde) ;
661  for (int i=0 ; i<n-4 ; i++)
662  for (int j=0 ; j<n ; j++)
663  res.set(i, j) = (tilde(i, j)+tilde(i+1, j)) ;
664 
665 // cout << "Apres comb. lin. : " << endl ;
666 // cout << res ;
667 // int fg ; cin >> fg ;
668 
669  return res ;
670  }
671 
672  // Cas ou le calcul a deja ete effectue :
673  else
674  return *tab[indice] ;
675 }
676 
677  //-------------------------
678  //- La routine a appeler ---
679  //---------------------------
680 
681 Matrice combinaison (const Matrice &source, int l, double echelle, int puis, int base_r) {
682 
683  // Routines de derivation
684  static Matrice (*combinaison[MAX_BASE])(const Matrice &, int, double, int) ;
685  static int nap = 0 ;
686 
687  // Premier appel
688  if (nap==0) {
689  nap = 1 ;
690  for (int i=0 ; i<MAX_BASE ; i++) {
691  combinaison[i] = _cl_pas_prevu ;
692  }
693  // Les routines existantes
694  combinaison[R_CHEB >> TRA_R] = _cl_r_cheb ;
695  combinaison[R_CHEBU >> TRA_R] = _cl_r_chebu ;
696  combinaison[R_CHEBP >> TRA_R] = _cl_r_chebp ;
697  combinaison[R_CHEBI >> TRA_R] = _cl_r_chebi ;
698  combinaison[R_JACO02 >> TRA_R] = _cl_r_jaco02 ;
699  }
700 
701  Matrice res(combinaison[base_r](source, l, echelle, puis)) ;
702  return res ;
703 }
704 
705 //--------------------------------------------------
706 // Version Tbl --> Tbl a 1D pour la source
707 //--------------------------------------------------
708 
709 
710 Tbl _cl_pas_prevu (const Tbl &source, int puis) {
711  cout << "Combinaison lineaire pas prevue..." << endl ;
712  cout << "source : " << &source << endl ;
713  cout << "dzpuis : " << puis << endl ;
714  abort() ;
715  exit(-1) ;
716  return source;
717 }
718 
719 
720 
721  //-------------------
722  //-- R_CHEB -------
723  //--------------------
724 
725 Tbl _cl_r_cheb (const Tbl &source, int) {
726  Tbl barre(source) ;
727  int n = source.get_dim(0) ;
728 
729  int dirac = 1 ;
730  for (int i=0 ; i<n-2 ; i++) {
731  barre.set(i) = ((1+dirac)*source(i)-source(i+2))
732  /(i+1) ;
733  if (i==0) dirac = 0 ;
734  }
735 
736  Tbl res(barre) ;
737  for (int i=0 ; i<n-4 ; i++)
738  res.set(i) = barre(i)-barre(i+2) ;
739  return res ;
740 }
741 
742 
743  //-------------------
744  //-- R_JACO02 ------
745  //-------------------
746 
747 Tbl _cl_r_jaco02 (const Tbl &source, int) {
748  Tbl barre(source) ;
749  int n = source.get_dim(0) ;
750 
751  for (int i=0 ; i<n ; i++) {
752  barre.set(i) = source(i) ;
753  }
754 
755  Tbl res(barre) ;
756  for (int i=0 ; i<n ; i++)
757  res.set(i) = barre(i);
758  return res ;
759 }
760 
761 
762  //-------------------
763  //-- R_CHEBP -----
764  //-------------------
765 
766 Tbl _cl_r_chebp (const Tbl &source, int) {
767  Tbl barre(source) ;
768  int n = source.get_dim(0) ;
769 
770  int dirac = 1 ;
771  for (int i=0 ; i<n-2 ; i++) {
772  barre.set(i) = (1+dirac)*source(i)-source(i+2) ;
773  if (i==0) dirac = 0 ;
774  }
775 
776  Tbl tilde(barre) ;
777  for (int i=0 ; i<n-4 ; i++)
778  tilde.set(i) = barre(i)-barre(i+2) ;
779 
780  Tbl res(tilde) ;
781  for (int i=0 ; i<n-4 ; i++)
782  res.set(i) = tilde(i)-tilde(i+1) ;
783 
784  return res ;
785 }
786 
787 
788  //-------------------
789  //-- R_CHEBI -----
790  //-------------------
791 
792 Tbl _cl_r_chebi (const Tbl &source, int) {
793  Tbl barre(source) ;
794  int n = source.get_dim(0) ;
795 
796  for (int i=0 ; i<n-2 ; i++)
797  barre.set(i) = source(i)-source(i+2) ;
798 
799  Tbl tilde(barre) ;
800  for (int i=0 ; i<n-4 ; i++)
801  tilde.set(i) = barre(i)-barre(i+2) ;
802 
803  Tbl res(tilde) ;
804  for (int i=0 ; i<n-4 ; i++)
805  res.set(i) = tilde(i)-tilde(i+1) ;
806 
807  return res ;
808 }
809 
810 
811  //-------------------
812  //-- R_CHEBU -----
813  //-------------------
814 
815 Tbl _cl_r_chebu (const Tbl &source, int puis) {
816 
817  int n=source.get_dim(0) ;
818  Tbl res(n) ;
819  res.set_etat_qcq() ;
820 
821  switch(puis) {
822  case 5 :
823  res = _cl_r_chebu_cinq(source) ;
824  break ;
825  case 4 :
826  res = _cl_r_chebu_quatre(source) ;
827  break ;
828  case 3 :
829  res = _cl_r_chebu_trois (source) ;
830  break ;
831  case 2 :
832  res = _cl_r_chebu_deux(source) ;
833  break ;
834 
835  default :
836  abort() ;
837  exit(-1) ;
838  }
839  return res ;
840 }
841 
842 // Cas dzpuis = 4 ;
843 Tbl _cl_r_chebu_quatre (const Tbl &source) {
844  Tbl barre(source) ;
845  int n = source.get_dim(0) ;
846 
847  int dirac = 1 ;
848  for (int i=0 ; i<n-2 ; i++) {
849  barre.set(i) = ((1+dirac)*source(i)-source(i+2)) ;
850  if (i==0) dirac = 0 ;
851  }
852 
853  Tbl tilde(barre) ;
854  for (int i=0 ; i<n-4 ; i++)
855  tilde.set(i) = (barre(i)-barre(i+2)) ;
856 
857  Tbl prime(tilde) ;
858  for (int i=0 ; i<n-4 ; i++)
859  prime.set(i) = (tilde(i)-tilde(i+1)) ;
860 
861  Tbl res(prime) ;
862  for (int i=0 ; i<n-4 ; i++)
863  res.set(i) = (prime(i)-prime(i+2)) ;
864 
865  return res ;
866 }
867 // cas dzpuis = 3
868 Tbl _cl_r_chebu_trois (const Tbl &source) {
869  Tbl barre(source) ;
870  int n = source.get_dim(0) ;
871 
872  int dirac = 1 ;
873  for (int i=0 ; i<n-2 ; i++) {
874  barre.set(i) = ((1+dirac)*source(i)-source(i+2)) ;
875  if (i==0) dirac = 0 ;
876  }
877 
878  Tbl tilde(barre) ;
879  for (int i=0 ; i<n-4 ; i++)
880  tilde.set(i) = (barre(i)-barre(i+2)) ;
881 
882  Tbl res(tilde) ;
883  for (int i=0 ; i<n-4 ; i++)
884  res.set(i) = (tilde(i)+tilde(i+1)) ;
885 
886  return res ;
887 }
888 
889 // Cas dzpuis = 2 ;
890 Tbl _cl_r_chebu_deux (const Tbl &source) {
891  Tbl barre(source) ;
892  int n = source.get_dim(0) ;
893 
894  int dirac = 1 ;
895  for (int i=0 ; i<n-2 ; i++) {
896  barre.set(i) = ((1+dirac)*source(i)-source(i+2)) ;
897  if (i==0) dirac = 0 ;
898  }
899 
900  Tbl tilde(barre) ;
901  for (int i=0 ; i<n-4 ; i++)
902  tilde.set(i) = (barre(i)-barre(i+2)) ;
903 
904  Tbl res(tilde) ;
905  for (int i=0 ; i<n-4 ; i++)
906  res.set(i) = (tilde(i)+tilde(i+1)) ;
907  return res ;
908 }
909 
910 // Cas dzpuis = 5 ;
911 Tbl _cl_r_chebu_cinq (const Tbl &source) {
912  Tbl barre(source) ;
913  int n = source.get_dim(0) ;
914 
915  int dirac = 1 ;
916  for (int i=0 ; i<n-2 ; i++) {
917  barre.set(i) = ((1+dirac)*source(i)-source(i+2)) ;
918  if (i==0) dirac = 0 ;
919  }
920 
921  Tbl tilde(barre) ;
922  for (int i=0 ; i<n-4 ; i++)
923  tilde.set(i) = (barre(i)-barre(i+2)) ;
924 
925  Tbl res(tilde) ;
926  for (int i=0 ; i<n-4 ; i++)
927  res.set(i) = (tilde(i)+tilde(i+1)) ;
928  return res ;
929 }
930 
931 
932  //----------------------------
933  //- Routine a appeler ---
934  //------------------------------
935 
936 Tbl combinaison (const Tbl &source, int puis, int base_r) {
937 
938  // Routines de derivation
939  static Tbl (*combinaison[MAX_BASE])(const Tbl &, int) ;
940  static int nap = 0 ;
941 
942  // Premier appel
943  if (nap==0) {
944  nap = 1 ;
945  for (int i=0 ; i<MAX_BASE ; i++) {
946  combinaison[i] = _cl_pas_prevu ;
947  }
948  // Les routines existantes
949  combinaison[R_CHEB >> TRA_R] = _cl_r_cheb ;
950  combinaison[R_CHEBU >> TRA_R] = _cl_r_chebu ;
951  combinaison[R_CHEBP >> TRA_R] = _cl_r_chebp ;
952  combinaison[R_CHEBI >> TRA_R] = _cl_r_chebi ;
953  combinaison[R_JACO02 >> TRA_R] = _cl_r_jaco02 ;
954  }
955 
956  Tbl res(combinaison[base_r](source, puis)) ;
957  return res ;
958 }
959 }
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Lorene
Lorene prototypes.
Definition: app_hor.h:64
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
Lorene::Matrice::get_dim
int get_dim(int i) const
Returns the dimension of the matrix.
Definition: matrice.C:260
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180