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