inte_kronrod_gsl.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2015, Jerry Gagelman
5  and Andrew W. Steiner
6 
7  This file is part of O2scl.
8 
9  O2scl is free software; you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation; either version 3 of the License, or
12  (at your option) any later version.
13 
14  O2scl is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
21 
22  -------------------------------------------------------------------
23 */
24 #ifndef O2SCL_INTE_GSL_KRONROD_H
25 #define O2SCL_INTE_GSL_KRONROD_H
26 
27 /** \file inte_kronrod_gsl.h
28  \brief File defining GSL-based integration coefficients, workspace,
29  and \ref o2scl::inte_kronrod_gsl
30 */
31 
32 #include <cmath>
33 #include <gsl/gsl_integration.h>
34 #include <o2scl/inte.h>
35 #include <o2scl/inte_gsl.h>
36 #include <o2scl/string_conv.h>
37 
38 /** \brief A namespace for the GSL adaptive Gauss-Kronrod integration
39  coefficients
40 
41  The storage scheme follows that of QUADPACK: symmetry about the
42  origin permits storing only those abscissae in the interval
43  \f$ [0, 1). \f$ Similiarly for the weights. The odd-indexed
44  abscissae \c xgk[1], \c xgk[3],... belong to the low-order Gauss
45  rule. For the full Gauss-Kronrod rule, the array \c wgk[] holds
46  the weights corresponding to the abscissae \c xgk[]. For the
47  low-order rule, \c wg[i] is the weight correpsonding to
48  \c xgk[2*i+1].
49 
50  <b>Documentation from GSL</b>: \n
51  Gauss quadrature weights and kronrod quadrature abscissae and
52  weights as evaluated with 80 decimal digit arithmetic by
53  L. W. Fullerton, Bell Labs, Nov. 1981.
54 */
56 
57  /** \brief Abscissae of the 15-point Kronrod rule */
58  static const double qk15_xgk[8]={
59  0.991455371120812639206854697526329,
60  0.949107912342758524526189684047851,
61  0.864864423359769072789712788640926,
62  0.741531185599394439863864773280788,
63  0.586087235467691130294144838258730,
64  0.405845151377397166906606412076961,
65  0.207784955007898467600689403773245,
66  0.000000000000000000000000000000000
67  };
68 
69  /** \brief Weights of the 7-point Gauss rule */
70  static const double qk15_wg[4] =
71  {
72  0.129484966168869693270611432679082,
73  0.279705391489276667901467771423780,
74  0.381830050505118944950369775488975,
75  0.417959183673469387755102040816327
76  };
77 
78  /** \brief Weights of the 15-point Kronrod rule */
79  static const double qk15_wgk[8] =
80  {
81  0.022935322010529224963732008058970,
82  0.063092092629978553290700663189204,
83  0.104790010322250183839876322541518,
84  0.140653259715525918745189590510238,
85  0.169004726639267902826583426598550,
86  0.190350578064785409913256402421014,
87  0.204432940075298892414161999234649,
88  0.209482141084727828012999174891714
89  };
90 
91  /** \brief Abscissae of the 21-point Kronrod rule */
92  static const double qk21_xgk[11] =
93  {
94  0.995657163025808080735527280689003,
95  0.973906528517171720077964012084452,
96  0.930157491355708226001207180059508,
97  0.865063366688984510732096688423493,
98  0.780817726586416897063717578345042,
99  0.679409568299024406234327365114874,
100  0.562757134668604683339000099272694,
101  0.433395394129247190799265943165784,
102  0.294392862701460198131126603103866,
103  0.148874338981631210884826001129720,
104  0.000000000000000000000000000000000
105  };
106 
107  /** \brief Weights of the 10-point Gauss rule */
108  static const double qk21_wg[5] =
109  {
110  0.066671344308688137593568809893332,
111  0.149451349150580593145776339657697,
112  0.219086362515982043995534934228163,
113  0.269266719309996355091226921569469,
114  0.295524224714752870173892994651338
115  };
116 
117  /** \brief Weights of the 21-point Kronrod rule */
118  static const double qk21_wgk[11] =
119  {
120  0.011694638867371874278064396062192,
121  0.032558162307964727478818972459390,
122  0.054755896574351996031381300244580,
123  0.075039674810919952767043140916190,
124  0.093125454583697605535065465083366,
125  0.109387158802297641899210590325805,
126  0.123491976262065851077958109831074,
127  0.134709217311473325928054001771707,
128  0.142775938577060080797094273138717,
129  0.147739104901338491374841515972068,
130  0.149445554002916905664936468389821
131  };
132 
133  /** \brief Abscissae of the 31-point Kronrod rule */
134  static const double qk31_xgk[16] =
135  {
136  0.998002298693397060285172840152271,
137  0.987992518020485428489565718586613,
138  0.967739075679139134257347978784337,
139  0.937273392400705904307758947710209,
140  0.897264532344081900882509656454496,
141  0.848206583410427216200648320774217,
142  0.790418501442465932967649294817947,
143  0.724417731360170047416186054613938,
144  0.650996741297416970533735895313275,
145  0.570972172608538847537226737253911,
146  0.485081863640239680693655740232351,
147  0.394151347077563369897207370981045,
148  0.299180007153168812166780024266389,
149  0.201194093997434522300628303394596,
150  0.101142066918717499027074231447392,
151  0.000000000000000000000000000000000
152  };
153 
154  /** \brief Weights of the 15-point Gauss rule */
155  static const double qk31_wg[8] =
156  {
157  0.030753241996117268354628393577204,
158  0.070366047488108124709267416450667,
159  0.107159220467171935011869546685869,
160  0.139570677926154314447804794511028,
161  0.166269205816993933553200860481209,
162  0.186161000015562211026800561866423,
163  0.198431485327111576456118326443839,
164  0.202578241925561272880620199967519
165  };
166 
167  /** \brief Weights of the 31-point Kronrod rule */
168  static const double qk31_wgk[16] =
169  {
170  0.005377479872923348987792051430128,
171  0.015007947329316122538374763075807,
172  0.025460847326715320186874001019653,
173  0.035346360791375846222037948478360,
174  0.044589751324764876608227299373280,
175  0.053481524690928087265343147239430,
176  0.062009567800670640285139230960803,
177  0.069854121318728258709520077099147,
178  0.076849680757720378894432777482659,
179  0.083080502823133021038289247286104,
180  0.088564443056211770647275443693774,
181  0.093126598170825321225486872747346,
182  0.096642726983623678505179907627589,
183  0.099173598721791959332393173484603,
184  0.100769845523875595044946662617570,
185  0.101330007014791549017374792767493
186  };
187 
188  /** \brief Abscissae of the 41-point Kronrod rule */
189  static const double qk41_xgk[21] =
190  {
191  0.998859031588277663838315576545863,
192  0.993128599185094924786122388471320,
193  0.981507877450250259193342994720217,
194  0.963971927277913791267666131197277,
195  0.940822633831754753519982722212443,
196  0.912234428251325905867752441203298,
197  0.878276811252281976077442995113078,
198  0.839116971822218823394529061701521,
199  0.795041428837551198350638833272788,
200  0.746331906460150792614305070355642,
201  0.693237656334751384805490711845932,
202  0.636053680726515025452836696226286,
203  0.575140446819710315342946036586425,
204  0.510867001950827098004364050955251,
205  0.443593175238725103199992213492640,
206  0.373706088715419560672548177024927,
207  0.301627868114913004320555356858592,
208  0.227785851141645078080496195368575,
209  0.152605465240922675505220241022678,
210  0.076526521133497333754640409398838,
211  0.000000000000000000000000000000000
212  };
213 
214  /** \brief Weights of the 20-point Gauss rule */
215  static const double qk41_wg[11] =
216  {
217  0.017614007139152118311861962351853,
218  0.040601429800386941331039952274932,
219  0.062672048334109063569506535187042,
220  0.083276741576704748724758143222046,
221  0.101930119817240435036750135480350,
222  0.118194531961518417312377377711382,
223  0.131688638449176626898494499748163,
224  0.142096109318382051329298325067165,
225  0.149172986472603746787828737001969,
226  0.152753387130725850698084331955098
227  };
228 
229  /** \brief Weights of the 41-point Kronrod rule */
230  static const double qk41_wgk[21] =
231  {
232  0.003073583718520531501218293246031,
233  0.008600269855642942198661787950102,
234  0.014626169256971252983787960308868,
235  0.020388373461266523598010231432755,
236  0.025882133604951158834505067096153,
237  0.031287306777032798958543119323801,
238  0.036600169758200798030557240707211,
239  0.041668873327973686263788305936895,
240  0.046434821867497674720231880926108,
241  0.050944573923728691932707670050345,
242  0.055195105348285994744832372419777,
243  0.059111400880639572374967220648594,
244  0.062653237554781168025870122174255,
245  0.065834597133618422111563556969398,
246  0.068648672928521619345623411885368,
247  0.071054423553444068305790361723210,
248  0.073030690332786667495189417658913,
249  0.074582875400499188986581418362488,
250  0.075704497684556674659542775376617,
251  0.076377867672080736705502835038061,
252  0.076600711917999656445049901530102
253  };
254 
255  /** \brief Abscissae of the 51-point Kronrod rule */
256  static const double qk51_xgk[26] =
257  {
258  0.999262104992609834193457486540341,
259  0.995556969790498097908784946893902,
260  0.988035794534077247637331014577406,
261  0.976663921459517511498315386479594,
262  0.961614986425842512418130033660167,
263  0.942974571228974339414011169658471,
264  0.920747115281701561746346084546331,
265  0.894991997878275368851042006782805,
266  0.865847065293275595448996969588340,
267  0.833442628760834001421021108693570,
268  0.797873797998500059410410904994307,
269  0.759259263037357630577282865204361,
270  0.717766406813084388186654079773298,
271  0.673566368473468364485120633247622,
272  0.626810099010317412788122681624518,
273  0.577662930241222967723689841612654,
274  0.526325284334719182599623778158010,
275  0.473002731445714960522182115009192,
276  0.417885382193037748851814394594572,
277  0.361172305809387837735821730127641,
278  0.303089538931107830167478909980339,
279  0.243866883720988432045190362797452,
280  0.183718939421048892015969888759528,
281  0.122864692610710396387359818808037,
282  0.061544483005685078886546392366797,
283  0.000000000000000000000000000000000
284  };
285 
286  /** \brief Weights of the 25-point Gauss rule */
287  static const double qk51_wg[13] =
288  {
289  0.011393798501026287947902964113235,
290  0.026354986615032137261901815295299,
291  0.040939156701306312655623487711646,
292  0.054904695975835191925936891540473,
293  0.068038333812356917207187185656708,
294  0.080140700335001018013234959669111,
295  0.091028261982963649811497220702892,
296  0.100535949067050644202206890392686,
297  0.108519624474263653116093957050117,
298  0.114858259145711648339325545869556,
299  0.119455763535784772228178126512901,
300  0.122242442990310041688959518945852,
301  0.123176053726715451203902873079050
302  };
303 
304  /** \brief Weights of the 51-point Kronrod rule */
305  static const double qk51_wgk[26] =
306  {
307  0.001987383892330315926507851882843,
308  0.005561932135356713758040236901066,
309  0.009473973386174151607207710523655,
310  0.013236229195571674813656405846976,
311  0.016847817709128298231516667536336,
312  0.020435371145882835456568292235939,
313  0.024009945606953216220092489164881,
314  0.027475317587851737802948455517811,
315  0.030792300167387488891109020215229,
316  0.034002130274329337836748795229551,
317  0.037116271483415543560330625367620,
318  0.040083825504032382074839284467076,
319  0.042872845020170049476895792439495,
320  0.045502913049921788909870584752660,
321  0.047982537138836713906392255756915,
322  0.050277679080715671963325259433440,
323  0.052362885806407475864366712137873,
324  0.054251129888545490144543370459876,
325  0.055950811220412317308240686382747,
326  0.057437116361567832853582693939506,
327  0.058689680022394207961974175856788,
328  0.059720340324174059979099291932562,
329  0.060539455376045862945360267517565,
330  0.061128509717053048305859030416293,
331  0.061471189871425316661544131965264,
332  0.061580818067832935078759824240066
333  };
334 
335  /** \brief Abscissae of the 61-point Kronrod rule */
336  static const double qk61_xgk[31] =
337  {
338  0.999484410050490637571325895705811,
339  0.996893484074649540271630050918695,
340  0.991630996870404594858628366109486,
341  0.983668123279747209970032581605663,
342  0.973116322501126268374693868423707,
343  0.960021864968307512216871025581798,
344  0.944374444748559979415831324037439,
345  0.926200047429274325879324277080474,
346  0.905573307699907798546522558925958,
347  0.882560535792052681543116462530226,
348  0.857205233546061098958658510658944,
349  0.829565762382768397442898119732502,
350  0.799727835821839083013668942322683,
351  0.767777432104826194917977340974503,
352  0.733790062453226804726171131369528,
353  0.697850494793315796932292388026640,
354  0.660061064126626961370053668149271,
355  0.620526182989242861140477556431189,
356  0.579345235826361691756024932172540,
357  0.536624148142019899264169793311073,
358  0.492480467861778574993693061207709,
359  0.447033769538089176780609900322854,
360  0.400401254830394392535476211542661,
361  0.352704725530878113471037207089374,
362  0.304073202273625077372677107199257,
363  0.254636926167889846439805129817805,
364  0.204525116682309891438957671002025,
365  0.153869913608583546963794672743256,
366  0.102806937966737030147096751318001,
367  0.051471842555317695833025213166723,
368  0.000000000000000000000000000000000
369  };
370 
371  /** \brief Weights of the 30-point Gauss rule */
372  static const double qk61_wg[15] =
373  {
374  0.007968192496166605615465883474674,
375  0.018466468311090959142302131912047,
376  0.028784707883323369349719179611292,
377  0.038799192569627049596801936446348,
378  0.048402672830594052902938140422808,
379  0.057493156217619066481721689402056,
380  0.065974229882180495128128515115962,
381  0.073755974737705206268243850022191,
382  0.080755895229420215354694938460530,
383  0.086899787201082979802387530715126,
384  0.092122522237786128717632707087619,
385  0.096368737174644259639468626351810,
386  0.099593420586795267062780282103569,
387  0.101762389748405504596428952168554,
388  0.102852652893558840341285636705415
389  };
390 
391  /** \brief Weights of the 61-point Kronrod rule */
392  static const double qk61_wgk[31] =
393  {
394  0.001389013698677007624551591226760,
395  0.003890461127099884051267201844516,
396  0.006630703915931292173319826369750,
397  0.009273279659517763428441146892024,
398  0.011823015253496341742232898853251,
399  0.014369729507045804812451432443580,
400  0.016920889189053272627572289420322,
401  0.019414141193942381173408951050128,
402  0.021828035821609192297167485738339,
403  0.024191162078080601365686370725232,
404  0.026509954882333101610601709335075,
405  0.028754048765041292843978785354334,
406  0.030907257562387762472884252943092,
407  0.032981447057483726031814191016854,
408  0.034979338028060024137499670731468,
409  0.036882364651821229223911065617136,
410  0.038678945624727592950348651532281,
411  0.040374538951535959111995279752468,
412  0.041969810215164246147147541285970,
413  0.043452539701356069316831728117073,
414  0.044814800133162663192355551616723,
415  0.046059238271006988116271735559374,
416  0.047185546569299153945261478181099,
417  0.048185861757087129140779492298305,
418  0.049055434555029778887528165367238,
419  0.049795683427074206357811569379942,
420  0.050405921402782346840893085653585,
421  0.050881795898749606492297473049805,
422  0.051221547849258772170656282604944,
423  0.051426128537459025933862879215781,
424  0.051494729429451567558340433647099
425  };
426 
427 }
428 
429 #ifndef DOXYGEN_NO_O2NS
430 namespace o2scl {
431 #endif
432 
433  /** \brief Integration workspace for the GSL integrators
434 
435  This is a simple rewrite of \c inte_gslgration_workspace
436  into a class.
437 
438  QUADPACK workspace documentation:
439  \verbatim
440  c parameters (meaning at output)
441  c limit - integer
442  c maximum number of error estimates the list
443  c can contain
444  c last - integer
445  c number of error estimates currently in the list
446  c maxerr - integer
447  c maxerr points to the nrmax-th largest error
448  c estimate currently in the list
449  c ermax - double precision
450  c nrmax-th largest error estimate
451  c ermax = elist(maxerr)
452  c elist - double precision
453  c vector of dimension last containing
454  c the error estimates
455  c iord - integer
456  c vector of dimension last, the first k elements
457  c of which contain pointers to the error
458  c estimates, such that
459  c elist(iord(1)),..., elist(iord(k))
460  c form a decreasing sequence, with
461  c k = last if last.le.(limit/2+2), and
462  c k = limit+1-last otherwise
463  c nrmax - integer
464  c maxerr = iord(nrmax)
465  \endverbatim
466  \verbatim
467  c alist - real
468  c vector of dimension at least limit, the first
469  c last elements of which are the left
470  c end points of the subintervals in the partition
471  c of the given integration range (a,b)
472  c blist - real
473  c vector of dimension at least limit, the first
474  c last elements of which are the right
475  c end points of the subintervals in the partition
476  c of the given integration range (a,b)
477  c rlist - real
478  c vector of dimension at least limit, the first
479  c last elements of which are the
480  c integral approximations on the subintervals
481  c elist - real
482  c vector of dimension at least limit, the first
483  c last elements of which are the moduli of the
484  c absolute error estimates on the subintervals
485  c iord - integer
486  c vector of dimension at least limit, the first k
487  c elements of which are pointers to the
488  c error estimates over the subintervals,
489  c such that elist(iord(1)), ...,
490  c elist(iord(k)) form a decreasing sequence,
491  c with k = last if last.le.(limit/2+2), and
492  c k = limit+1-last otherwise
493  c last - integer
494  c number of subintervals actually produced in the
495  c subdivision process
496  \endverbatim
497 
498  */
500 
501  public:
502 
504 
506 
507  /// Maximum number of subintervals allocated
508  size_t limit;
509  /// Current number of subintervals being used
510  size_t size;
511  /// Counter for extrapolation routine
512  size_t nrmax;
513  /// Index of current subinterval
514  size_t i;
515  /// Depth of subdivisions reached
517  /// Left endpoints of subintervals
518  double *alist;
519  /// Right endpoints of subintervals
520  double *blist;
521  /// Integral approximations on subintervals
522  double *rlist;
523  /// Integral error estimates
524  double *elist;
525  /// Linear ordering vector for sort routine
526  size_t *order;
527  /// Numbers of subdivisions made
528  size_t *level;
529 
530  /// Allocate a workspace
531  int allocate(size_t sz);
532 
533  /// Free allocated workspace memory
534  int free();
535 
536  /** \brief Initialize the workspace for an integration with limits
537  \c a and \c b.
538  */
539  int initialise(double a, double b);
540 
541  /** \brief Update the workspace with the result and error
542  from the first integration
543  */
544  int set_initial_result(double result, double error);
545 
546  /** \brief Retrieve the ith result from the
547  workspace stack
548 
549  The workspace variable \c i is used to specify which
550  interval is requested.
551  */
552  int retrieve(double *a, double *b, double *r, double *e) const;
553 
554  /** \brief Sort the workspace stack
555 
556  This routine maintains the descending ordering in the list of
557  the local error estimated resulting from the interval
558  subdivision process. at each call two error estimates are
559  inserted using the sequential search method, top-down for the
560  largest error estimate and bottom-up for the smallest error
561  estimate.
562 
563  Originally written in QUADPACK by R. Piessens and E. de
564  Doncker, translated into C for GSL by Brian Gough, and then
565  rewritten for \o2.
566  */
567  int qpsrt();
568 
569  /** \brief Determine which new subinterval to add to the workspace
570  stack and perform update
571  */
572  int update(double a1, double b1, double area1, double error1,
573  double a2, double b2, double area2, double error2);
574 
575  /// Add up all of the contributions to construct the final result
576  double sum_results();
577 
578  /** \brief Test whether the proposed subdivision falls
579  before floating-point precision
580  */
581  int subinterval_too_small(double a1, double a2, double b2);
582 
583  /// Push a new interval to the workspace stack
584  void append_interval(double a1, double b1, double area1, double error1);
585 
586  };
587 
588  /** \brief Basic Gauss-Kronrod integration class (GSL)
589 
590  This class provides the basic Gauss-Kronrod integration function
591  and integration workspace for some of the GSL-based integration
592  classes.
593 
594  The main function of interest is \ref set_rule(), which
595  sets the integration rule for the GSL integration classes
596  which inherit from this base class. The argument to
597  set rule should be selected from the following list:
598  - 1: 15-point Gauss-Kronrod integration rule (default)
599  - 2: 21-point rule
600  - 3: 31-point rule
601  - 4: 41-point rule
602  - 5: 51-point rule
603  - 6: 61-point rule
604  Any value other than 1-6 forces the error handler to be
605  called. All classes default to the 15-point rule,
606  except for \ref inte_qags_gsl which defaults to
607  the 21-point rule for singularities.
608 
609  The integration coefficients for use with this class and
610  associated children are stored in the \ref o2scl_inte_gk_coeffs
611  namespace.
612 
613  \future Convert 'double *' objects to ubvector
614  */
615  template<class func_t> class inte_kronrod_gsl : public inte_gsl,
616  public inte<func_t> {
617 
618  protected:
619 
620  /// The integration workspace
622 
623  /// Size of Gauss-Kronrod arrays
624  int n_gk;
625 
626  /// Gauss-Kronrod abscissae pointer
627  const double *x_gk;
628 
629  /// Gauss weight pointer
630  const double *w_g;
631 
632  /// Gauss-Kronrod weight pointer
633  const double *w_gk;
634 
635  /// Scratch space
636  double *f_v1;
637 
638  /// Scratch space
639  double *f_v2;
640 
641  public:
642 
643  inte_kronrod_gsl() {
644  w=new inte_workspace_gsl;
645  w->allocate(1000);
646  n_gk=0;
647  set_rule(1);
648  }
649 
650  ~inte_kronrod_gsl() {
651  w->free();
652  delete w;
653  if (n_gk>0) {
654  delete[] f_v1;
655  delete[] f_v2;
656  }
657  }
658 
659  /** \brief Get the Gauss-Kronrod integration rule
660 
661  This returns the index of the GSL integration rule
662  a number between 1 and 6 (inclusive)
663  */
664  int get_rule() {
665  switch (n_gk) {
666  case 8:
667  return 1;
668  case 11:
669  return 2;
670  case 16:
671  return 3;
672  case 21:
673  return 4;
674  case 26:
675  return 5;
676  case 31:
677  return 6;
678  }
679  return 0;
680  }
681 
682  /// Set the Gauss-Kronrod integration rule to be used
683  void set_rule(int rule) {
684 
685  using namespace o2scl_inte_gk_coeffs;
686 
687  if (n_gk > 0) {
688  delete[] f_v1;
689  delete[] f_v2;
690  }
691 
692  switch (rule) {
693  case 1:
694  n_gk=8;
695  x_gk=qk15_xgk;
696  w_g=qk15_wg;
697  w_gk=qk15_wgk;
698  break;
699  case 2:
700  n_gk=11;
701  x_gk=qk21_xgk;
702  w_g=qk21_wg;
703  w_gk=qk21_wgk;
704  break;
705  case 3:
706  n_gk=16;
707  x_gk=qk31_xgk;
708  w_g=qk31_wg;
709  w_gk=qk31_wgk;
710  break;
711  case 4:
712  n_gk=21;
713  x_gk=qk41_xgk;
714  w_g=qk41_wg;
715  w_gk=qk41_wgk;
716  break;
717  case 5:
718  n_gk=26;
719  x_gk=qk51_xgk;
720  w_g=qk51_wg;
721  w_gk=qk51_wgk;
722  break;
723  case 6:
724  n_gk=31;
725  x_gk=qk61_xgk;
726  w_g=qk61_wg;
727  w_gk=qk61_wgk;
728  break;
729  default:
730  O2SCL_ERR("Invalid rule.",exc_einval);
731  break;
732  }
733 
734  f_v1=new double[n_gk];
735  f_v2=new double[n_gk];
736 
737  return;
738  }
739 
740  /** \brief Set the limit for the number of subdivisions of
741  the integration region (default 1000)
742 
743  If the value of \c size is zero, the error handler will
744  be called.
745  */
746  int set_limit(size_t lim) {
747  if (lim==0) {
748  O2SCL_ERR2("Requested zero limit in ",
749  "inte_kronrod_gsl::set_limit().",exc_einval);
750  }
751  w->free();
752  w->allocate(lim);
753  return 0;
754  }
755 
756  /** \brief The base Gauss-Kronrod integration function template
757 
758  Given abcissas and weights, this performs the integration of
759  \c func between \c a and \c b, providing a result with
760  uncertainties.
761 
762  The Gauss-Kronrod rule uses \f$ 2m+1 \f$ abscissae \f$ x_1,
763  \ldots, x_{2m+1} \in (a, b) \f$ to estimate the integral of a
764  function as a linear combination of values,
765  \f[
766  \int_a^b f~dx \; \approx \;
767  Q_m^{GK}f \; \equiv \;
768  \sum_{k=1}^{2m+1} \beta_k f(x_k),
769  \f]
770  where the weights \f$ \beta_1,\ldots, \beta_{2m+1} \f$ are
771  intrinsic to the abscissae. The data are designed so that the
772  even-indexed abscissae yield a coarser estimate,
773  \f[
774  Q_m^{G}f \; \equiv \;
775  \sum_{j=1}^{m} \alpha_j f(x_{2j}),
776  \f]
777  and their difference \f$ |Q_m^Gf - Q_m^{GK}f| \f$ is the "raw"
778  error estimate. The various quantities that the function
779  computes are
780  \f{eqnarray*}{
781  \mathtt{result} &=& Q_m^{GK}f, \\
782  \mathtt{resabs} &=& Q_m^{GK}|f|, \\
783  \mathtt{resasc} &=& Q_m^{GK}(|f| - \mu),
784  \quad \mu \equiv \frac{Q_m^{GK}f}{b - a}.
785  \f}
786  The "absolute" error \c abserr is computed from the raw error
787  value using the function \ref inte_gsl::rescale_error.
788 
789  This function is designed for use with the values given in the
790  o2scl_inte_gk_coeffs namespace.
791 
792  This function never calls the error handler.
793 
794  \future This function, in principle, is an integration
795  object in and of itself, and could be implemented separately
796  as an object of type \ref o2scl::inte.
797  */
798  template<class func2_t> void gauss_kronrod_base
799  (func2_t &func, double a, double b, double *result,
800  double *abserr, double *resabs, double *resasc) {
801 
802  const double center=0.5*(a+b);
803  const double half_length=0.5*(b-a);
804  const double abs_half_length=fabs (half_length);
805 
806  double f_center;
807  f_center=func(center);
808 
809  double result_gauss=0.0;
810  double result_kronrod=f_center*this->w_gk[this->n_gk-1];
811 
812  double result_abs=fabs (result_kronrod);
813  double result_asc=0.0;
814  double mean=0.0, err=0.0;
815 
816  if (this->n_gk % 2 == 0) {
817  result_gauss=f_center*this->w_g[this->n_gk / 2-1];
818  }
819 
820  // Sum up results from left half of interval
821  for (int j=0; j < (this->n_gk-1) / 2; j++) {
822 
823  const int jtw=j*2+1;
824  const double abscissa=half_length*this->x_gk[jtw];
825 
826  double fval1, fval2;
827  fval1=func(center-abscissa);
828  fval2=func(center+abscissa);
829 
830  const double fsum=fval1+fval2;
831  this->f_v1[jtw]=fval1;
832  this->f_v2[jtw]=fval2;
833  result_gauss+=this->w_g[j]*fsum;
834  result_kronrod+=this->w_gk[jtw]*fsum;
835  result_abs+=this->w_gk[jtw]*(fabs(fval1)+fabs(fval2));
836  }
837 
838  // Sum up results from right half of interval
839  for (int j=0; j < this->n_gk / 2; j++) {
840 
841  int jtwm1=j*2;
842  const double abscissa=half_length*this->x_gk[jtwm1];
843 
844  double fval1, fval2;
845  fval1=func(center-abscissa);
846  fval2=func(center+abscissa);
847 
848  this->f_v1[jtwm1]=fval1;
849  this->f_v2[jtwm1]=fval2;
850  result_kronrod+=this->w_gk[jtwm1]*(fval1+fval2);
851  result_abs+=this->w_gk[jtwm1]*(fabs(fval1)+fabs(fval2));
852  }
853 
854  mean=result_kronrod*0.5;
855 
856  result_asc=this->w_gk[this->n_gk-1]*fabs(f_center-mean);
857 
858  for (int j=0;j<this->n_gk-1;j++) {
859  result_asc+=this->w_gk[j]*(fabs(this->f_v1[j]-mean)+
860  fabs(this->f_v2[j]-mean));
861  }
862 
863  /* Scale by the width of the integration region */
864 
865  err=(result_kronrod-result_gauss)*half_length;
866 
867  result_kronrod*=half_length;
868  result_abs*=abs_half_length;
869  result_asc*=abs_half_length;
870 
871  *result=result_kronrod;
872  *resabs=result_abs;
873  *resasc=result_asc;
874  *abserr=rescale_error(err,result_abs,result_asc);
875 
876  return;
877  }
878 
879  /** \brief Integration wrapper for user-specified function type
880  */
881  virtual void gauss_kronrod
882  (func_t &func, double a, double b,
883  double *result, double *abserr, double *resabs, double *resasc) {
884  return gauss_kronrod_base<func_t>
885  (func,a,b,result,abserr,resabs,resasc);
886 
887  }
888 
889  };
890 
891 #ifndef DOXYGEN_NO_O2NS
892 }
893 #endif
894 
895 #endif
static const double qk51_xgk[26]
Abscissae of the 51-point Kronrod rule.
int set_limit(size_t lim)
Set the limit for the number of subdivisions of the integration region (default 1000) ...
double * f_v1
Scratch space.
static const double qk41_wg[11]
Weights of the 20-point Gauss rule.
Integration workspace for the GSL integrators.
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
double * f_v2
Scratch space.
double * alist
Left endpoints of subintervals.
double * blist
Right endpoints of subintervals.
static const double qk31_wg[8]
Weights of the 15-point Gauss rule.
static const double qk41_wgk[21]
Weights of the 41-point Kronrod rule.
size_t * level
Numbers of subdivisions made.
invalid argument supplied by user
Definition: err_hnd.h:59
int n_gk
Size of Gauss-Kronrod arrays.
static const double qk21_wgk[11]
Weights of the 21-point Kronrod rule.
size_t * order
Linear ordering vector for sort routine.
static const double qk41_xgk[21]
Abscissae of the 41-point Kronrod rule.
size_t i
Index of current subinterval.
static const double qk61_xgk[31]
Abscissae of the 61-point Kronrod rule.
int get_rule()
Get the Gauss-Kronrod integration rule.
static const double qk21_xgk[11]
Abscissae of the 21-point Kronrod rule.
static const double qk15_wgk[8]
Weights of the 15-point Kronrod rule.
static const double qk61_wg[15]
Weights of the 30-point Gauss rule.
Base integration class [abstract base].
Definition: inte.h:44
double * rlist
Integral approximations on subintervals.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
A namespace for the GSL adaptive Gauss-Kronrod integration coefficients.
int allocate(size_t sz)
Allocate a workspace.
static const double qk61_wgk[31]
Weights of the 61-point Kronrod rule.
static const double qk31_xgk[16]
Abscissae of the 31-point Kronrod rule.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
size_t size
Current number of subintervals being used.
static const double qk15_wg[4]
Weights of the 7-point Gauss rule.
const double * w_gk
Gauss-Kronrod weight pointer.
static const double qk15_xgk[8]
Abscissae of the 15-point Kronrod rule.
GSL integration base.
Definition: inte_gsl.h:62
double * elist
Integral error estimates.
Basic Gauss-Kronrod integration class (GSL)
size_t nrmax
Counter for extrapolation routine.
static const double qk51_wgk[26]
Weights of the 51-point Kronrod rule.
inte_workspace_gsl * w
The integration workspace.
void set_rule(int rule)
Set the Gauss-Kronrod integration rule to be used.
static const double qk21_wg[5]
Weights of the 10-point Gauss rule.
static const double qk31_wgk[16]
Weights of the 31-point Kronrod rule.
size_t limit
Maximum number of subintervals allocated.
static const double qk51_wg[13]
Weights of the 25-point Gauss rule.
size_t maximum_level
Depth of subdivisions reached.
const double * x_gk
Gauss-Kronrod abscissae pointer.
const double * w_g
Gauss weight pointer.
int free()
Free allocated workspace memory.

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).