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
23char 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
136namespace Lorene {
137Matrice _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
153Matrice _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
221Matrice _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
287Matrice _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
349Matrice _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
406Matrice _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
436Matrice _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
499Matrice _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
559Matrice _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
617Matrice _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
681Matrice 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
710Tbl _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
725Tbl _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
747Tbl _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
766Tbl _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
792Tbl _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
815Tbl _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 ;
843Tbl _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
868Tbl _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 ;
890Tbl _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 ;
911Tbl _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
936Tbl 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}
int get_dim(int i) const
Returns the dimension of the matrix.
Definition matrice.C:260
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64