LORENE
prepa_pvect_r.C
1/*
2 * Methods preparing the operators of Ope_pois_vect_r for inversion.
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 prepa_pvect_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/prepa_pvect_r.C,v 1.4 2014/10/13 08:53:34 j_novak Exp $" ;
29
30/*
31 * $Id: prepa_pvect_r.C,v 1.4 2014/10/13 08:53:34 j_novak Exp $
32 * $Log: prepa_pvect_r.C,v $
33 * Revision 1.4 2014/10/13 08:53:34 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.3 2014/10/06 15:16:13 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.2 2004/05/17 15:42:21 j_novak
40 * The method 1 of vector Poisson eq. solves directly for F^r.
41 * Some bugs were corrected in the operator poisson_vect.
42 *
43 * Revision 1.1 2004/05/10 15:28:22 j_novak
44 * First version of functions for the solution of the r-component of the
45 * vector Poisson equation.
46 *
47 *
48 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/prepa_pvect_r.C,v 1.4 2014/10/13 08:53:34 j_novak Exp $
49 *
50 */
51
52//fichiers includes
53#include <cstdlib>
54
55#include "matrice.h"
56#include "type_parite.h"
57
58/*
59 * Fonctions supprimant le nombre de colonnes (les premieres)
60 et de lignes (les dernieres) a l'operateur renvoye par laplacien_mat, de facon
61 a ce qu'il ne soit plus degenere. Ceci doit etre fait apres les combinaisons
62 lineaires. La mise a bandes et la decomposition LU sont egalement effectuees ici
63
64
65 Entree : lap : resultat de laplacien_mat
66 l : associe a lap
67 puis : puissance dans la ZEC
68 base_r : base de developpement
69
70 Sortie : renvoie un operateur non degenere ....
71 */
72
73namespace Lorene {
74Matrice _nondeg_pvect_r_pas_prevu(const Matrice &, int , double, int) ;
75Matrice _nondeg_pvect_r_cheb (const Matrice&, int, double, int) ;
76Matrice _nondeg_pvect_r_chebp (const Matrice&, int, double, int) ;
77Matrice _nondeg_pvect_r_chebi (const Matrice&, int, double, int) ;
78Matrice _nondeg_pvect_r_chebu (const Matrice&, int, double, int) ;
79
80
81 //------------------------------------
82 // Routine pour les cas non prevus --
83 //-----------------------------------
84
85Matrice _nondeg_pvect_r_pas_prevu(const Matrice &lap, int l, double echelle, int puis) {
86 cout << "Construction non degeneree pas prevue..." << endl ;
87 cout << "l : " << l << endl ;
88 cout << "lap : " << lap << endl ;
89 cout << "echelle : " << echelle << endl ;
90 cout << " puis : " << puis << endl ;
91 abort() ;
92 exit(-1) ;
93 Matrice res(1, 1) ;
94 return res;
95}
96
97
98
99 //-------------------
100 //-- R_CHEB -------
101 //--------------------
102
103Matrice _nondeg_pvect_r_cheb (const Matrice &lap, int l, double echelle, int) {
104
105
106 int n = lap.get_dim(0) ;
107
108 const int nmax = 200 ; // Nombre de Matrices stockees
109 static Matrice* tab[nmax] ; // les matrices calculees
110 static int nb_dejafait = 0 ; // nbre de matrices calculees
111 static int l_dejafait[nmax] ;
112 static int nr_dejafait[nmax] ;
113 static double vieux_echelle = 0;
114
115 // Si on a change l'echelle : on detruit tout :
116 if (vieux_echelle != echelle) {
117 for (int i=0 ; i<nb_dejafait ; i++) {
118 l_dejafait[i] = -1 ;
119 nr_dejafait[i] = -1 ;
120 delete tab[i] ;
121 }
122 vieux_echelle = echelle ;
123 nb_dejafait = 0 ;
124 }
125
126 int indice = -1 ;
127
128 // On determine si la matrice a deja ete calculee :
129 for (int conte=0 ; conte<nb_dejafait ; conte ++)
130 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
131 indice = conte ;
132
133 // Calcul a faire :
134 if (indice == -1) {
135 if (nb_dejafait >= nmax) {
136 cout << "_nondeg_pvect_r_cheb : trop de matrices" << endl ;
137 abort() ;
138 exit (-1) ;
139 }
140
141
142 l_dejafait[nb_dejafait] = l ;
143 nr_dejafait[nb_dejafait] = n ;
144
145
146 //assert (l<n) ;
147
148 Matrice res(n-2, n-2) ;
149 res.set_etat_qcq() ;
150 for (int i=0 ; i<n-2 ; i++)
151 for (int j=0 ; j<n-2 ; j++)
152 res.set(i, j) = lap(i, j+2) ;
153
154 res.set_band(2, 2) ;
155 res.set_lu() ;
156 tab[nb_dejafait] = new Matrice(res) ;
157 nb_dejafait ++ ;
158 return res ;
159 }
160
161 // Cas ou le calcul a deja ete effectue :
162 else
163 return *tab[indice] ;
164}
165
166
167
168
169 //------------------
170 //-- R_CHEBP ----
171 //------------------
172
173Matrice _nondeg_pvect_r_chebp (const Matrice &lap, int l, double, int) {
174
175 int n = lap.get_dim(0) ;
176
177
178 const int nmax = 200 ; // Nombre de Matrices stockees
179 static Matrice* tab[nmax] ; // les matrices calculees
180 static int nb_dejafait = 0 ; // nbre de matrices calculees
181 static int l_dejafait[nmax] ;
182 static int nr_dejafait[nmax] ;
183
184 int indice = -1 ;
185
186 // On determine si la matrice a deja ete calculee :
187 for (int conte=0 ; conte<nb_dejafait ; conte ++)
188 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
189 indice = conte ;
190
191 // Calcul a faire :
192 if (indice == -1) {
193 if (nb_dejafait >= nmax) {
194 cout << "_nondeg_pvect_r_chebp : trop de matrices" << endl ;
195 abort() ;
196 exit (-1) ;
197 }
198
199
200 l_dejafait[nb_dejafait] = l ;
201 nr_dejafait[nb_dejafait] = n ;
202
203 assert (div(l, 2).rem == 1) ;
204
205 if (l==1) {
206 Matrice res(n-1, n-1) ;
207 res.set_etat_qcq() ;
208 for (int i=0 ; i<n-1 ; i++)
209 for (int j=0 ; j<n-1 ; j++)
210 res.set(i, j) = lap(i, j+1) ;
211 res.set_band(3, 0) ;
212 res.set_lu() ;
213 tab[nb_dejafait] = new Matrice(res) ;
214 nb_dejafait ++ ;
215 return res ;
216 }
217 else {
218 Matrice res(n-2, n-2) ;
219 res.set_etat_qcq() ;
220 for (int i=0 ;i<n-2 ; i++)
221 for (int j=0 ; j<n-2 ; j++)
222 res.set(i, j) = lap(i, j+2) ;
223
224 res.set_band(2, 1) ;
225 res.set_lu() ;
226 tab[nb_dejafait] = new Matrice(res) ;
227 nb_dejafait ++ ;
228 return res ;
229 }
230 }
231 // Cas ou le calcul a deja ete effectue :
232 else
233 return *tab[indice] ;
234}
235
236
237
238
239 //-------------------
240 //-- R_CHEBI -----
241 //-------------------
242
243Matrice _nondeg_pvect_r_chebi (const Matrice &lap, int l, double, int) {
244
245 int n = lap.get_dim(0) ;
246
247 const int nmax = 200 ; // Nombre de Matrices stockees
248 static Matrice* tab[nmax] ; // les matrices calculees
249 static int nb_dejafait = 0 ; // nbre de matrices calculees
250 static int l_dejafait[nmax] ;
251 static int nr_dejafait[nmax] ;
252
253 int indice = -1 ;
254
255 // On determine si la matrice a deja ete calculee :
256 for (int conte=0 ; conte<nb_dejafait ; conte ++)
257 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
258 indice = conte ;
259
260 // Calcul a faire :
261 if (indice == -1) {
262 if (nb_dejafait >= nmax) {
263 cout << "_nondeg_pvect_r_chebi : trop de matrices" << endl ;
264 abort() ;
265 exit (-1) ;
266 }
267
268
269 l_dejafait[nb_dejafait] = l ;
270 nr_dejafait[nb_dejafait] = n ;
271
272
273 assert (div(l, 2).rem == 0) ;
274 // assert (l<=2*n-1) ;
275
276 if (l<=2) { //### to be checked!!!
277 Matrice res(n-1, n-1) ;
278 res.set_etat_qcq() ;
279 for (int i=0 ; i<n-1 ; i++)
280 for (int j=0 ; j<n-1 ; j++)
281 res.set(i, j) = lap(i, j+1) ;
282 res.set_band(3, 0) ;
283 res.set_lu() ;
284 tab[nb_dejafait] = new Matrice(res) ;
285 nb_dejafait ++ ;
286 return res ;
287 }
288 else {
289 Matrice res(n-2, n-2) ;
290 res.set_etat_qcq() ;
291 for (int i=0 ;i<n-2 ; i++)
292 for (int j=0 ; j<n-2 ; j++)
293 res.set(i, j) = lap(i, j+2) ;
294
295 res.set_band(2, 1) ;
296 res.set_lu() ;
297 tab[nb_dejafait] = new Matrice(res) ;
298 nb_dejafait ++ ;
299 return res ;
300 }
301 }
302 // Cas ou le calcul a deja ete effectue :
303 else
304 return *tab[indice] ;
305}
306
307
308
309
310 //-------------------
311 //-- R_CHEBU -----
312 //-------------------
313
314
315Matrice _nondeg_pvect_r_chebu (const Matrice &lap, int l, double, int puis) {
316
317 if (puis != 4) {
318 cout << "_ope_pvect_r_mat_r_chebu : only the case dzpuis = 4 "
319 << '\n' << "is implemented! \n"
320 << "dzpuis = " << puis << endl ;
321 abort() ;
322 }
323 int n = lap.get_dim(0) ;
324
325 const int nmax = 200; // Nombre de Matrices stockees
326 static Matrice* tab[nmax] ; // les matrices calculees
327 static int nb_dejafait = 0 ; // nbre de matrices calculees
328 static int l_dejafait[nmax] ;
329 static int nr_dejafait[nmax] ;
330
331 int indice = -1 ;
332
333 // On determine si la matrice a deja ete calculee :
334 for (int conte=0 ; conte<nb_dejafait ; conte ++)
335 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
336 indice = conte ;
337
338 // Calcul a faire :
339 if (indice == -1) {
340 if (nb_dejafait >= nmax) {
341 cout << "_nondeg_pvect_r_chebu : trop de matrices" << endl ;
342 abort() ;
343 exit (-1) ;
344 }
345
346 l_dejafait[nb_dejafait] = l ;
347 nr_dejafait[nb_dejafait] = n ;
348
349 Matrice res(n-3, n-3) ;
350 res.set_etat_qcq() ;
351 for (int i=0 ;i<n-3 ; i++)
352 for (int j=0 ; j<n-3 ; j++)
353 res.set(i, j) = lap(i, j+3) ;
354
355 res.set_band(2, 1) ;
356 res.set_lu() ;
357 tab[nb_dejafait] = new Matrice(res) ;
358 nb_dejafait ++ ;
359 return res ;
360
361 }
362 // Cas ou le calcul a deja ete effectue :
363 else
364 return *tab[indice] ;
365}
366
367
368 //-------------------
369 //-- Fonction ----
370 //-------------------
371
372
373Matrice nondeg_pvect_r(const Matrice &lap, int l, double echelle, int puis, int base_r)
374{
375
376 // Routines de derivation
377 static Matrice (*nondeg_pvect_r[MAX_BASE])(const Matrice&, int, double, int) ;
378 static int nap = 0 ;
379
380 // Premier appel
381 if (nap==0) {
382 nap = 1 ;
383 for (int i=0 ; i<MAX_BASE ; i++) {
384 nondeg_pvect_r[i] = _nondeg_pvect_r_pas_prevu ;
385 }
386 // Les routines existantes
387 nondeg_pvect_r[R_CHEB >> TRA_R] = _nondeg_pvect_r_cheb ;
388 nondeg_pvect_r[R_CHEBU >> TRA_R] = _nondeg_pvect_r_chebu ;
389 nondeg_pvect_r[R_CHEBP >> TRA_R] = _nondeg_pvect_r_chebp ;
390 nondeg_pvect_r[R_CHEBI >> TRA_R] = _nondeg_pvect_r_chebi ;
391 }
392
393 Matrice res(nondeg_pvect_r[base_r](lap, l, echelle, puis)) ;
394 return res ;
395}
396
397}
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