LORENE
cl_ptens_rr.C
1/*
2 * Methods for linear combinations for Ope_pois_tens_rr
3 *
4 * (see file ope_elementary.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2004 Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char cl_ptens_rr_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/cl_ptens_rr.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $" ;
29
30/*
31 * $Id: cl_ptens_rr.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
32 * $Log: cl_ptens_rr.C,v $
33 * Revision 1.3 2014/10/13 08:53:34 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.2 2014/10/06 15:16:13 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.1 2004/12/23 16:30:15 j_novak
40 * New files and class for the solution of the rr component of the tensor Poisson
41 * equation.
42 *
43 *
44 *
45 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/cl_ptens_rr.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
46 *
47 */
48
49//fichiers includes
50#include <cstdlib>
51
52#include "matrice.h"
53#include "type_parite.h"
54
55/* FONCTIONS FAISANT DES COMBINAISONS LINEAIRES DES LIGNES.
56 *
57 * Entree :
58 * La Matrice de l'operateur
59 * l : nbre quantique
60 * puis = puissance dans la ZEC
61 * La base de developpement en R
62 *
63 * Sortie :
64 * Renvoie la matrice apres Combinaison lineaire
65 *
66 */
67
68namespace Lorene {
69Matrice _cl_ptens_rr_pas_prevu (const Matrice&, int, double, int) ;
70Matrice _cl_ptens_rr_cheb (const Matrice&, int, double, int) ;
71Matrice _cl_ptens_rr_chebi (const Matrice&, int, double, int) ;
72Matrice _cl_ptens_rr_chebu (const Matrice&, int, double, int) ;
73Matrice _cl_ptens_rr_chebp (const Matrice&, int, double, int) ;
74
75// Version Matrice --> Matrice
76Matrice _cl_ptens_rr_pas_prevu (const Matrice &source, int l, double echelle, int puis) {
77 cout << "Combinaison lineaire pas prevu..." << endl ;
78 cout << "Source : " << source << endl ;
79 cout << "l : " << l << endl ;
80 cout << "dzpuis : " << puis << endl ;
81 cout << "Echelle : " << echelle << endl ;
82 abort() ;
83 exit(-1) ;
84 return source;
85}
86
87
88 //-------------------
89 //-- R_CHEB ------
90 //-------------------
91
92Matrice _cl_ptens_rr_cheb (const Matrice &source, int l, double echelle, int) {
93 int n = source.get_dim(0) ;assert (n == source.get_dim(1)) ;
94
95
96 const int nmax = 100 ; // Nombre de Matrices stockees
97 static Matrice* tab[nmax] ; // les matrices calculees
98 static int nb_dejafait = 0 ; // nbre de matrices calculees
99 static int l_dejafait[nmax] ;
100 static int nr_dejafait[nmax] ;
101 static double vieux_echelle = 0 ;
102
103 // Si on a change l'echelle : on detruit tout :
104 if (vieux_echelle != echelle) {
105 for (int i=0 ; i<nb_dejafait ; i++) {
106 l_dejafait[i] = -1 ;
107 nr_dejafait[i] = -1 ;
108 delete tab[i] ;
109 }
110 nb_dejafait = 0 ;
111 vieux_echelle = echelle ;
112 }
113
114 int indice = -1 ;
115
116 // On determine si la matrice a deja ete calculee :
117 for (int conte=0 ; conte<nb_dejafait ; conte ++)
118 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
119 indice = conte ;
120
121 // Calcul a faire :
122 if (indice == -1) {
123 if (nb_dejafait >= nmax) {
124 cout << "_cl_ptens_rr_cheb : trop de matrices" << endl ;
125 abort() ;
126 exit (-1) ;
127 }
128
129 l_dejafait[nb_dejafait] = l ;
130 nr_dejafait[nb_dejafait] = n ;
131
132 Matrice barre(source) ;
133 int dirac = 1 ;
134 for (int i=0 ; i<n-2 ; i++) {
135 for (int j=i ; j<(n>(i+7)? i+7 : n) ; j++)
136 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j))
137 /(i+1) ;
138 if (i==0) dirac = 0 ;
139 }
140
141 Matrice res(barre) ;
142 for (int i=0 ; i<n-4 ; i++)
143 for (int j=i ; j<(n>(i+5)? i+5 : n) ; j++)
144 res.set(i, j) = barre(i, j)-barre(i+2, j) ;
145 tab[nb_dejafait] = new Matrice(res) ;
146 nb_dejafait ++ ;
147 return res ;
148 }
149
150 // Cas ou le calcul a deja ete effectue :
151 else
152 return *tab[indice] ;
153}
154
155 //-------------------
156 //-- R_CHEBP -----
157 //-------------------
158
159
160Matrice _cl_ptens_rr_chebp (const Matrice &source, int l, double, int) {
161
162 int n = source.get_dim(0) ;
163 assert (n == source.get_dim(1)) ;
164
165 const int nmax = 100 ; // Nombre de Matrices stockees
166 static Matrice* tab[nmax] ; // les matrices calculees
167 static int nb_dejafait = 0 ; // nbre de matrices calculees
168 static int l_dejafait[nmax] ;
169 static int nr_dejafait[nmax] ;
170
171 int indice = -1 ;
172
173 // On determine si la matrice a deja ete calculee :
174 for (int conte=0 ; conte<nb_dejafait ; conte ++)
175 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
176 indice = conte ;
177
178 // Calcul a faire :
179 if (indice == -1) {
180 if (nb_dejafait >= nmax) {
181 cout << "_cl_ptens_rr_chebp : trop de matrices" << endl ;
182 abort() ;
183 exit (-1) ;
184 }
185
186 l_dejafait[nb_dejafait] = l ;
187 nr_dejafait[nb_dejafait] = n ;
188
189 Matrice barre(source) ;
190
191 int dirac = 1 ;
192 for (int i=0 ; i<n-2 ; i++) {
193 for (int j=0 ; j<n ; j++)
194 barre.set(i, j) = (1+dirac)*source(i, j)-source(i+2, j) ;
195 if (i==0) dirac = 0 ;
196 }
197
198 Matrice tilde(barre) ;
199 for (int i=0 ; i<n-4 ; i++)
200 for (int j=0 ; j<n ; j++)
201 tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
202
203 Matrice res(tilde) ;
204 for (int i=0 ; i<n-4 ; i++)
205 for (int j=0 ; j<n ; j++)
206 res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
207 tab[nb_dejafait] = new Matrice(res) ;
208 nb_dejafait ++ ;
209 return res ;
210 }
211
212 // Cas ou le calcul a deja ete effectue :
213 else
214 return *tab[indice] ;
215}
216
217 //-------------------
218 //-- R_CHEBI -----
219 //-------------------
220
221
222Matrice _cl_ptens_rr_chebi (const Matrice &source, int l, double, int) {
223 int n = source.get_dim(0) ;
224 assert (n == source.get_dim(1)) ;
225
226
227 const int nmax = 100 ; // Nombre de Matrices stockees
228 static Matrice* tab[nmax] ; // les matrices calculees
229 static int nb_dejafait = 0 ; // nbre de matrices calculees
230 static int l_dejafait[nmax] ;
231 static int nr_dejafait[nmax] ;
232
233 int indice = -1 ;
234
235 // On determine si la matrice a deja ete calculee :
236 for (int conte=0 ; conte<nb_dejafait ; conte ++)
237 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
238 indice = conte ;
239
240 // Calcul a faire :
241 if (indice == -1) {
242 if (nb_dejafait >= nmax) {
243 cout << "_cl_ptens_rr_chebi : trop de matrices" << endl ;
244 abort() ;
245 exit (-1) ;
246 }
247
248 l_dejafait[nb_dejafait] = l ;
249 nr_dejafait[nb_dejafait] = n ;
250
251 Matrice barre(source) ;
252
253 for (int i=0 ; i<n-2 ; i++)
254 for (int j=0 ; j<n ; j++)
255 barre.set(i, j) = source(i, j)-source(i+2, j) ;
256
257 Matrice tilde(barre) ;
258 for (int i=0 ; i<n-4 ; i++)
259 for (int j=0 ; j<n ; j++)
260 tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
261
262 Matrice res(tilde) ;
263 for (int i=0 ; i<n-4 ; i++)
264 for (int j=0 ; j<n ; j++)
265 res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
266 tab[nb_dejafait] = new Matrice(res) ;
267 nb_dejafait ++ ;
268 return res ;
269 }
270
271 // Cas ou le calcul a deja ete effectue :
272 else
273 return *tab[indice] ;
274}
275 //-------------------
276 //-- R_CHEBU -----
277 //-------------------
278
279Matrice _cl_ptens_rr_chebu (const Matrice &source, int l, double, int puis) {
280 int n = source.get_dim(0) ;
281 assert (n == source.get_dim(1)) ;
282 if (puis != 4) {
283 cout << "_ope_ptens_rr_mat_r_chebu : only the case dzpuis = 4 "
284 << '\n' << "is implemented! \n"
285 << "dzpuis = " << puis << endl ;
286 abort() ;
287 }
288
289 const int nmax = 200 ; // Nombre de Matrices stockees
290 static Matrice* tab[nmax] ; // les matrices calculees
291 static int nb_dejafait = 0 ; // nbre de matrices calculees
292 static int l_dejafait[nmax] ;
293 static int nr_dejafait[nmax] ;
294
295 int indice = -1 ;
296
297 // On determine si la matrice a deja ete calculee :
298 for (int conte=0 ; conte<nb_dejafait ; conte ++)
299 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
300 indice = conte ;
301
302 // Calcul a faire :
303 if (indice == -1) {
304 if (nb_dejafait >= nmax) {
305 cout << "_cl_ptens_rr_chebu_quatre : trop de matrices" << endl ;
306 abort() ;
307 exit (-1) ;
308 }
309
310 l_dejafait[nb_dejafait] = l ;
311 nr_dejafait[nb_dejafait] = n ;
312
313 Matrice barre(source) ;
314
315 int dirac = 1 ;
316 for (int i=0 ; i<n-2 ; i++) {
317 for (int j=0 ; j<n ; j++)
318 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
319 if (i==0) dirac = 0 ;
320 }
321
322 Matrice tilde(barre) ;
323 for (int i=0 ; i<n-4 ; i++)
324 for (int j=0 ; j<n ; j++)
325 tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
326
327 Matrice prime(tilde) ;
328 for (int i=0 ; i<n-4 ; i++)
329 for (int j=0 ; j<n ; j++)
330 prime.set(i, j) = (tilde(i, j)-tilde(i+1, j)) ;
331
332 Matrice res(prime) ;
333 for (int i=0 ; i<n-4 ; i++)
334 for (int j=0 ; j<n ; j++)
335 res.set(i, j) = (prime(i, j)-prime(i+2, j)) ;
336 tab[nb_dejafait] = new Matrice(res) ;
337 nb_dejafait ++ ;
338 return res ;
339 }
340
341 // Cas ou le calcul a deja ete effectue :
342 else
343 return *tab[indice] ;
344}
345
346
347 //-------------------------
348 //- La routine a appeler ---
349 //---------------------------
350
351Matrice cl_ptens_rr(const Matrice &source, int l, double echelle,
352 int puis, int base_r) {
353
354 // Routines de derivation
355 static Matrice (*combinaison[MAX_BASE])(const Matrice &, int, double, int) ;
356 static int nap = 0 ;
357
358 // Premier appel
359 if (nap==0) {
360 nap = 1 ;
361 for (int i=0 ; i<MAX_BASE ; i++) {
362 combinaison[i] = _cl_ptens_rr_pas_prevu ;
363 }
364 // Les routines existantes
365 combinaison[R_CHEB >> TRA_R] = _cl_ptens_rr_cheb ;
366 combinaison[R_CHEBU >> TRA_R] = _cl_ptens_rr_chebu ;
367 combinaison[R_CHEBP >> TRA_R] = _cl_ptens_rr_chebp ;
368 combinaison[R_CHEBI >> TRA_R] = _cl_ptens_rr_chebi ;
369 }
370
371 Matrice res(combinaison[base_r](source, l, echelle, puis)) ;
372 return res ;
373}
374
375}
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_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