LORENE
cl_pvect_r.C
1/*
2 * Methods for linear combinations for Ope_pois_vect_r
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_pvect_r_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/cl_pvect_r.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $" ;
29
30/*
31 * $Id: cl_pvect_r.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
32 * $Log: cl_pvect_r.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/05/10 15:28:22 j_novak
40 * First version of functions for the solution of the r-component of the
41 * vector Poisson equation.
42 *
43 *
44 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/cl_pvect_r.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
45 *
46 */
47
48//fichiers includes
49#include <cstdlib>
50
51#include "matrice.h"
52#include "type_parite.h"
53
54/* FONCTIONS FAISANT DES COMBINAISONS LINEAIRES DES LIGNES.
55 *
56 * Entree :
57 * La Matrice de l'operateur
58 * l : nbre quantique
59 * puis = puissance dans la ZEC
60 * La base de developpement en R
61 *
62 * Sortie :
63 * Renvoie la matrice apres Combinaison lineaire
64 *
65 */
66
67namespace Lorene {
68Matrice _cl_pvect_r_pas_prevu (const Matrice&, int, double, int) ;
69Matrice _cl_pvect_r_cheb (const Matrice&, int, double, int) ;
70Matrice _cl_pvect_r_chebi (const Matrice&, int, double, int) ;
71Matrice _cl_pvect_r_chebu (const Matrice&, int, double, int) ;
72Matrice _cl_pvect_r_chebp (const Matrice&, int, double, int) ;
73
74// Version Matrice --> Matrice
75Matrice _cl_pvect_r_pas_prevu (const Matrice &source, int l, double echelle, int puis) {
76 cout << "Combinaison lineaire pas prevu..." << endl ;
77 cout << "Source : " << source << endl ;
78 cout << "l : " << l << endl ;
79 cout << "dzpuis : " << puis << endl ;
80 cout << "Echelle : " << echelle << endl ;
81 abort() ;
82 exit(-1) ;
83 return source;
84}
85
86
87 //-------------------
88 //-- R_CHEB ------
89 //-------------------
90
91Matrice _cl_pvect_r_cheb (const Matrice &source, int l, double echelle, int) {
92 int n = source.get_dim(0) ;assert (n == source.get_dim(1)) ;
93
94
95 const int nmax = 100 ; // Nombre de Matrices stockees
96 static Matrice* tab[nmax] ; // les matrices calculees
97 static int nb_dejafait = 0 ; // nbre de matrices calculees
98 static int l_dejafait[nmax] ;
99 static int nr_dejafait[nmax] ;
100 static double vieux_echelle = 0 ;
101
102 // Si on a change l'echelle : on detruit tout :
103 if (vieux_echelle != echelle) {
104 for (int i=0 ; i<nb_dejafait ; i++) {
105 l_dejafait[i] = -1 ;
106 nr_dejafait[i] = -1 ;
107 delete tab[i] ;
108 }
109 nb_dejafait = 0 ;
110 vieux_echelle = echelle ;
111 }
112
113 int indice = -1 ;
114
115 // On determine si la matrice a deja ete calculee :
116 for (int conte=0 ; conte<nb_dejafait ; conte ++)
117 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
118 indice = conte ;
119
120 // Calcul a faire :
121 if (indice == -1) {
122 if (nb_dejafait >= nmax) {
123 cout << "_cl_pvect_r_cheb : trop de matrices" << endl ;
124 abort() ;
125 exit (-1) ;
126 }
127
128 l_dejafait[nb_dejafait] = l ;
129 nr_dejafait[nb_dejafait] = n ;
130
131 Matrice barre(source) ;
132 int dirac = 1 ;
133 for (int i=0 ; i<n-2 ; i++) {
134 for (int j=i ; j<(n>(i+7)? i+7 : n) ; j++)
135 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j))
136 /(i+1) ;
137 if (i==0) dirac = 0 ;
138 }
139
140 Matrice res(barre) ;
141 for (int i=0 ; i<n-4 ; i++)
142 for (int j=i ; j<(n>(i+5)? i+5 : n) ; j++)
143 res.set(i, j) = barre(i, j)-barre(i+2, j) ;
144 tab[nb_dejafait] = new Matrice(res) ;
145 nb_dejafait ++ ;
146 return res ;
147 }
148
149 // Cas ou le calcul a deja ete effectue :
150 else
151 return *tab[indice] ;
152}
153
154 //-------------------
155 //-- R_CHEBP -----
156 //-------------------
157
158
159Matrice _cl_pvect_r_chebp (const Matrice &source, int l, double, int) {
160
161 int n = source.get_dim(0) ;
162 assert (n == source.get_dim(1)) ;
163
164 const int nmax = 100 ; // Nombre de Matrices stockees
165 static Matrice* tab[nmax] ; // les matrices calculees
166 static int nb_dejafait = 0 ; // nbre de matrices calculees
167 static int l_dejafait[nmax] ;
168 static int nr_dejafait[nmax] ;
169
170 int indice = -1 ;
171
172 // On determine si la matrice a deja ete calculee :
173 for (int conte=0 ; conte<nb_dejafait ; conte ++)
174 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
175 indice = conte ;
176
177 // Calcul a faire :
178 if (indice == -1) {
179 if (nb_dejafait >= nmax) {
180 cout << "_cl_pvect_r_chebp : trop de matrices" << endl ;
181 abort() ;
182 exit (-1) ;
183 }
184
185 l_dejafait[nb_dejafait] = l ;
186 nr_dejafait[nb_dejafait] = n ;
187
188 Matrice barre(source) ;
189
190 int dirac = 1 ;
191 for (int i=0 ; i<n-2 ; i++) {
192 for (int j=0 ; j<n ; j++)
193 barre.set(i, j) = (1+dirac)*source(i, j)-source(i+2, j) ;
194 if (i==0) dirac = 0 ;
195 }
196
197 Matrice tilde(barre) ;
198 for (int i=0 ; i<n-4 ; i++)
199 for (int j=0 ; j<n ; j++)
200 tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
201
202 Matrice res(tilde) ;
203 for (int i=0 ; i<n-4 ; i++)
204 for (int j=0 ; j<n ; j++)
205 res.set(i, j) = tilde(i, j)-tilde(i+1, 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 //-- R_CHEBI -----
218 //-------------------
219
220
221Matrice _cl_pvect_r_chebi (const Matrice &source, int l, double, int) {
222 int n = source.get_dim(0) ;
223 assert (n == source.get_dim(1)) ;
224
225
226 const int nmax = 100 ; // Nombre de Matrices stockees
227 static Matrice* tab[nmax] ; // les matrices calculees
228 static int nb_dejafait = 0 ; // nbre de matrices calculees
229 static int l_dejafait[nmax] ;
230 static int nr_dejafait[nmax] ;
231
232 int indice = -1 ;
233
234 // On determine si la matrice a deja ete calculee :
235 for (int conte=0 ; conte<nb_dejafait ; conte ++)
236 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
237 indice = conte ;
238
239 // Calcul a faire :
240 if (indice == -1) {
241 if (nb_dejafait >= nmax) {
242 cout << "_cl_pvect_r_chebi : trop de matrices" << endl ;
243 abort() ;
244 exit (-1) ;
245 }
246
247 l_dejafait[nb_dejafait] = l ;
248 nr_dejafait[nb_dejafait] = n ;
249
250 Matrice barre(source) ;
251
252 for (int i=0 ; i<n-2 ; i++)
253 for (int j=0 ; j<n ; j++)
254 barre.set(i, j) = source(i, j)-source(i+2, j) ;
255
256 Matrice tilde(barre) ;
257 for (int i=0 ; i<n-4 ; i++)
258 for (int j=0 ; j<n ; j++)
259 tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
260
261 Matrice res(tilde) ;
262 for (int i=0 ; i<n-4 ; i++)
263 for (int j=0 ; j<n ; j++)
264 res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
265 tab[nb_dejafait] = new Matrice(res) ;
266 nb_dejafait ++ ;
267 return res ;
268 }
269
270 // Cas ou le calcul a deja ete effectue :
271 else
272 return *tab[indice] ;
273}
274 //-------------------
275 //-- R_CHEBU -----
276 //-------------------
277
278Matrice _cl_pvect_r_chebu (const Matrice &source, int l, double, int puis) {
279 int n = source.get_dim(0) ;
280 assert (n == source.get_dim(1)) ;
281 if (puis != 4) {
282 cout << "_ope_pvect_r_mat_r_chebu : only the case dzpuis = 4 "
283 << '\n' << "is implemented! \n"
284 << "dzpuis = " << puis << endl ;
285 abort() ;
286 }
287
288 const int nmax = 200 ; // Nombre de Matrices stockees
289 static Matrice* tab[nmax] ; // les matrices calculees
290 static int nb_dejafait = 0 ; // nbre de matrices calculees
291 static int l_dejafait[nmax] ;
292 static int nr_dejafait[nmax] ;
293
294 int indice = -1 ;
295
296 // On determine si la matrice a deja ete calculee :
297 for (int conte=0 ; conte<nb_dejafait ; conte ++)
298 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
299 indice = conte ;
300
301 // Calcul a faire :
302 if (indice == -1) {
303 if (nb_dejafait >= nmax) {
304 cout << "_cl_pvect_r_chebu_quatre : trop de matrices" << endl ;
305 abort() ;
306 exit (-1) ;
307 }
308
309 l_dejafait[nb_dejafait] = l ;
310 nr_dejafait[nb_dejafait] = n ;
311
312 Matrice barre(source) ;
313
314 int dirac = 1 ;
315 for (int i=0 ; i<n-2 ; i++) {
316 for (int j=0 ; j<n ; j++)
317 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
318 if (i==0) dirac = 0 ;
319 }
320
321 Matrice tilde(barre) ;
322 for (int i=0 ; i<n-4 ; i++)
323 for (int j=0 ; j<n ; j++)
324 tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
325
326 Matrice prime(tilde) ;
327 for (int i=0 ; i<n-4 ; i++)
328 for (int j=0 ; j<n ; j++)
329 prime.set(i, j) = (tilde(i, j)-tilde(i+1, j)) ;
330
331 Matrice res(prime) ;
332 for (int i=0 ; i<n-4 ; i++)
333 for (int j=0 ; j<n ; j++)
334 res.set(i, j) = (prime(i, j)-prime(i+2, j)) ;
335 tab[nb_dejafait] = new Matrice(res) ;
336 nb_dejafait ++ ;
337 return res ;
338 }
339
340 // Cas ou le calcul a deja ete effectue :
341 else
342 return *tab[indice] ;
343}
344
345
346 //-------------------------
347 //- La routine a appeler ---
348 //---------------------------
349
350Matrice cl_pvect_r(const Matrice &source, int l, double echelle,
351 int puis, int base_r) {
352
353 // Routines de derivation
354 static Matrice (*combinaison[MAX_BASE])(const Matrice &, int, double, int) ;
355 static int nap = 0 ;
356
357 // Premier appel
358 if (nap==0) {
359 nap = 1 ;
360 for (int i=0 ; i<MAX_BASE ; i++) {
361 combinaison[i] = _cl_pvect_r_pas_prevu ;
362 }
363 // Les routines existantes
364 combinaison[R_CHEB >> TRA_R] = _cl_pvect_r_cheb ;
365 combinaison[R_CHEBU >> TRA_R] = _cl_pvect_r_chebu ;
366 combinaison[R_CHEBP >> TRA_R] = _cl_pvect_r_chebp ;
367 combinaison[R_CHEBI >> TRA_R] = _cl_pvect_r_chebi ;
368 }
369
370 Matrice res(combinaison[base_r](source, l, echelle, puis)) ;
371 return res ;
372}
373
374}
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