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
139namespace Lorene {
140Matrice _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
156Matrice _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
224Matrice _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
290Matrice _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
352Matrice _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
409Matrice _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
439Matrice _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
502Matrice _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
562Matrice _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
620Matrice _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
684Matrice 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
713Tbl _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
728Tbl _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
750Tbl _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
769Tbl _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
795Tbl _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
818Tbl _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 ;
846Tbl _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
871Tbl _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 ;
893Tbl _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 ;
914Tbl _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
939Tbl 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::Matrice
Matrix handling.
Definition matrice.h:152
Lorene::Matrice::get_dim
int get_dim(int i) const
Returns the dimension of the matrix.
Definition matrice.C:263
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i]).
Definition tbl.h:403
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
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition type_parite.h:188
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
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene
Lorene prototypes.
Definition app_hor.h:67