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 
28 char 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 
69 namespace Lorene {
70 Matrice _nondeg_ptens_rr_pas_prevu(const Matrice &, int , double, int) ;
71 Matrice _nondeg_ptens_rr_cheb (const Matrice&, int, double, int) ;
72 Matrice _nondeg_ptens_rr_chebp (const Matrice&, int, double, int) ;
73 Matrice _nondeg_ptens_rr_chebi (const Matrice&, int, double, int) ;
74 Matrice _nondeg_ptens_rr_chebu (const Matrice&, int, double, int) ;
75 
76 
77  //------------------------------------
78  // Routine pour les cas non prevus --
79  //-----------------------------------
80 
81 Matrice _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 
99 Matrice _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 
169 Matrice _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 
239 Matrice _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 
297 Matrice _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 
355 Matrice 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.
Definition: type_parite.h:144
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Lorene prototypes.
Definition: app_hor.h:64