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 char 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 
73 namespace Lorene {
74 Matrice _nondeg_pvect_r_pas_prevu(const Matrice &, int , double, int) ;
75 Matrice _nondeg_pvect_r_cheb (const Matrice&, int, double, int) ;
76 Matrice _nondeg_pvect_r_chebp (const Matrice&, int, double, int) ;
77 Matrice _nondeg_pvect_r_chebi (const Matrice&, int, double, int) ;
78 Matrice _nondeg_pvect_r_chebu (const Matrice&, int, double, int) ;
79 
80 
81  //------------------------------------
82  // Routine pour les cas non prevus --
83  //-----------------------------------
84 
85 Matrice _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 
103 Matrice _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 
173 Matrice _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 
243 Matrice _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 
315 Matrice _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 
373 Matrice 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.
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