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