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