LORENE
cl_pvect_r.C
1 /*
2  * Methods for linear combinations for Ope_pois_vect_r
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 cl_pvect_r_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/cl_pvect_r.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $" ;
29 
30 /*
31  * $Id: cl_pvect_r.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
32  * $Log: cl_pvect_r.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/05/10 15:28:22 j_novak
40  * First version of functions for the solution of the r-component of the
41  * vector Poisson equation.
42  *
43  *
44  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/cl_pvect_r.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 /* FONCTIONS FAISANT DES COMBINAISONS LINEAIRES DES LIGNES.
55  *
56  * Entree :
57  * La Matrice de l'operateur
58  * l : nbre quantique
59  * puis = puissance dans la ZEC
60  * La base de developpement en R
61  *
62  * Sortie :
63  * Renvoie la matrice apres Combinaison lineaire
64  *
65  */
66 
67 namespace Lorene {
68 Matrice _cl_pvect_r_pas_prevu (const Matrice&, int, double, int) ;
69 Matrice _cl_pvect_r_cheb (const Matrice&, int, double, int) ;
70 Matrice _cl_pvect_r_chebi (const Matrice&, int, double, int) ;
71 Matrice _cl_pvect_r_chebu (const Matrice&, int, double, int) ;
72 Matrice _cl_pvect_r_chebp (const Matrice&, int, double, int) ;
73 
74 // Version Matrice --> Matrice
75 Matrice _cl_pvect_r_pas_prevu (const Matrice &source, int l, double echelle, int puis) {
76  cout << "Combinaison lineaire pas prevu..." << endl ;
77  cout << "Source : " << source << endl ;
78  cout << "l : " << l << endl ;
79  cout << "dzpuis : " << puis << endl ;
80  cout << "Echelle : " << echelle << endl ;
81  abort() ;
82  exit(-1) ;
83  return source;
84 }
85 
86 
87  //-------------------
88  //-- R_CHEB ------
89  //-------------------
90 
91 Matrice _cl_pvect_r_cheb (const Matrice &source, int l, double echelle, int) {
92  int n = source.get_dim(0) ;assert (n == source.get_dim(1)) ;
93 
94 
95  const int nmax = 100 ; // Nombre de Matrices stockees
96  static Matrice* tab[nmax] ; // les matrices calculees
97  static int nb_dejafait = 0 ; // nbre de matrices calculees
98  static int l_dejafait[nmax] ;
99  static int nr_dejafait[nmax] ;
100  static double vieux_echelle = 0 ;
101 
102  // Si on a change l'echelle : on detruit tout :
103  if (vieux_echelle != echelle) {
104  for (int i=0 ; i<nb_dejafait ; i++) {
105  l_dejafait[i] = -1 ;
106  nr_dejafait[i] = -1 ;
107  delete tab[i] ;
108  }
109  nb_dejafait = 0 ;
110  vieux_echelle = echelle ;
111  }
112 
113  int indice = -1 ;
114 
115  // On determine si la matrice a deja ete calculee :
116  for (int conte=0 ; conte<nb_dejafait ; conte ++)
117  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
118  indice = conte ;
119 
120  // Calcul a faire :
121  if (indice == -1) {
122  if (nb_dejafait >= nmax) {
123  cout << "_cl_pvect_r_cheb : trop de matrices" << endl ;
124  abort() ;
125  exit (-1) ;
126  }
127 
128  l_dejafait[nb_dejafait] = l ;
129  nr_dejafait[nb_dejafait] = n ;
130 
131  Matrice barre(source) ;
132  int dirac = 1 ;
133  for (int i=0 ; i<n-2 ; i++) {
134  for (int j=i ; j<(n>(i+7)? i+7 : n) ; j++)
135  barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j))
136  /(i+1) ;
137  if (i==0) dirac = 0 ;
138  }
139 
140  Matrice res(barre) ;
141  for (int i=0 ; i<n-4 ; i++)
142  for (int j=i ; j<(n>(i+5)? i+5 : n) ; j++)
143  res.set(i, j) = barre(i, j)-barre(i+2, j) ;
144  tab[nb_dejafait] = new Matrice(res) ;
145  nb_dejafait ++ ;
146  return res ;
147  }
148 
149  // Cas ou le calcul a deja ete effectue :
150  else
151  return *tab[indice] ;
152 }
153 
154  //-------------------
155  //-- R_CHEBP -----
156  //-------------------
157 
158 
159 Matrice _cl_pvect_r_chebp (const Matrice &source, int l, double, int) {
160 
161  int n = source.get_dim(0) ;
162  assert (n == source.get_dim(1)) ;
163 
164  const int nmax = 100 ; // Nombre de Matrices stockees
165  static Matrice* tab[nmax] ; // les matrices calculees
166  static int nb_dejafait = 0 ; // nbre de matrices calculees
167  static int l_dejafait[nmax] ;
168  static int nr_dejafait[nmax] ;
169 
170  int indice = -1 ;
171 
172  // On determine si la matrice a deja ete calculee :
173  for (int conte=0 ; conte<nb_dejafait ; conte ++)
174  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
175  indice = conte ;
176 
177  // Calcul a faire :
178  if (indice == -1) {
179  if (nb_dejafait >= nmax) {
180  cout << "_cl_pvect_r_chebp : trop de matrices" << endl ;
181  abort() ;
182  exit (-1) ;
183  }
184 
185  l_dejafait[nb_dejafait] = l ;
186  nr_dejafait[nb_dejafait] = n ;
187 
188  Matrice barre(source) ;
189 
190  int dirac = 1 ;
191  for (int i=0 ; i<n-2 ; i++) {
192  for (int j=0 ; j<n ; j++)
193  barre.set(i, j) = (1+dirac)*source(i, j)-source(i+2, j) ;
194  if (i==0) dirac = 0 ;
195  }
196 
197  Matrice tilde(barre) ;
198  for (int i=0 ; i<n-4 ; i++)
199  for (int j=0 ; j<n ; j++)
200  tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
201 
202  Matrice res(tilde) ;
203  for (int i=0 ; i<n-4 ; i++)
204  for (int j=0 ; j<n ; j++)
205  res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
206  tab[nb_dejafait] = new Matrice(res) ;
207  nb_dejafait ++ ;
208  return res ;
209  }
210 
211  // Cas ou le calcul a deja ete effectue :
212  else
213  return *tab[indice] ;
214 }
215 
216  //-------------------
217  //-- R_CHEBI -----
218  //-------------------
219 
220 
221 Matrice _cl_pvect_r_chebi (const Matrice &source, int l, double, int) {
222  int n = source.get_dim(0) ;
223  assert (n == source.get_dim(1)) ;
224 
225 
226  const int nmax = 100 ; // Nombre de Matrices stockees
227  static Matrice* tab[nmax] ; // les matrices calculees
228  static int nb_dejafait = 0 ; // nbre de matrices calculees
229  static int l_dejafait[nmax] ;
230  static int nr_dejafait[nmax] ;
231 
232  int indice = -1 ;
233 
234  // On determine si la matrice a deja ete calculee :
235  for (int conte=0 ; conte<nb_dejafait ; conte ++)
236  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
237  indice = conte ;
238 
239  // Calcul a faire :
240  if (indice == -1) {
241  if (nb_dejafait >= nmax) {
242  cout << "_cl_pvect_r_chebi : trop de matrices" << endl ;
243  abort() ;
244  exit (-1) ;
245  }
246 
247  l_dejafait[nb_dejafait] = l ;
248  nr_dejafait[nb_dejafait] = n ;
249 
250  Matrice barre(source) ;
251 
252  for (int i=0 ; i<n-2 ; i++)
253  for (int j=0 ; j<n ; j++)
254  barre.set(i, j) = source(i, j)-source(i+2, j) ;
255 
256  Matrice tilde(barre) ;
257  for (int i=0 ; i<n-4 ; i++)
258  for (int j=0 ; j<n ; j++)
259  tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
260 
261  Matrice res(tilde) ;
262  for (int i=0 ; i<n-4 ; i++)
263  for (int j=0 ; j<n ; j++)
264  res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
265  tab[nb_dejafait] = new Matrice(res) ;
266  nb_dejafait ++ ;
267  return res ;
268  }
269 
270  // Cas ou le calcul a deja ete effectue :
271  else
272  return *tab[indice] ;
273 }
274  //-------------------
275  //-- R_CHEBU -----
276  //-------------------
277 
278 Matrice _cl_pvect_r_chebu (const Matrice &source, int l, double, int puis) {
279  int n = source.get_dim(0) ;
280  assert (n == source.get_dim(1)) ;
281  if (puis != 4) {
282  cout << "_ope_pvect_r_mat_r_chebu : only the case dzpuis = 4 "
283  << '\n' << "is implemented! \n"
284  << "dzpuis = " << puis << endl ;
285  abort() ;
286  }
287 
288  const int nmax = 200 ; // Nombre de Matrices stockees
289  static Matrice* tab[nmax] ; // les matrices calculees
290  static int nb_dejafait = 0 ; // nbre de matrices calculees
291  static int l_dejafait[nmax] ;
292  static int nr_dejafait[nmax] ;
293 
294  int indice = -1 ;
295 
296  // On determine si la matrice a deja ete calculee :
297  for (int conte=0 ; conte<nb_dejafait ; conte ++)
298  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
299  indice = conte ;
300 
301  // Calcul a faire :
302  if (indice == -1) {
303  if (nb_dejafait >= nmax) {
304  cout << "_cl_pvect_r_chebu_quatre : trop de matrices" << endl ;
305  abort() ;
306  exit (-1) ;
307  }
308 
309  l_dejafait[nb_dejafait] = l ;
310  nr_dejafait[nb_dejafait] = n ;
311 
312  Matrice barre(source) ;
313 
314  int dirac = 1 ;
315  for (int i=0 ; i<n-2 ; i++) {
316  for (int j=0 ; j<n ; j++)
317  barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
318  if (i==0) dirac = 0 ;
319  }
320 
321  Matrice tilde(barre) ;
322  for (int i=0 ; i<n-4 ; i++)
323  for (int j=0 ; j<n ; j++)
324  tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
325 
326  Matrice prime(tilde) ;
327  for (int i=0 ; i<n-4 ; i++)
328  for (int j=0 ; j<n ; j++)
329  prime.set(i, j) = (tilde(i, j)-tilde(i+1, j)) ;
330 
331  Matrice res(prime) ;
332  for (int i=0 ; i<n-4 ; i++)
333  for (int j=0 ; j<n ; j++)
334  res.set(i, j) = (prime(i, j)-prime(i+2, j)) ;
335  tab[nb_dejafait] = new Matrice(res) ;
336  nb_dejafait ++ ;
337  return res ;
338  }
339 
340  // Cas ou le calcul a deja ete effectue :
341  else
342  return *tab[indice] ;
343 }
344 
345 
346  //-------------------------
347  //- La routine a appeler ---
348  //---------------------------
349 
350 Matrice cl_pvect_r(const Matrice &source, int l, double echelle,
351  int puis, int base_r) {
352 
353  // Routines de derivation
354  static Matrice (*combinaison[MAX_BASE])(const Matrice &, int, double, int) ;
355  static int nap = 0 ;
356 
357  // Premier appel
358  if (nap==0) {
359  nap = 1 ;
360  for (int i=0 ; i<MAX_BASE ; i++) {
361  combinaison[i] = _cl_pvect_r_pas_prevu ;
362  }
363  // Les routines existantes
364  combinaison[R_CHEB >> TRA_R] = _cl_pvect_r_cheb ;
365  combinaison[R_CHEBU >> TRA_R] = _cl_pvect_r_chebu ;
366  combinaison[R_CHEBP >> TRA_R] = _cl_pvect_r_chebp ;
367  combinaison[R_CHEBI >> TRA_R] = _cl_pvect_r_chebi ;
368  }
369 
370  Matrice res(combinaison[base_r](source, l, echelle, puis)) ;
371  return res ;
372 }
373 
374 }
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Lorene
Lorene prototypes.
Definition: app_hor.h:64
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
Lorene::Matrice::get_dim
int get_dim(int i) const
Returns the dimension of the matrix.
Definition: matrice.C:260
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180