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
28char 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
75namespace Lorene {
76Matrice _ope_pvect_r_mat_pas_prevu(int, int, double, int) ;
77Matrice _ope_pvect_r_mat_r_chebp(int, int, double, int) ;
78Matrice _ope_pvect_r_mat_r_chebi(int, int, double, int) ;
79Matrice _ope_pvect_r_mat_r_chebu(int, int, double, int) ;
80Matrice _ope_pvect_r_mat_r_cheb(int, int, double, int) ;
81
82 //-----------------------------------
83 // Routine pour les cas non prevus --
84 //-----------------------------------
85
86Matrice _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
104Matrice _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
193Matrice _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
280Matrice _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
373Matrice _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
502Matrice 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}
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64