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 char 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  //------------------------------------
91 namespace Lorene {
92 Tbl _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 
105 Tbl _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 
144 Tbl _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 
214 Tbl _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 
253 Tbl _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 
329 Tbl _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 
388 Tbl _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 
464 Tbl _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 
526 Tbl _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 
610 Tbl _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 
673 Tbl _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 
752 Tbl _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 
813 Tbl _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 
890 Tbl _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 
949 Tbl 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.
Definition: type_parite.h:144
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define MAX_DAL
Nombre max d'operateurs (pour l'instant)
Definition: type_parite.h:274
#define O2DEGE_SMALL
Operateur du deuxieme ordre degenere .
Definition: type_parite.h:280
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define O2NOND_LARGE
Operateur du deuxieme ordre non degenere .
Definition: type_parite.h:286
#define O2DEGE_LARGE
Operateur du deuxieme ordre degenere .
Definition: type_parite.h:282
#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_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Lorene prototypes.
Definition: app_hor.h:64