LORENE
ope_pvect_r_mat.C
1 /*
2  * Builds the operator 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 ope_pvect_r_mat_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/ope_pvect_r_mat.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $" ;
29 
30 /*
31  * $Id: ope_pvect_r_mat.C,v 1.3 2014/10/13 08:53:34 j_novak Exp $
32  * $Log: ope_pvect_r_mat.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/ope_pvect_r_mat.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 #include "proto.h"
54 
55 /*
56  * Routine caluclant l'operateur suivant :
57  *
58  * -R_CHEB : r^2d2sdx2+2*rdsdx-l*(l+1)Id
59  *
60  * -R_CHEBP et R_CHEBI : d2sdx2+2/r dsdx-l(l+1)/r^2
61  *
62  * -R_CHEBU : d2sdx2-l(l+1)/x^2
63  *
64  * Entree :
65  * -n nbre de points en r
66  * -l voire operateur.
67  * -echelle utile uniquement pour R_CHEB : represente beta/alpha
68  * (cf mapping)
69  * - puis : exposant de multiplication dans la ZEC
70  * - base_r : base de developpement
71  * Sortie :
72  * La fonction renvoie la matrice.
73  */
74 
75 namespace Lorene {
76 Matrice _ope_pvect_r_mat_pas_prevu(int, int, double, int) ;
77 Matrice _ope_pvect_r_mat_r_chebp(int, int, double, int) ;
78 Matrice _ope_pvect_r_mat_r_chebi(int, int, double, int) ;
79 Matrice _ope_pvect_r_mat_r_chebu(int, int, double, int) ;
80 Matrice _ope_pvect_r_mat_r_cheb(int, int, double, int) ;
81 
82  //-----------------------------------
83  // Routine pour les cas non prevus --
84  //-----------------------------------
85 
86 Matrice _ope_pvect_r_mat_pas_prevu(int n, int l, double echelle, int puis) {
87  cout << "laplacien vect_r pas prevu..." << endl ;
88  cout << "n : " << n << endl ;
89  cout << "l : " << l << endl ;
90  cout << "puissance : " << puis << endl ;
91  cout << "echelle : " << echelle << endl ;
92  abort() ;
93  exit(-1) ;
94  Matrice res(1, 1) ;
95  return res;
96 }
97 
98 
99  //-------------------------
100  //---- CAS R_CHEBP ----
101  //-------------------------
102 
103 
104 Matrice _ope_pvect_r_mat_r_chebp (int n, int l, double, int) {
105 
106  const int nmax = 100 ;// Nombre de Matrices stockees
107  static Matrice* tab[nmax] ; // les matrices calculees
108  static int nb_dejafait = 0 ; // nbre de matrices calculees
109  static int l_dejafait[nmax] ;
110  static int nr_dejafait[nmax] ;
111 
112  int indice = -1 ;
113 
114  // On determine si la matrice a deja ete calculee :
115  for (int conte=0 ; conte<nb_dejafait ; conte ++)
116  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
117  indice = conte ;
118 
119  // Calcul a faire :
120  if (indice == -1) {
121  if (nb_dejafait >= nmax) {
122  cout << "_ope_pvect_r_mat_r_chebp : trop de matrices" << endl ;
123  abort() ;
124  exit (-1) ;
125  }
126 
127 
128  l_dejafait[nb_dejafait] = l ;
129  nr_dejafait[nb_dejafait] = n ;
130 
131  Matrice dd(n, n) ;
132  dd.set_etat_qcq() ;
133  Matrice xd(n, n) ;
134  xd.set_etat_qcq() ;
135  Matrice xx(n, n) ;
136  xx.set_etat_qcq() ;
137 
138  double* vect = new double[n] ;
139 
140  for (int i=0 ; i<n ; i++) {
141  for (int j=0 ; j<n ; j++)
142  vect[j] = 0 ;
143  vect[i] = 1 ;
144  d2sdx2_1d (n, &vect, R_CHEBP) ;
145 
146  for (int j=0 ; j<n ; j++)
147  dd.set(j, i) = vect[j] ;
148  }
149 
150  for (int i=0 ; i<n ; i++) {
151  for (int j=0 ; j<n ; j++)
152  vect[j] = 0 ;
153  vect[i] = 1 ;
154  sxdsdx_1d (n, &vect, R_CHEBP) ;
155  for (int j=0 ; j<n ; j++)
156  xd.set(j, i) = vect[j] ;
157 
158  }
159 
160  for (int i=0 ; i<n ; i++) {
161  for (int j=0 ; j<n ; j++)
162  vect[j] = 0 ;
163  vect[i] = 1 ;
164  sx2_1d (n, &vect, R_CHEBP) ;
165  for (int j=0 ; j<n ; j++)
166  xx.set(j, i) = vect[j] ;
167  }
168 
169  delete [] vect ;
170 
171  Matrice res(n, n) ;
172  if (l != 0)
173  res = dd+4*xd+(2-l*(l+1))*xx ;
174  else
175  res = dd + 2*xd - 2*xx ;
176  tab[nb_dejafait] = new Matrice(res) ;
177  nb_dejafait ++ ;
178  return res ;
179  }
180 
181  // Cas ou le calcul a deja ete effectue :
182  else
183  return *tab[indice] ;
184 }
185 
186 
187 
188  //------------------------
189  //-- CAS R_CHEBI ----
190  //------------------------
191 
192 
193 Matrice _ope_pvect_r_mat_r_chebi (int n, int l, double, int) {
194 
195  const int nmax = 100 ;// Nombre de Matrices stockees
196  static Matrice* tab[nmax] ; // les matrices calculees
197  static int nb_dejafait = 0 ; // nbre de matrices calculees
198  static int l_dejafait[nmax] ;
199  static int nr_dejafait[nmax] ;
200 
201  int indice = -1 ;
202 
203  // On determine si la matrice a deja ete calculee :
204  for (int conte=0 ; conte<nb_dejafait ; conte ++)
205  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
206  indice = conte ;
207 
208  // Calcul a faire :
209  if (indice == -1) {
210  if (nb_dejafait >= nmax) {
211  cout << "_ope_pvect_r_mat_r_chebi : trop de matrices" << endl ;
212  abort() ;
213  exit (-1) ;
214  }
215 
216 
217  l_dejafait[nb_dejafait] = l ;
218  nr_dejafait[nb_dejafait] = n ;
219 
220  Matrice dd(n, n) ;
221  dd.set_etat_qcq() ;
222  Matrice xd(n, n) ;
223  xd.set_etat_qcq() ;
224  Matrice xx(n, n) ;
225  xx.set_etat_qcq() ;
226 
227  double* vect = new double[n] ;
228 
229  for (int i=0 ; i<n ; i++) {
230  for (int j=0 ; j<n ; j++)
231  vect[j] = 0 ;
232  vect[i] = 1 ;
233  d2sdx2_1d (n, &vect, R_CHEBI) ; // appel dans le cas impair
234  for (int j=0 ; j<n ; j++)
235  dd.set(j, i) = vect[j] ;
236  }
237 
238  for (int i=0 ; i<n ; i++) {
239  for (int j=0 ; j<n ; j++)
240  vect[j] = 0 ;
241  vect[i] = 1 ;
242  sxdsdx_1d (n, &vect, R_CHEBI) ;
243  for (int j=0 ; j<n ; j++)
244  xd.set(j, i) = vect[j] ;
245  }
246 
247  for (int i=0 ; i<n ; i++) {
248  for (int j=0 ; j<n ; j++)
249  vect[j] = 0 ;
250  vect[i] = 1 ;
251  sx2_1d (n, &vect, R_CHEBI) ;
252  for (int j=0 ; j<n ; j++)
253  xx.set(j, i) = vect[j] ;
254  }
255 
256  delete [] vect ;
257 
258  Matrice res(n, n) ;
259  if (l != 0)
260  res = dd+4*xd+(2-l*(l+1))*xx ;
261  else
262  res = dd + 2*xd - 2*xx ;
263  tab[nb_dejafait] = new Matrice(res) ;
264  nb_dejafait ++ ;
265  return res ;
266  }
267 
268  // Cas ou le calcul a deja ete effectue :
269  else
270  return *tab[indice] ;
271 }
272 
273 
274 
275 
276  //------------------------
277  //---- CAS R_CHEBU ----
278  //------------------------
279 
280 Matrice _ope_pvect_r_mat_r_chebu( int n, int l, double, int puis) {
281 
282  if (puis != 4) {
283  cout << "_ope_pvect_r_mat_r_chebu : only the case dzpuis = 4 "
284  << '\n' << "is implemented! \n"
285  << "dzpuis = " << puis << endl ;
286  abort() ;
287  }
288 
289  const int nmax = 200 ;// Nombre de Matrices stockees
290  static Matrice* tab[nmax] ; // les matrices calculees
291  static int nb_dejafait = 0 ; // nbre de matrices calculees
292  static int l_dejafait[nmax] ;
293  static int nr_dejafait[nmax] ;
294 
295  int indice = -1 ;
296 
297  // On determine si la matrice a deja ete calculee :
298  for (int conte=0 ; conte<nb_dejafait ; conte ++)
299  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
300  indice = conte ;
301 
302  // Calcul a faire :
303  if (indice == -1) {
304  if (nb_dejafait >= nmax) {
305  cout << "_ope_pvect_r_mat_r_chebu : trop de matrices" << endl ;
306  abort() ;
307  exit (-1) ;
308  }
309 
310  l_dejafait[nb_dejafait] = l ;
311  nr_dejafait[nb_dejafait] = n ;
312 
313  Matrice dd(n, n) ;
314  dd.set_etat_qcq() ;
315  Matrice xd(n, n) ;
316  xd.set_etat_qcq() ;
317  Matrice xx(n, n) ;
318  xx.set_etat_qcq() ;
319 
320  double* vect = new double[n] ;
321 
322  for (int i=0 ; i<n ; i++) {
323  for (int j=0 ; j<n ; j++)
324  vect[j] = 0 ;
325  vect[i] = 1 ;
326  d2sdx2_1d (n, &vect, R_CHEBU) ; // appel dans le cas unsurr
327  for (int j=0 ; j<n ; j++)
328  dd.set(j, i) = vect[j] ;
329  }
330 
331  for (int i=0 ; i<n ; i++) {
332  for (int j=0 ; j<n ; j++)
333  vect[j] = 0 ;
334  vect[i] = 1 ;
335  dsdx_1d (n, &vect, R_CHEBU) ; // appel dans le cas unsurr
336  sxm1_1d_cheb (n, vect) ;
337  for (int j=0 ; j<n ; j++)
338  xd.set(j, i) = vect[j] ;
339  }
340 
341  for (int i=0 ; i<n ; i++) {
342  for (int j=0 ; j<n ; j++)
343  vect[j] = 0 ;
344  vect[i] = 1 ;
345  sx2_1d (n, &vect, R_CHEBU) ;
346  for (int j=0 ; j<n ; j++)
347  xx.set(j, i) = vect[j] ;
348  }
349 
350  delete [] vect ;
351 
352  Matrice res(n, n) ;
353  if (l == 0)
354  res = dd-2*xx ;
355  else
356  res = dd - 2*xd + (2 -l*(l+1))*xx ;
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  //---- CAS R_CHEB ----
370  //-----------------------
371 
372 
373 Matrice _ope_pvect_r_mat_r_cheb (int n, int l, double echelle, int) {
374 
375  const int nmax = 200 ;// Nombre de Matrices stockees
376  static Matrice* tab[nmax] ; // les matrices calculees
377  static int nb_dejafait = 0 ; // nbre de matrices calculees
378  static int l_dejafait[nmax] ;
379  static int nr_dejafait[nmax] ;
380  static double vieux_echelle = 0;
381 
382  // Si on a change l'echelle : on detruit tout :
383  if (vieux_echelle != echelle) {
384  for (int i=0 ; i<nb_dejafait ; i++) {
385  l_dejafait[i] = -1 ;
386  nr_dejafait[i] = -1 ;
387  delete tab[i] ;
388  }
389 
390  nb_dejafait = 0 ;
391  vieux_echelle = echelle ;
392  }
393 
394  int indice = -1 ;
395 
396  // On determine si la matrice a deja ete calculee :
397  for (int conte=0 ; conte<nb_dejafait ; conte ++)
398  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
399  indice = conte ;
400 
401  // Calcul a faire :
402  if (indice == -1) {
403  if (nb_dejafait >= nmax) {
404  cout << "_ope_pvect_r_mat_r_cheb : trop de matrices" << endl ;
405  abort() ;
406  exit (-1) ;
407  }
408 
409 
410  l_dejafait[nb_dejafait] = l ;
411  nr_dejafait[nb_dejafait] = n ;
412 
413  Matrice dd(n, n) ;
414  dd.set_etat_qcq() ;
415  Matrice xd(n, n) ;
416  xd.set_etat_qcq() ;
417  Matrice xx(n, n) ;
418  xx.set_etat_qcq() ;
419 
420  double* vect = new double[n] ;
421 
422  for (int i=0 ; i<n ; i++) {
423  for (int j=0 ; j<n ; j++)
424  vect[j] = 0 ;
425  vect[i] = 1 ;
426  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
427  for (int j=0 ; j<n ; j++)
428  dd.set(j, i) = vect[j]*echelle*echelle ;
429  }
430 
431  for (int i=0 ; i<n ; i++) {
432  for (int j=0 ; j<n ; j++)
433  vect[j] = 0 ;
434  vect[i] = 1 ;
435  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
436  multx_1d (n, &vect, R_CHEB) ;
437  for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
438  dd.set(j, i) += 2*echelle*vect[j] ;
439  }
440 
441  for (int i=0 ; i<n ; i++) {
442  for (int j=0 ; j<n ; j++)
443  vect[j] = 0 ;
444  vect[i] = 1 ;
445  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
446  multx_1d (n, &vect, R_CHEB) ;
447  multx_1d (n, &vect, R_CHEB) ;
448  for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
449  dd.set(j, i) += vect[j] ;
450  }
451 
452  for (int i=0 ; i<n ; i++) {
453  for (int j=0 ; j<n ; j++)
454  vect[j] = 0 ;
455  vect[i] = 1 ;
456  sxdsdx_1d (n, &vect, R_CHEB) ;
457  for (int j=0 ; j<n ; j++)
458  xd.set(j, i) = vect[j]*echelle ;
459  }
460 
461  for (int i=0 ; i<n ; i++) {
462  for (int j=0 ; j<n ; j++)
463  vect[j] = 0 ;
464  vect[i] = 1 ;
465  sxdsdx_1d (n, &vect, R_CHEB) ;
466  multx_1d (n, &vect, R_CHEB) ;
467  for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
468  xd.set(j, i) += vect[j] ;
469  }
470 
471  for (int i=0 ; i<n ; i++) {
472  for (int j=0 ; j<n ; j++)
473  vect[j] = 0 ;
474  vect[i] = 1 ;
475  sx2_1d (n, &vect, R_CHEB) ;
476  for (int j=0 ; j<n ; j++)
477  xx.set(j, i) = vect[j] ;
478  }
479 
480  delete [] vect ;
481 
482  Matrice res(n, n) ;
483  if (l == 0)
484  res = dd+2*xd-2*xx ;
485  else
486  res = dd + 4*xd + (2 - l*(l+1))*xx ;
487  tab[nb_dejafait] = new Matrice(res) ;
488  nb_dejafait ++ ;
489  return res ;
490  }
491 
492  // Cas ou le calcul a deja ete effectue :
493  else
494  return *tab[indice] ;
495 }
496 
497 
498  //----------------------------
499  //--- La routine a appeler ---
500  //----------------------------
501 
502 Matrice ope_pvect_r_mat(int n, int l, double echelle, int puis, int base_r)
503 {
504 
505  // Routines de derivation
506  static Matrice (*ope_pvect_r_mat[MAX_BASE])(int, int, double, int) ;
507  static int nap = 0 ;
508 
509  // Premier appel
510  if (nap==0) {
511  nap = 1 ;
512  for (int i=0 ; i<MAX_BASE ; i++) {
513  ope_pvect_r_mat[i] = _ope_pvect_r_mat_pas_prevu ;
514  }
515  // Les routines existantes
516  ope_pvect_r_mat[R_CHEB >> TRA_R] = _ope_pvect_r_mat_r_cheb ;
517  ope_pvect_r_mat[R_CHEBU >> TRA_R] = _ope_pvect_r_mat_r_chebu ;
518  ope_pvect_r_mat[R_CHEBP >> TRA_R] = _ope_pvect_r_mat_r_chebp ;
519  ope_pvect_r_mat[R_CHEBI >> TRA_R] = _ope_pvect_r_mat_r_chebi ;
520  }
521 
522  Matrice res(ope_pvect_r_mat[base_r](n, l, echelle, puis)) ;
523  return res ;
524 }
525 
526 }
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
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