LORENE
dal_inverse.C
1/*
2 * Copyright (c) 2000-2001 Jerome Novak
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23
24
25/*
26 * $Id: dal_inverse.C,v 1.10 2016/12/05 16:18:09 j_novak Exp $
27 * $Log: dal_inverse.C,v $
28 * Revision 1.10 2016/12/05 16:18:09 j_novak
29 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30 *
31 * Revision 1.9 2014/10/13 08:53:28 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.8 2014/10/06 15:16:08 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.7 2008/08/27 08:51:15 jl_cornou
38 * Added Jacobi(0,2) polynomials
39 *
40 * Revision 1.6 2008/02/18 13:53:43 j_novak
41 * Removal of special indentation instructions.
42 *
43 * Revision 1.5 2004/10/05 15:44:21 j_novak
44 * Minor speed enhancements.
45 *
46 * Revision 1.4 2002/01/03 15:30:28 j_novak
47 * Some comments modified.
48 *
49 * Revision 1.3 2002/01/03 13:18:41 j_novak
50 * Optimization: the members set(i,j) and operator(i,j) of class Matrice are
51 * now defined inline. Matrice is a friend class of Tbl.
52 *
53 * Revision 1.2 2002/01/02 14:07:57 j_novak
54 * Dalembert equation is now solved in the shells. However, the number of
55 * points in theta and phi must be the same in each domain. The solver is not
56 * completely tested (beta version!).
57 *
58 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
59 * LORENE
60 *
61 * Revision 1.1 2000/12/04 16:37:03 novak
62 * Initial revision
63 *
64 *
65 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/dal_inverse.C,v 1.10 2016/12/05 16:18:09 j_novak Exp $
66 *
67 */
68
69//fichiers includes
70#include <cmath>
71
72//Headers LORENE
73#include "param.h"
74#include "matrice.h"
75#include "proto.h"
76
77/***************************************************************************
78 *
79 * Set of routines to invert the dalembertian operator given
80 * by get_operateur. The matrix is put into a banded form, following
81 * the type of the operator (type_dal) and the odd/even decomposition
82 * base (base_r).
83 * type_dal is supposed to be given by get_operateur (see the various cases
84 * there and in type_parite.h).
85 * part is a boolean saying wether one is looking for a particular sol.
86 * of the system defined by operateur (i.e. X such as operateur*X = source),
87 * part = true, or a homogeneous one (part = false).
88 *
89 ***************************************************************************/
90
91 //------------------------------------
92 // Routine pour les cas non prevus --
93 //------------------------------------
94namespace Lorene {
95Tbl _dal_inverse_pas_prevu (const Matrice&, const Tbl&, const bool) {
96 cout << " Inversion du dalembertien pas prevue ..... : "<< endl ;
97 abort() ;
98 exit(-1) ;
99 Tbl res(1) ;
100 return res;
101}
102
103
104 //-------------------
105 //-- R_CHEB -----
106 //-------------------
107
108Tbl _dal_inverse_r_cheb_o2d_s(const Matrice &op, const Tbl &source,
109 const bool part) {
110
111 // Operator and source are copied and prepared
112 Matrice barre(op) ;
113 int nr = op.get_dim(0) ;
114 Tbl aux(source) ;
115
116 // Operator is put into banded form (changing the image base)
117
118 int dirac = 2 ; // Don't forget the factor 2 for T_0!!
119 for (int i=0; i<nr-4; i++) {
120 int nrmin = (i>1 ? i-1 : 0) ;
121 int nrmax = (i<nr-9 ? i+10 : nr) ;
122 for (int j=nrmin; j<nrmax; j++)
123 barre.set(i,j) = (op(i+2,j) - dirac*op(i,j))/(i+1) ;
124 if (part)
125 aux.set(i) = (source(i+2) - dirac*source(i))/(i+1) ;
126 if (i==0) dirac = 1 ;
127 }
128 for (int i=0; i<nr-4; i++) {
129 int nrmin = (i>1 ? i-1 : 0) ;
130 int nrmax = (i<nr-9 ? i+10 : nr) ;
131 for (int j=nrmin; j<nrmax; j++)
132 barre.set(i,j) = barre(i+2,j) - barre(i,j) ;
133 if (part) aux.set(i) = aux(i+2) - aux(i) ;
134 }
135
136 // 6 over-diagonals and 2 under...
137 barre.set_band(5,1) ;
138
139 // LU decomposition
140 barre.set_lu() ;
141
142 // Inversion using LAPACK
143
144 return barre.inverse(aux) ;
145}
146
147Tbl _dal_inverse_r_cheb_o2d_l(const Matrice &op, const Tbl &source,
148 const bool part) {
149
150 // Operator and source are copied and prepared
151 Matrice barre(op) ;
152 int nr = op.get_dim(0) ;
153 Tbl aux(source) ;
154
155 // Operator is put into banded form (changing the image base)
156
157 int dirac = 2 ; // Don't forget the factor 2 for T_0!!
158 for (int i=0; i<nr-4; i++) {
159 int nrmin = (i>1 ? i-1 : 0) ;
160 int nrmax = (i<nr-9 ? i+10 : nr) ;
161 for (int j=nrmin; j<nrmax; j++)
162 barre.set(i,j) = (op(i+2,j) - dirac*op(i,j))/(i+1) ;
163 if (part)
164 aux.set(i) = (source(i+2) - dirac*source(i))/(i+1) ;
165 if (i==0) dirac = 1 ;
166 }
167 for (int i=0; i<nr-4; i++) {
168 int nrmin = (i>1 ? i-1 : 0) ;
169 int nrmax = (i<nr-9 ? i+10 : nr) ;
170 for (int j=nrmin; j<nrmax; j++)
171 barre.set(i,j) = barre(i+2,j) - barre(i,j) ;
172 if (part) aux.set(i) = aux(i+2) - aux(i) ;
173 }
174
175 // In this case the time-step is too large for the number of points
176 // (or the number of points too large for the time-step!)
177 // the matrix is ill-conditionned: the last lines are put as first
178 // ones and all the others are shifted.
179
180 double temp1, temp2 ;
181 temp1 = aux(nr-1) ;
182 temp2 = aux(nr-2) ;
183 for (int i=nr-3; i>=0; i--) {
184 int nrmin = (i>1 ? i-1 : 0) ;
185 int nrmax = (i<nr-9 ? i+10 : nr) ;
186 for (int j=nrmin; j<nrmax; j++)
187 barre.set(i+2,j) = barre(i,j) ;
188 aux.set(i+2) = aux(i) ;
189 }
190 aux.set(0) = temp2 ;
191 aux.set(1) = temp1 ;
192
193 barre.set(0,0) = 0.5 ;
194 barre.set(0,1) = 1. ;
195 barre.set(0,2) = 1. ;
196 barre.set(0,3) = 1. ;
197 barre.set(0,4) = 0. ;
198
199 barre.set(1,0) = 0. ;
200 barre.set(1,1) = 0.5 ;
201 barre.set(1,2) = 1. ;
202 barre.set(1,3) = -1. ;
203 barre.set(1,4) = 1. ;
204 barre.set(1,5) = 0. ;
205
206 // 3 over diagonals and 3 under ...
207 barre.set_band(3,3) ;
208
209 // LU decomposition
210 barre.set_lu() ;
211
212 // Inversion using LAPACK
213
214 return barre.inverse(aux) ;
215}
216
217Tbl _dal_inverse_r_cheb_o2_s(const Matrice &op, const Tbl &source,
218 const bool part) {
219
220 // Operator and source are copied and prepared
221 Matrice barre(op) ;
222 int nr = op.get_dim(0) ;
223 Tbl aux(source) ;
224
225 // Operator is put into banded form (changing the image base)
226
227 int dirac = 2 ; // Don't forget the factor 2 for T_0!!
228 for (int i=0; i<nr-4; i++) {
229 int nrmin = (i>2 ? i-2 : 0) ;
230 int nrmax = (i<nr-10 ? i+11 : nr) ;
231 for (int j=nrmin; j<nrmax; j++)
232 barre.set(i,j) = (op(i+2,j) - dirac*op(i,j))/(i+1) ;
233 if (part)
234 aux.set(i) = (source(i+2) - dirac*source(i))/(i+1) ;
235 if (i==0) dirac = 1 ;
236 }
237 for (int i=0; i<nr-4; i++) {
238 int nrmin = (i>2 ? i-2 : 0) ;
239 int nrmax = (i<nr-10 ? i+11 : nr) ;
240 for (int j=nrmin; j<nrmax; j++)
241 barre.set(i,j) = barre(i+2,j) - barre(i,j) ;
242 if (part) aux.set(i) = aux(i+2) - aux(i) ;
243 }
244
245 // 6 over-diagonals and 2 under...
246 barre.set_band(6,2) ;
247
248 // LU decomposition
249 barre.set_lu() ;
250
251 // Inversion using LAPACK
252
253 return barre.inverse(aux) ;
254}
255
256Tbl _dal_inverse_r_cheb_o2_l(const Matrice &op, const Tbl &source,
257 const bool part) {
258
259 // Operator and source are copied and prepared
260 Matrice barre(op) ;
261 int nr = op.get_dim(0) ;
262 Tbl aux(source) ;
263
264 // Operator is put into banded form (changing the image base)
265
266 int dirac = 2 ; // Don't forget the factor 2 for T_0!!
267 for (int i=0; i<nr-4; i++) {
268 int nrmin = (i>2 ? i-2 : 0) ;
269 int nrmax = (i<nr-10 ? i+11 : nr) ;
270 for (int j=nrmin; j<nrmax; j++)
271 barre.set(i,j) = (op(i+2,j) - dirac*op(i,j))/(i+1) ;
272 if (part)
273 aux.set(i) = (source(i+2) - dirac*source(i))/(i+1) ;
274 if (i==0) dirac = 1 ;
275 }
276 for (int i=0; i<nr-4; i++) {
277 int nrmin = (i>2 ? i-2 : 0) ;
278 int nrmax = (i<nr-10 ? i+11 : nr) ;
279 for (int j=nrmin; j<nrmax; j++)
280 barre.set(i,j) = barre(i+2,j) - barre(i,j) ;
281 if (part) aux.set(i) = aux(i+2) - aux(i) ;
282 }
283
284 // In this case the time-step is too large for the number of points
285 // (or the number of points too large for the time-step!)
286 // the matrix is ill-conditionned: the last lines are put as first
287 // ones and all the others are shifted.
288
289 double temp1, temp2 ;
290 temp1 = aux(nr-1) ;
291 temp2 = aux(nr-2) ;
292 for (int i=nr-3; i>=0; i--) {
293 int nrmin = (i>2 ? i-2 : 0) ;
294 int nrmax = (i<nr-10 ? i+11 : nr) ;
295 for (int j=nrmin; j<nrmax; j++)
296 barre.set(i+2,j) = barre(i,j) ;
297 aux.set(i+2) = aux(i) ;
298 }
299 aux.set(0) = temp2 ;
300 aux.set(1) = temp1 ;
301
302 barre.set(0,0) = 0.5 ;
303 barre.set(0,1) = 1. ;
304 barre.set(0,2) = 1. ;
305 barre.set(0,3) = 1. ;
306 barre.set(0,4) = 1. ;
307 barre.set(0,5) = 0. ;
308
309 barre.set(1,0) = 0. ;
310 barre.set(1,1) = 0.5 ;
311 barre.set(1,2) = -1. ;
312 barre.set(1,3) = 1. ;
313 barre.set(1,4) = -1. ;
314 barre.set(1,5) = 1. ;
315 barre.set(1,6) = 0. ;
316
317 // 4 over diagonals and 4 under ...
318 barre.set_band(4,4) ;
319
320 // LU decomposition
321 barre.set_lu() ;
322
323 // Inversion using LAPACK
324
325 return barre.inverse(aux) ;
326}
327
328 //-------------------
329 //-- R_CHEBP -----
330 //-------------------
331
332Tbl _dal_inverse_r_chebp_o2d_s(const Matrice &op, const Tbl &source,
333 const bool part) {
334
335 // Operator and source are copied and prepared
336 Matrice barre(op) ;
337 int nr = op.get_dim(0) ;
338 Tbl aux(nr) ;
339 if (part) {
340 aux = source ;
341 aux.set(nr-1) = 0. ;
342 }
343 else {
344 aux.annule_hard() ;
345 aux.set(nr-1) = 1. ;
346 }
347
348 // Operator is put into banded form (changing the image base)
349
350 int dirac = 2 ; // Don't forget the factor 2 for T_0!!
351 for (int i=0; i<nr-4; i++) {
352 int nrmax = (i<nr-7 ? i+8 : nr) ;
353 for (int j=i; j<nrmax; j++)
354 barre.set(i,j) = (op(i+2,j) - dirac*op(i,j))/(i+1) ;
355 if (part)
356 aux.set(i) = (source(i+2) - dirac*source(i))/(i+1) ;
357 if (i==0) dirac = 1 ;
358 }
359 for (int i=0; i<nr-4; i++) {
360 int nrmax = (i<nr-7 ? i+8 : nr) ;
361 for (int j=i; j<nrmax; j++) barre.set(i,j) = barre(i+1,j) - barre(i,j) ;
362 if (part) aux.set(i) = aux(i+1) - aux(i) ;
363 }
364 if (fabs(barre(nr-5,nr-1)) >= 1.e-16) {
365 if (fabs(barre(nr-5,nr-1)) > fabs(barre(nr-2,nr-1))) {
366 double lambda = barre(nr-2,nr-1)/barre(nr-5,nr-1) ;
367 for (int j=nr-5; j<nr; j++) barre.set(nr-5,j) = barre(nr-5,j)*lambda
368 -barre(nr-2,j) ;
369 if (part) aux.set(nr-5) = aux(nr-5)*lambda - aux(nr-2) ;
370 }
371 else {
372 double lambda = barre(nr-5,nr-1)/barre(nr-2,nr-1) ;
373 for (int j=nr-5; j<nr; j++) barre.set(nr-5,j) -= lambda*barre(nr-2,j) ;
374 if (part) aux.set(nr-5) -= lambda*aux(nr-2) ;
375 }
376 }
377
378 // 3 over-diagonals and 0 under...
379 barre.set_band(3,0) ;
380
381 // LU decomposition
382 barre.set_lu() ;
383
384 // Inversion using LAPACK
385
386 return barre.inverse(aux) ;
387}
388
389
390
391Tbl _dal_inverse_r_chebp_o2d_l(const Matrice &op, const Tbl &source,
392 const bool part) {
393
394 // Operator and source are copied and prepared
395 Matrice barre(op) ;
396 int nr = op.get_dim(0) ;
397 Tbl aux(nr) ;
398 if (part) {
399 aux = source ;
400 aux.set(nr-1) = 0. ;
401 }
402 else {
403 aux.annule_hard() ;
404 aux.set(0) = 1. ;
405 }
406
407 // Operator is put into banded form (changing the image base)
408
409 int dirac = 2 ; // Don't forget the factor 2 for T_0!!
410 for (int i=0; i<nr-4; i++) {
411 int nrmax = (i<nr-7 ? i+8 : nr) ;
412 for (int j=i; j<nrmax; j++)
413 barre.set(i,j) = (op(i+2,j) - dirac*op(i,j))/(i+1) ;
414 if (part)
415 aux.set(i) = (source(i+2) - dirac*source(i))/(i+1) ;
416 if (i==0) dirac = 1 ;
417 }
418 for (int i=0; i<nr-4; i++) {
419 int nrmax = (i<nr-7 ? i+8 : nr) ;
420 for (int j=i; j<nrmax; j++) barre.set(i,j) = barre(i+1,j) - barre(i,j) ;
421 if (part) aux.set(i) = aux(i+1) - aux(i) ;
422 }
423 if (fabs(barre(nr-5,nr-1)) >= 1.e-16) {
424 if (fabs(barre(nr-5,nr-1)) > fabs(barre(nr-2,nr-1))) {
425 double lambda = barre(nr-2,nr-1)/barre(nr-5,nr-1) ;
426 for (int j=nr-5; j<nr; j++) barre.set(nr-5,j) = barre(nr-5,j)*lambda
427 -barre(nr-2,j) ;
428 if (part) aux.set(nr-5) = aux(nr-5)*lambda - aux(nr-2) ;
429 }
430 else {
431 double lambda = barre(nr-5,nr-1)/barre(nr-2,nr-1) ;
432 for (int j=nr-5; j<nr; j++) barre.set(nr-5,j) -= lambda*barre(nr-2,j) ;
433 if (part) aux.set(nr-5) -= lambda*aux(nr-2) ;
434 }
435 }
436
437 // In this case the time-step is too large for the number of points
438 // (or the number of points too large for the time-step!)
439 // the matrix is ill-conditionned: the last line is put as the first
440 // one and all the others are shifted.
441
442 for (int i=nr-2; i>=0; i--) {
443 for (int j=i; j<((i+5 > nr) ? nr : i+5); j++)
444 barre.set(i+1,j) = barre(i,j) ;
445 if (part) aux.set(i+1) = aux(i) ;
446 }
447 barre.set(0,0) = 0.5 ;
448 barre.set(0,1) = 1. ;
449 barre.set(0,2) = 1. ;
450 barre.set(0,3) = 0. ;
451
452 if (part) aux.set(0) = 0 ;
453
454 // 2 over diagonals and one under ...
455 barre.set_band(2,1) ;
456
457 // LU decomposition
458 barre.set_lu() ;
459
460 // Inversion using LAPACK
461
462 return barre.inverse(aux);
463}
464
465
466
467Tbl _dal_inverse_r_chebp_o2_s(const Matrice &op, const Tbl &source,
468 const bool part) {
469
470 // Operator and source are copied and prepared
471 Matrice barre(op) ;
472 int nr = op.get_dim(0) ;
473 Tbl aux(nr-1) ;
474 if (part) {
475 aux.set_etat_qcq() ;
476 for (int i=nr-4; i<nr-1; i++) aux.set(i) = source(i) ;
477 }
478 else {
479 aux.annule_hard() ;
480 aux.set(nr-2) = 1. ;
481 }
482
483 // Operator is put into banded form (changing the image base ...)
484
485 int dirac = 2 ;// Don't forget the factor 2 for T_0!!
486 for (int i=0; i<nr-4; i++) {
487 int nrmax = (i<nr-7 ? i+8 : nr) ;
488 for (int j=i; j<nrmax; j++)
489 barre.set(i,j) = (op(i+2,j) - dirac*op(i,j))/(i+1) ;
490 if (part)
491 aux.set(i) = (source(i+2) - dirac*source(i))/(i+1) ;
492 if (i==0) dirac = 1 ;
493 }
494 for (int i=0; i<nr-4; i++) {
495 int nrmax = (i<nr-7 ? i+8 : nr) ;
496 for (int j=i; j<nrmax; j++) barre.set(i,j) = barre(i+1,j) - barre(i,j) ;
497 if (part) aux.set(i) = aux(i+1) - aux(i) ;
498 }
499
500 // ... and changing the starting base (the last column is quit) ...
501
502 Matrice tilde(nr-1,nr-1) ;
503 tilde.set_etat_qcq() ;
504 for (int i=0; i<nr-1; i++)
505 for (int j=0; j<nr-1;j++)
506 tilde.set(i,j) = barre(i,j+1) + barre(i,j) ;
507
508 // 3 over-diagonals and 1 under...
509 tilde.set_band(3,1) ;
510
511 // LU decomposition
512 tilde.set_lu() ;
513
514 // Inversion using LAPACK
515 Tbl res0(tilde.inverse(aux)) ;
516 Tbl res(nr) ;
517 res.set_etat_qcq() ;
518
519 // ... finally, one has to recover the original starting base
520 res.set(0) = res0(0) ;
521 res.set(nr-1) = res0(nr-2) ;
522 for (int i=1; i<nr-1; i++) res.set(i) = res0(i-1) + res0(i);
523
524 return res ;
525
526
527}
528
529Tbl _dal_inverse_r_chebp_o2_l(const Matrice &op, const Tbl &source,
530 const bool part) {
531
532 // Operator and source are copied and prepared
533 Matrice barre(op) ;
534 int nr = op.get_dim(0) ;
535 Tbl aux(nr-1) ;
536 if (part) {
537 aux.set_etat_qcq() ;
538 for (int i=nr-4; i<nr-1; i++) aux.set(i) = source(i) ;
539 }
540 else {
541 aux.annule_hard() ;
542 aux.set(0) = 1 ;
543 }
544
545 // Operator is put into banded form (changing the image base ...)
546
547 int dirac = 2 ;// Don't forget the factor 2 for T_0!!
548 for (int i=0; i<nr-4; i++) {
549 int nrmax = (i<nr-7 ? i+8 : nr) ;
550 for (int j=i; j<nrmax; j++) {
551 barre.set(i,j) = (op(i+2,j) - dirac*op(i,j))/(i+1) ;
552 }
553 if (part)
554 aux.set(i) = (source(i+2) - dirac*source(i))/(i+1) ;
555 if (i==0) dirac = 1 ;
556 }
557 for (int i=0; i<nr-4; i++) {
558 int nrmax = (i<nr-7 ? i+8 : nr) ;
559 for (int j=i; j<nrmax; j++) barre.set(i,j) = barre(i+1,j) - barre(i,j) ;
560 if (part) aux.set(i) = aux(i+1) - aux(i) ;
561 }
562
563 // ... and changing the starting base (the last column is quit) ...
564
565 Matrice tilde(nr-1,nr-1) ;
566 tilde.set_etat_qcq() ;
567 for (int i=0; i<nr-1; i++)
568 for (int j=0; j<nr-1;j++)
569 tilde.set(i,j) = barre(i,j+1) + barre(i,j) ;
570
571 // In this case the time-step is too large for the number of points
572 // (or the number of points too large for the time-step!)
573 // the matrix is ill-conditionned: the last line is put as the first
574 // one and all the others are shifted.
575
576 for (int i=nr-3; i>=0; i--) {
577 for (int j=((i>0) ? i-1 : 0); j<((i+5 > nr-1) ? nr-1 : i+5); j++)
578 tilde.set(i+1,j) = tilde(i,j) ;
579 if (part) aux.set(i+1) = aux(i) ;
580 }
581 tilde.set(0,0) = 0.5 ;
582 tilde.set(0,1) = 1. ;
583 tilde.set(0,2) = 1. ;
584 tilde.set(0,3) = 0. ;
585
586 if (part) aux.set(0) = 0 ;
587
588 // 2 over-diagonals and 2 under...
589 tilde.set_band(2,2) ;
590
591 // LU decomposition
592 tilde.set_lu() ;
593
594 // Inversion using LAPACK
595 Tbl res0(tilde.inverse(aux)) ;
596 Tbl res(nr) ;
597 res.set_etat_qcq() ;
598
599 // ... finally, one has to recover the original starting base
600
601 res.set(0) = res0(0) ;
602 res.set(nr-1) = res0(nr-2) ;
603 for (int i=1; i<nr-1; i++) res.set(i) = res0(i-1) + res0(i);
604
605 return res ;
606
607
608}
609 //-------------------
610 //-- R_CHEBI -----
611 //-------------------
612
613Tbl _dal_inverse_r_chebi_o2d_s(const Matrice &op, const Tbl &source,
614 const bool part) {
615
616 // Operator and source are copied and prepared
617 int nr = op.get_dim(0) ;
618 Matrice barre(nr-1,nr-1) ;
619 barre.set_etat_qcq() ;
620 for (int i=0; i<nr-1; i++)
621 for (int j=0; j<nr-1; j++)
622 barre.set(i,j) = op(i,j) ;
623 Tbl aux(nr-1) ;
624 if (part) {
625 aux.set_etat_qcq() ;
626 for (int i=nr-4; i<nr-1; i++) aux.set(i) = source(i) ;
627 }
628 else {
629 aux.annule_hard() ;
630 aux.set(nr-2) = 1 ;
631 }
632
633 // Operator is put into banded form (changing the image base )
634 // Last column is quit.
635
636 for (int i=0; i<nr-4; i++) {
637 for (int j=i; j<nr-1; j++) {
638 barre.set(i,j) = (op(i+1,j) - op(i,j))/(i+1) ;
639 }
640 if (part) aux.set(i) = (source(i+1) - source(i))/(i+1) ;
641 }
642 for (int i=0; i<nr-5; i++) {
643 for (int j=i; j<nr-1; j++) barre.set(i,j) = barre(i+2,j) - barre(i,j) ;
644 if (part) aux.set(i) = aux(i+2) - aux(i) ;
645 }
646 if (fabs(barre(nr-6,nr-2)) >= 1.e-16) {
647 if (fabs(barre(nr-6,nr-2)) > fabs(barre(nr-3,nr-2))) {
648 double lambda = barre(nr-3,nr-2)/barre(nr-6,nr-2) ;
649 for (int j=0; j<nr-1; j++) barre.set(nr-6,j) = barre(nr-6,j)*lambda
650 -barre(nr-3,j) ;
651 if (part) aux.set(nr-6) = aux(nr-6)*lambda - aux(nr-3) ;
652 }
653 else {
654 double lambda = barre(nr-6,nr-2)/barre(nr-3,nr-2) ;
655 for (int j=0; j<nr-1; j++) barre.set(nr-6,j) -= lambda*barre(nr-3,j) ;
656 if (part) aux.set(nr-6) -= lambda*aux(nr-3) ;
657 }
658 }
659
660 // 3 over-diagonals and 0 under...
661 barre.set_band(3,0) ;
662
663 // LU decomposition
664 barre.set_lu() ;
665
666 // Inversion using LAPACK
667 Tbl res0(barre.inverse(aux)) ;
668 Tbl res(nr) ;
669 res.set_etat_qcq() ;
670 for (int i=0; i<nr-1; i++) res.set(i) = res0(i) ;
671 res.set(nr-1) = 0 ;
672
673 return res ;
674}
675
676Tbl _dal_inverse_r_chebi_o2d_l(const Matrice &op, const Tbl &source,
677 const bool part) {
678
679 // Operator and source are copied and prepared
680 Matrice barre(op) ;
681 int nr = op.get_dim(0) ;
682 Tbl aux(nr-1) ;
683 if (part) {
684 aux.set_etat_qcq() ;
685 for (int i=0; i<nr-2; i++) aux.set(i) = source(i) ;
686 aux.set(nr-2) = 0 ;
687 }
688 else {
689 aux.annule_hard() ;
690 aux.set(0) = 1 ;
691 }
692 // Operator is put into banded form (changing the image base)
693 // Last column is quit.
694
695 for (int i=0; i<nr-4; i++) {
696 for (int j=i; j<nr-1; j++) {
697 barre.set(i,j) = (op(i+1,j) - op(i,j))/(i+1) ;
698 }
699 if (part) aux.set(i) = (source(i+1) - source(i))/(i+1) ;
700 }
701 for (int i=0; i<nr-5; i++) {
702 for (int j=i; j<nr-1; j++) barre.set(i,j) = barre(i+2,j) - barre(i,j) ;
703 if (part) aux.set(i) = aux(i+2) - aux(i) ;
704 }
705 if (fabs(barre(nr-6,nr-2)) >= 1.e-16) {
706 if (fabs(barre(nr-6,nr-2)) > fabs(barre(nr-3,nr-2))) {
707 double lambda = barre(nr-3,nr-2)/barre(nr-6,nr-2) ;
708 for (int j=0; j<nr-1; j++) barre.set(nr-6,j) = barre(nr-6,j)*lambda
709 -barre(nr-3,j) ;
710 if (part) aux.set(nr-6) = aux(nr-6)*lambda - aux(nr-3) ;
711 }
712 else {
713 double lambda = barre(nr-6,nr-2)/barre(nr-3,nr-2) ;
714 for (int j=0; j<nr-1; j++) barre.set(nr-6,j) -= lambda*barre(nr-3,j) ;
715 aux.set(nr-6) -= lambda*aux(nr-3) ;
716 }
717 }
718
719 // In this case the time-step is too large for the number of points
720 // (or the number of points too large for the time-step!)
721 // the matrix is ill-conditionned: the last line is put as the first
722 // one and all the others are shifted.
723
724 Matrice tilde(nr-1,nr-1) ;
725 tilde.set_etat_qcq() ;
726 for (int i=nr-3; i>=0; i--) {
727 for (int j=0; j<nr-1; j++)
728 tilde.set(i+1,j) = barre(i,j) ;
729 if (part) aux.set(i+1) = aux(i) ;
730 }
731 tilde.set(0,0) = 1. ;
732 tilde.set(0,1) = 1. ;
733 tilde.set(0,2) = 1. ;
734 tilde.set(0,3) = 0. ;
735
736 if (part) aux.set(0) = 0 ;
737
738
739 // 2 over-diagonals and 1 under...
740 tilde.set_band(2,1) ;
741
742 // LU decomposition
743 tilde.set_lu() ;
744
745 // Inversion using LAPACK
746 Tbl res0(tilde.inverse(aux)) ;
747 Tbl res(nr) ;
748 res.set_etat_qcq() ;
749 for (int i=0; i<nr-1; i++) res.set(i) = res0(i) ;
750 res.set(nr-1) = 0 ;
751
752 return res ;
753}
754
755Tbl _dal_inverse_r_chebi_o2_s(const Matrice &op, const Tbl &source,
756 const bool part) {
757
758 // Operator and source are copied and prepared
759 Matrice barre(op) ;
760 int nr = op.get_dim(0) ;
761 Tbl aux(nr-2) ;
762 if (part) {
763 aux.set_etat_qcq() ;
764 aux.set(nr-4) = source(nr-4) ;
765 aux.set(nr-3) = 0 ;
766 }
767 else {
768 aux.annule_hard() ;
769 aux.set(nr-3) = 1. ;
770 }
771
772 // Operator is put into banded form (changing the image base ...)
773
774 for (int i=0; i<nr-4; i++) {
775 for (int j=i; j<nr; j++) {
776 barre.set(i,j) = (op(i+1,j) - op(i,j))/(i+1) ;
777 }
778 if (part)
779 aux.set(i) = (source(i+1) - source(i))/(i+1) ;
780 }
781 for (int i=0; i<nr-5; i++) {
782 for (int j=i; j<nr; j++) barre.set(i,j) = barre(i+2,j) - barre(i,j) ;
783 if (part) aux.set(i) = aux(i+2) - aux(i) ;
784 }
785
786 // ... and changing the starting base (first and last columns are quit)...
787
788 Matrice tilde(nr-2,nr-2) ;
789 tilde.set_etat_qcq() ;
790 for (int i=0; i<nr-2; i++)
791 for (int j=0; j<nr-2;j++)
792 tilde.set(i,j) = barre(i,j+1)*(2*j+1) + barre(i,j)*(2*j+3) ;
793
794 // 3 over-diagonals and 1 under...
795 tilde.set_band(3,1) ;
796
797 // LU decomposition
798 tilde.set_lu() ;
799
800 // Inversion using LAPACK
801 Tbl res0(tilde.inverse(aux)) ;
802 Tbl res(nr) ;
803 res.set_etat_qcq() ;
804
805 // ... finally, one has to recover the original starting base
806
807 res.set(0) = 3*res0(0) ;
808 for (int i=1; i<nr-2; i++) res.set(i) = res0(i-1)*(2*i-1)
809 + res0(i)*(2*i+3) ;
810 res.set(nr-2) = res0(nr-3)*(2*nr-5) ;
811 res.set(nr-1) = 0 ;
812
813 return res ;
814}
815
816Tbl _dal_inverse_r_chebi_o2_l(const Matrice &op, const Tbl &source,
817 const bool part) {
818
819 // Operator and source are copied and prepared
820 Matrice barre(op) ;
821 int nr = op.get_dim(0) ;
822 Tbl aux(nr-2) ;
823 if (part) {
824 aux.set_etat_qcq() ;
825 aux.set(nr-4) = source(nr-4) ;
826 aux.set(nr-3) = 0 ;
827 }
828 else {
829 aux.annule_hard() ;
830 aux.set(0) = 1. ;
831 }
832
833 // Operator is put into banded form (changing the image base ...)
834
835 for (int i=0; i<nr-4; i++) {
836 for (int j=i; j<nr; j++) {
837 barre.set(i,j) = (op(i+1,j) - op(i,j))/(i+1) ;
838 }
839 if (part)
840 aux.set(i) = (source(i+1) - source(i))/(i+1) ;
841 }
842 for (int i=0; i<nr-5; i++) {
843 for (int j=i; j<nr; j++) barre.set(i,j) = barre(i+2,j) - barre(i,j) ;
844 if (part) aux.set(i) = aux(i+2) - aux(i) ;
845 }
846
847 // ... and changing the starting base (first and last columns are quit) ...
848
849 Matrice tilde(nr-2,nr-2) ;
850 tilde.set_etat_qcq() ;
851 for (int i=0; i<nr-2; i++)
852 for (int j=0; j<nr-2;j++)
853 tilde.set(i,j) = barre(i,j+1)*(2*j+1) + barre(i,j)*(2*j+3) ;
854
855 // In this case the time-step is too large for the number of points
856 // (or the number of points too large for the time-step!)
857 // the matrix is ill-conditionned: the last line is put as the first
858 // one and all the others are shifted.
859
860 for (int i=nr-4; i>=0; i--) {
861 for (int j=((i>0) ? i-1 : 0); j<((i+5 > nr-2) ? nr-2 : i+5); j++)
862 tilde.set(i+1,j) = tilde(i,j) ;
863 if (part) aux.set(i+1) = aux(i) ;
864 }
865 tilde.set(0,0) = 1. ;
866 tilde.set(0,1) = 1. ;
867 tilde.set(0,2) = 1. ;
868 tilde.set(0,3) = 0. ;
869
870 if (part) aux.set(0) = 0 ;
871
872
873 // 2 over-diagonals and 2 under...
874 tilde.set_band(2,2) ;
875
876 // LU decomposition
877 tilde.set_lu() ;
878
879 // Inversion using LAPACK
880 Tbl res0(tilde.inverse(aux)) ;
881 Tbl res(nr) ;
882 res.set_etat_qcq() ;
883 // ... finally, one has to recover the original starting base
884 res.set(0) = 3*res0(0) ;
885 for (int i=1; i<nr-2; i++) res.set(i) = res0(i-1)*(2*i-1)
886 + res0(i)*(2*i+3) ;
887 res.set(nr-2) = res0(nr-3)*(2*nr-5) ;
888 res.set(nr-1) = 0 ;
889
890 return res ;
891}
892
893Tbl _dal_inverse_r_jaco02(const Matrice &op, const Tbl &source,
894 const bool part) {
895
896 // Operator and source are copied and prepared
897 Matrice barre(op) ;
898 int nr = op.get_dim(0) ;
899 Tbl aux(nr) ;
900 if (part) {
901 aux.set_etat_qcq() ;
902 aux.set(nr-2) = source(nr-2) ;
903 aux.set(nr-1) = 0 ;
904 }
905 else {
906 aux.annule_hard() ;
907 aux.set(0) = 1. ;
908 }
909
910 // Operator is put into banded form (changing the image base ...)
911
912 for (int i=0; i<nr; i++) {
913 for (int j=0; j<nr; j++) {
914 barre.set(i,j) = (op(i,j)) ;
915 }
916 if (part)
917 aux.set(i) = (source(i));
918 }
919 for (int i=0; i<nr; i++) {
920 for (int j=0; j<nr; j++) barre.set(i,j) = barre(i,j) ;
921 if (part) aux.set(i) = aux(i) ;
922 }
923
924 // ... and changing the starting base (first and last columns are quit) ...
925
926 Matrice tilde(nr,nr) ;
927 tilde.set_etat_qcq() ;
928 for (int i=0; i<nr; i++)
929 for (int j=0; j<nr;j++)
930 tilde.set(i,j) = barre(i,j) ;
931
932 // LU decomposition
933 tilde.set_lu() ;
934
935 // Inversion using LAPACK
936 Tbl res0(tilde.inverse(aux)) ;
937 Tbl res(nr) ;
938 res.set_etat_qcq() ;
939 // ... finally, one has to recover the original starting base
940 for (int i=0; i<nr; i++) res.set(i) = res0(i) ;
941
942 return res ;
943}
944
945
946
947 //----------------------------
948 //-- Fonction a appeler ---
949 //----------------------------
950
951
952Tbl dal_inverse(const int& base_r, const int& type_dal, const
953 Matrice& operateur, const Tbl& source, const bool part) {
954
955 // Routines de derivation
956 static Tbl (*dal_inverse[MAX_BASE][MAX_DAL])(const Matrice&, const Tbl&,
957 const bool) ;
958 static int nap = 0 ;
959
960 // Premier appel
961 if (nap==0) {
962 nap = 1 ;
963 for (int i=0 ; i<MAX_DAL ; i++) {
964 for (int j=0; j<MAX_BASE; j++)
965 dal_inverse[i][j] = _dal_inverse_pas_prevu ;
966 }
967 // Les routines existantes
968// dal_inverse[R_CHEB >> TRA_R] = _dal_inverse_r_cheb ; not good!!
969// dal_inverse[ORDRE1_SMALL][R_CHEB >> TRA_R] =
970// _dal_inverse_r_cheb_o1_s ;
971// dal_inverse[ORDRE1_LARGE][R_CHEB >> TRA_R] =
972// _dal_inverse_r_cheb_o1_l ;
973 dal_inverse[O2DEGE_SMALL][R_CHEB >> TRA_R] =
974 _dal_inverse_r_cheb_o2d_s ;
975 dal_inverse[O2DEGE_LARGE][R_CHEB >> TRA_R] =
976 _dal_inverse_r_cheb_o2d_l ;
977 dal_inverse[O2NOND_SMALL][R_CHEB >> TRA_R] =
978 _dal_inverse_r_cheb_o2_s ;
979 dal_inverse[O2NOND_LARGE][R_CHEB >> TRA_R] =
980 _dal_inverse_r_cheb_o2_l ;
981// dal_inverse[R_CHEB >> TRA_R] = _dal_inverse_r_cheb ; not good!!
982// dal_inverse[ORDRE1_SMALL][R_CHEBP >> TRA_R] =
983// _dal_inverse_r_chebp_o1_s ;
984// dal_inverse[ORDRE1_LARGE][R_CHEBP >> TRA_R] =
985// _dal_inverse_r_chebp_o1_l ;
986 dal_inverse[O2DEGE_SMALL][R_CHEBP >> TRA_R] =
987 _dal_inverse_r_chebp_o2d_s ;
988 dal_inverse[O2DEGE_LARGE][R_CHEBP >> TRA_R] =
989 _dal_inverse_r_chebp_o2d_l ;
990 dal_inverse[O2NOND_SMALL][R_CHEBP >> TRA_R] =
991 _dal_inverse_r_chebp_o2_s ;
992 dal_inverse[O2NOND_LARGE][R_CHEBP >> TRA_R] =
993 _dal_inverse_r_chebp_o2_l ;
994// dal_inverse[ORDRE1_SMALL][R_CHEBI >> TRA_R] =
995// _dal_inverse_r_chebi_o1_s ;
996// dal_inverse[ORDRE1_LARGE][R_CHEBI >> TRA_R] =
997// _dal_inverse_r_chebi_o1_l ;
998 dal_inverse[O2DEGE_SMALL][R_CHEBI >> TRA_R] =
999 _dal_inverse_r_chebi_o2d_s ;
1000 dal_inverse[O2DEGE_LARGE][R_CHEBI >> TRA_R] =
1001 _dal_inverse_r_chebi_o2d_l ;
1002 dal_inverse[O2NOND_SMALL][R_CHEBI >> TRA_R] =
1003 _dal_inverse_r_chebi_o2_s ;
1004 dal_inverse[O2NOND_LARGE][R_CHEBI >> TRA_R] =
1005 _dal_inverse_r_chebi_o2_l ;
1006// Only one routine pour Jacobi(0,2) polynomials
1007 dal_inverse[O2DEGE_SMALL][R_JACO02 >> TRA_R] =
1008 _dal_inverse_r_jaco02 ;
1009 dal_inverse[O2DEGE_LARGE][R_JACO02 >> TRA_R] =
1010 _dal_inverse_r_jaco02 ;
1011 dal_inverse[O2NOND_SMALL][R_JACO02 >> TRA_R] =
1012 _dal_inverse_r_jaco02 ;
1013 dal_inverse[O2NOND_LARGE][R_JACO02 >> TRA_R] =
1014 _dal_inverse_r_jaco02 ;
1015}
1016
1017 return dal_inverse[type_dal][base_r](operateur, source, part) ;
1018}
1019}
Lorene::Matrice
Matrix handling.
Definition matrice.h:152
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i]).
Definition tbl.h:403
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition type_parite.h:144
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition type_parite.h:188
MAX_DAL
#define MAX_DAL
Nombre max d'operateurs (pour l'instant).
Definition type_parite.h:274
O2DEGE_SMALL
#define O2DEGE_SMALL
Operateur du deuxieme ordre degenere .
Definition type_parite.h:280
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
O2NOND_LARGE
#define O2NOND_LARGE
Operateur du deuxieme ordre non degenere .
Definition type_parite.h:286
O2DEGE_LARGE
#define O2DEGE_LARGE
Operateur du deuxieme ordre degenere .
Definition type_parite.h:282
O2NOND_SMALL
#define O2NOND_SMALL
Operateur du deuxieme ordre non degenere .
Definition type_parite.h:284
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex).
Definition type_parite.h:158
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene
Lorene prototypes.
Definition app_hor.h:67