cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
integrate.h
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 
4 #ifndef INTEGRATE_H_
5 #define INTEGRATE_H_
6 
7 #include "vectorize.h"
8 
9 // Define integration methods
11 
12 static const int numPoints = 32;
13 
14 ALIGNED(CD_ALIGN) static const double qg32_w[numPoints] = {
15  .35093050047350483e-2, .81371973654528350e-2, .12696032654631030e-1, .17136931456510717e-1,
16  .21417949011113340e-1, .25499029631188088e-1, .29342046739267774e-1, .32911111388180923e-1,
17  .36172897054424253e-1, .39096947893535153e-1, .41655962113473378e-1, .43826046502201906e-1,
18  .45586939347881942e-1, .46922199540402283e-1, .47819360039637430e-1, .48270044257363900e-1,
19  .48270044257363900e-1, .47819360039637430e-1, .46922199540402283e-1, .45586939347881942e-1,
20  .43826046502201906e-1, .41655962113473378e-1, .39096947893535153e-1, .36172897054424253e-1,
21  .32911111388180923e-1, .29342046739267774e-1, .25499029631188088e-1, .21417949011113340e-1,
22  .17136931456510717e-1, .12696032654631030e-1, .81371973654528350e-2, .35093050047350483e-2
23 };
24 
25 ALIGNED(CD_ALIGN) static const double qg32_a[numPoints] = {
26  -.498631930924740780, -.49280575577263417, -.4823811277937532200, -.46745303796886984000,
27  -.448160577883026060, -.42468380686628499, -.3972418979839712000, -.36609105937014484000,
28  -.331522133465107600, -.29385787862038116, -.2534499544661147000, -.21067563806531767000,
29  -.165934301141063820, -.11964368112606854, -.7223598079139825e-1, -.24153832843869158e-1,
30  .24153832843869158e-1, .7223598079139825e-1, .11964368112606854, .165934301141063820,
31  .21067563806531767000, .2534499544661147000, .29385787862038116, .331522133465107600,
32  .36609105937014484000, .3972418979839712000, .42468380686628499, .448160577883026060,
33  .46745303796886984000, .4823811277937532200, .49280575577263417, .498631930924740780
34 };
35 
36 // define an integrator class. Currently hard-wired to 32-point Gaussian
37 template<typename Integrand, methods Method>
39 {
40  const double *p_w, *p_a;
41 public:
42  Integrand func;
43  explicit Integrator( const Integrand& fun ) : p_w(qg32_w), p_a(qg32_a), func(fun) {}
44  double sum(double min, double max) const
45  {
46  if( Method == Gaussian32 )
47  {
48  double a = 0.5*(max+min);
49  double b = max-min;
50  double total = 0.;
51  for( long i=0; i < numPoints; i++ )
52  total += b * p_w[i] * func(a + b*p_a[i]);
53  return total;
54  }
55  else
56  TotalInsanity();
57  }
58 };
59 
60 // define a vectorized integrator class. Currently hard-wired to 32-point Gaussian
61 template<typename Integrand, methods Method>
63 {
64  const double *p_w, *p_a;
65 public:
66  Integrand func;
67  explicit VecIntegrator( const Integrand& fun ) : p_w(qg32_w), p_a(qg32_a), func(fun) {}
68  double sum(double min, double max) const
69  {
70  if( Method == VecGaussian32 )
71  {
72  double a = 0.5*(max+min);
73  double b = max-min;
74  ALIGNED(CD_ALIGN) double x[numPoints], y[numPoints];
75  for( long i=0; i < numPoints; i++ )
76  x[i] = a + b*p_a[i];
77  func(x, y, numPoints);
78  return b * reduce_ab(p_w, y, 0, numPoints);
79  }
80  else
81  TotalInsanity();
82  }
83 };
84 
91 double qg32( double, double, double(*)(double) );
92 /* declar of optimize_func, the last arg, changed from double(*)() to above,
93  * seemed to fix flags that were raised */
94 
95 // Open and closed integration routines, developed on qromb etc. in
96 // Press et al.
97 namespace integrate
98 {
99  template <class T>
100  class Trapezium
101  {
102  T m_f;
103  double m_a, m_b, m_sum;
104  long m_n;
105  public:
106  static const int NREF=2;
107  Trapezium(const T& f, double a, double b) : m_f(f), m_a(a), m_b(b), m_sum(0.0), m_n(0) {}
108  void step()
109  {
110  if (m_n == 0)
111  {
112  m_sum = 0.5*(m_b-m_a)*(m_f(m_a)+m_f(m_b));
113  m_n = 1;
114  }
115  else
116  {
117  double rn = 1.0/m_n;
118  double si = 0.0;
119  for (int i=0; i<m_n; ++i)
120  {
121  double x = (m_a*(m_n-i-0.5)+m_b*(i+0.5))*rn;
122  si += m_f(x);
123  }
124  m_sum = (m_sum+si*(m_b-m_a)*rn)/NREF;
125  m_n *= NREF;
126  }
127  }
128  double sum() const
129  {
130  return m_sum;
131  }
132  long evals() const
133  {
134  return m_n;
135  }
136  };
137 
138  template <class T>
139  class Midpoint
140  {
141  T m_f;
142  double m_a, m_b, m_sum;
143  long m_n;
144  public:
145  static const int NREF=3;
146  Midpoint(const T& f, double a, double b) : m_f(f), m_a(a), m_b(b), m_sum(0.0), m_n(0) {}
147  void step()
148  {
149  if (m_n == 0)
150  {
151  m_sum = (m_b-m_a)*m_f(0.5*(m_a+m_b));
152  m_n = 1;
153  }
154  else
155  {
156  double rn = 1.0/m_n;
157  const double sixth = 1./6.;
158  double si = 0.0;
159  for (int i=0; i<m_n; ++i)
160  {
161  double x1 = (m_a*(m_n-i-sixth)+m_b*(i+sixth))*rn;
162  si += m_f(x1);
163  double x2 = (m_a*(m_n-i-1+sixth)+m_b*(i+1-sixth))*rn;
164  si += m_f(x2);
165  }
166  m_sum = (m_sum+si*(m_b-m_a)*rn)/NREF;
167  m_n *= NREF;
168  }
169  }
170  double sum() const
171  {
172  return m_sum;
173  }
174  long evals() const
175  {
176  return m_n;
177  }
178  };
179 
180  template <class T>
181  class Romberg
182  {
183  T m_f;
184  double m_sum, m_dy;
185  public:
186  explicit Romberg(const T& f) : m_f(f), m_sum(0.0), m_dy(0.0) {}
187  void update(double eps);
188  double sum() const
189  {
190  return m_sum;
191  }
192  double error() const
193  {
194  return m_dy;
195  }
196  long evals() const
197  {
198  return m_f.evals();
199  }
200  };
201 
202  template <class T>
203  inline void Romberg<T>::update(double eps)
204  {
205  const int itmax=40/m_f.NREF, npt=5;
206  double d1[npt], d0[npt-1];
207  double y = 0.0;
208  const double w1 = m_f.NREF*m_f.NREF, w2 = 1.0/w1;
209  for (int i=0; i<npt; ++i)
210  {
211  d1[i] = 0.;
212  }
213  for (int i=0; i<itmax; ++i)
214  {
215  m_f.step();
216  y = m_f.sum();
217 
218  int l = (i<npt-1) ? i : npt-1;
219  for (int m=0; m<l; ++m)
220  {
221  d0[m] = d1[m];
222  }
223  d1[0] = y;
224  double fr = 1.0;
225  for (int m=0; m<l; ++m)
226  {
227  d1[m+1] = (d1[m]-d0[m]*fr)/(w1-fr);
228  y += d1[m+1];
229  fr *= w2;
230  }
231 
232  m_dy = fabs(m_sum-y);
233  m_sum = y;
234  //printf("Level %d value %g->%g error %g samples %ld\n",l,d1[0],m_sum,m_dy,m_f.evals());
235  if ( i > 2 && m_dy <= eps*fabs(y) )
236  {
237  return;
238  }
239  }
240  return;
241  }
242 
243  template <class T>
244  class Simple
245  {
246  T m_f;
247  double m_sum, m_dy;
248  public:
249  explicit Simple(const T& f) : m_f(f), m_sum(0.0), m_dy(0.0) {}
250  void update(double eps);
251  double sum() const
252  {
253  return m_sum;
254  }
255  double error() const
256  {
257  return m_dy;
258  }
259  long evals() const
260  {
261  return m_f.evals();
262  }
263  };
264 
265  template <class T>
266  inline void Simple<T>::update(double eps)
267  {
268  const int itmax=40/m_f.NREF;
269  double coarse=0.,fine=0.;
270  for (int i=0; i<itmax; ++i)
271  {
272  coarse = fine;
273  m_f.step();
274  fine = m_f.sum();
275 
276  m_dy = fabs(coarse-fine);
277  m_sum = fine;
278  // printf("Level %d value %g error %g samples %ld\n",i,m_sum,m_dy,m_f.evals());
279  if ( i > 2 && m_dy <= eps*fabs(fine) )
280  {
281  return;
282  }
283  }
284  return;
285  }
286 }
287 #endif /* INTEGRATE_H_ */
Integrand func
Definition: integrate.h:66
long evals() const
Definition: integrate.h:196
Midpoint(const T &f, double a, double b)
Definition: integrate.h:146
static double x2[63]
static const int numPoints
Definition: integrate.h:12
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
static double x1[83]
long evals() const
Definition: integrate.h:174
const double * p_w
Definition: integrate.h:40
void update(double eps)
Definition: integrate.h:203
ALIGNED(CD_ALIGN) static const double qg32_w[numPoints]
Romberg(const T &f)
Definition: integrate.h:186
double sum() const
Definition: integrate.h:128
long evals() const
Definition: integrate.h:259
double sum() const
Definition: integrate.h:188
methods
Definition: integrate.h:10
Trapezium(const T &f, double a, double b)
Definition: integrate.h:107
const double * p_a
Definition: integrate.h:40
double sum(double min, double max) const
Definition: integrate.h:68
static const int NREF
Definition: integrate.h:145
double sum() const
Definition: integrate.h:170
long max(int a, long b)
Definition: cddefines.h:817
const double * p_a
Definition: integrate.h:64
long min(int a, long b)
Definition: cddefines.h:762
double sum(double min, double max) const
Definition: integrate.h:44
double sum() const
Definition: integrate.h:251
long evals() const
Definition: integrate.h:132
Integrand func
Definition: integrate.h:42
void update(double eps)
Definition: integrate.h:266
const double * p_w
Definition: integrate.h:64
double reduce_ab(const double *a, const double *b, long ilo, long ihi)
double qg32(double, double, double(*)(double))
Definition: service.cpp:1175
static const int NREF
Definition: integrate.h:106
double error() const
Definition: integrate.h:255
VecIntegrator(const Integrand &fun)
Definition: integrate.h:67
Simple(const T &f)
Definition: integrate.h:249
#define CD_ALIGN
Definition: cpu.h:156
double error() const
Definition: integrate.h:192
Integrator(const Integrand &fun)
Definition: integrate.h:43