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  //------------------------------------
94 namespace Lorene {
95 Tbl _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 
108 Tbl _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 
147 Tbl _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 
217 Tbl _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 
256 Tbl _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 
332 Tbl _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 
391 Tbl _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 
467 Tbl _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 
529 Tbl _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 
613 Tbl _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 
676 Tbl _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 
755 Tbl _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 
816 Tbl _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 
893 Tbl _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 
952 Tbl 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 }
#define MAX_DAL
Nombre max d&#39;operateurs (pour l&#39;instant)
Definition: type_parite.h:274
#define O2DEGE_LARGE
Operateur du deuxieme ordre degenere .
Definition: type_parite.h:282
Lorene prototypes.
Definition: app_hor.h:67
#define O2NOND_LARGE
Operateur du deuxieme ordre non degenere .
Definition: type_parite.h:286
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define O2NOND_SMALL
Operateur du deuxieme ordre non degenere .
Definition: type_parite.h:284
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
#define O2DEGE_SMALL
Operateur du deuxieme ordre degenere .
Definition: type_parite.h:280
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166