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