55 #ifndef O2SCL_CUBATURE_H 56 #define O2SCL_CUBATURE_H 60 #include <boost/numeric/ublas/vector.hpp> 62 #ifndef O2SCL_CLENCURT_H 63 #define O2SCL_CLENCURT_H 64 #include <o2scl/clencurt.h> 66 #include <o2scl/err_hnd.h> 67 #include <o2scl/vector.h> 132 template<
class func_t>
138 typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
169 double errMax(
const std::vector<esterr> &ee) {
171 for (
size_t k = 0; k < ee.size(); ++k) {
172 if (ee[k].err > errmax) errmax = ee[k].err;
210 for (
size_t i = 0; i < h.
dim; ++i) {
218 template<
class vec_t>
223 h.
data.resize(dim*2);
225 for (
size_t i = 0; i < dim; ++i) {
226 h.
data[i] = center[i];
227 h.
data[i + dim] = halfwidth[i];
229 h.
vol = compute_vol(h);
239 h.
data.resize(dim*2);
241 for (
size_t i = 0; i < dim; ++i) {
243 h.
data[i + dim] = dat[i+dim];
245 h.
vol = compute_vol(h);
251 template<
class vec_t>
252 void make_hypercube_range
253 (
size_t dim,
const vec_t &xmin,
const vec_t &xmax,
hypercube &h) {
255 make_hypercube(dim,xmin,xmax,h);
256 for (
size_t i = 0; i < dim; ++i) {
257 h.
data[i] = 0.5 * (xmin[i] + xmax[i]);
258 h.
data[i + dim] = 0.5 * (xmax[i] - xmin[i]);
260 h.
vol = compute_vol(h);
305 std::vector<esterr>
ee;
329 R.
h.
data[d + dim] *= 0.5;
331 make_hypercube2(dim, R.
h.
data,R2.
h);
395 r.
pts.resize((num_regions
427 if (rule15gauss_evalError(r, R[0].fdim, f, nR, R)) {
431 if (rule75genzmalik_evalError(r, R[0].fdim, f, nR, R)) {
435 for (iR = 0; iR < nR; ++iR) {
436 R[iR].errmax = errMax(R[iR].ee);
450 #if defined(__GNUC__) && \ 451 ((__GNUC__ == 3 && __GNUC_MINOR__ >= 4) || __GNUC__ > 3) 452 return __builtin_ctz(~n);
454 const size_t bits[256] = {
455 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
456 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
457 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
458 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
459 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
460 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
461 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
462 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 7,
463 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
464 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
465 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
466 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
467 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
468 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
469 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
470 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 8,
473 while ((n & 0xff) == 0xff) {
477 return bit + bits[n & 0xff];
489 std::vector<double> &p,
size_t p_ix,
490 const std::vector<double> &c,
size_t c_ix,
491 const std::vector<double> &r,
size_t r_ix) {
499 for (i = 0; i < dim; ++i) {
500 p[p_ix+i] = c[c_ix+i] + r[r_ix+i];
507 for(
size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
518 p[p_ix+d] = (signs & mask) ? c[c_ix+d] - r[r_ix+d] :
519 c[c_ix+d] + r[r_ix+d];
527 std::vector<double> &p,
size_t p_ix,
528 const std::vector<double> &c,
size_t c_ix,
529 const std::vector<double> &r,
size_t r_ix) {
531 for (
size_t i = 0; i < dim - 1; ++i) {
532 p[p_ix+i] = c[c_ix+i] - r[r_ix+i];
533 for (
size_t j = i + 1; j < dim; ++j) {
534 p[p_ix+j] = c[c_ix+j] - r[r_ix+j];
535 for(
size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
537 p[p_ix+i] = c[c_ix+i] + r[r_ix+i];
538 for(
size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
540 p[p_ix+j] = c[c_ix+j] + r[r_ix+j];
541 for(
size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
543 p[p_ix+i] = c[c_ix+i] - r[r_ix+i];
544 for(
size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
548 p[p_ix+j] = c[c_ix+j];
551 p[p_ix+i] = c[c_ix+i];
559 (ubvector &pts,
size_t pts_ix,
size_t dim,
560 std::vector<double> &p,
size_t p_ix,
561 const std::vector<double> &c,
size_t c_ix,
562 const std::vector<double> &r1,
size_t r1_ix,
563 const std::vector<double> &r2,
size_t r2_ix) {
565 for(
size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
568 for (
size_t i = 0; i < dim; i++) {
569 p[p_ix+i] = c[c_ix+i] - r1[r1_ix+i];
570 for(
size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
573 p[p_ix+i] = c[c_ix+i] + r1[r1_ix+i];
574 for(
size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
577 p[p_ix+i] = c[c_ix+i] - r2[r2_ix+i];
578 for(
size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
581 p[p_ix+i] = c[c_ix+i] + r2[r2_ix+i];
582 for(
size_t k=0;k<dim;k++) pts[pts_ix+k]=p[p_ix+k];
585 p[p_ix+i] = c[c_ix+i];
647 std::vector<double>
p;
661 #ifdef O2SCL_NEVER_DEFINED 667 int rule75genzmalik_evalError
668 (
rule &runder,
size_t fdim, func_t &f,
size_t nR,
669 std::vector<region> &R) {
672 const double lambda2 = 0.3585685828003180919906451539079374954541;
673 const double lambda4 = 0.9486832980505137995996680633298155601160;
674 const double lambda5 = 0.6882472016116852977216287342936235251269;
675 const double weight2 = 980. / 6561.;
676 const double weight4 = 200. / 19683.;
677 const double weightE2 = 245. / 486.;
678 const double weightE4 = 25. / 729.;
679 const double ratio = (lambda2 * lambda2) / (lambda4 * lambda4);
682 size_t i, j, iR, dim = runder.
dim;
685 ubvector &pts2=runder.
pts;
687 alloc_rule_pts(runder, nR);
690 for (iR = 0; iR < nR; ++iR) {
691 std::vector<double> ¢er2=R[iR].h.data;
693 for (i = 0; i < dim; ++i) {
694 r->
p[i] = center2[i];
697 for (i = 0; i < dim; ++i) {
698 r->
p[i+2*dim] = center2[i+dim] * lambda2;
700 for (i = 0; i < dim; ++i) {
701 r->
p[i+dim] = center2[i+dim] * lambda4;
706 evalR0_0fs4d2(runder.
pts,npts*dim, dim, r->
p,0,center2,0,
707 r->
p,2*dim,r->
p,dim);
708 npts += num0_0(dim) + 2 * numR0_0fs(dim);
711 evalRR0_0fs2(runder.
pts,npts*dim,dim,r->
p,0,center2,0,r->
p,dim);
712 npts += numRR0_0fs(dim);
715 for (i = 0; i < dim; ++i) {
716 r->
p[i+dim] = center2[i+dim] * lambda5;
718 evalR_Rfs2(runder.
pts,npts*dim,dim,r->
p,0,center2,0,r->
p,dim);
719 npts += numR_Rfs(dim);
723 if (f(dim, npts, &(pts2[0]), fdim, vals)) {
731 for (i = 0; i < dim * nR; ++i) pts2[i] = 0;
733 for (j = 0; j < fdim; ++j) {
735 const double *v = vals + j;
737 for (iR = 0; iR < nR; ++iR) {
738 double result, res5th;
739 double val0, sum2=0, sum3=0, sum4=0, sum5=0;
749 for (k = 0; k < dim; ++k) {
750 double v0 = v[fdim*(k0 + 4*k)];
751 double v1 = v[fdim*((k0 + 4*k) + 1)];
752 double v2 = v[fdim*((k0 + 4*k) + 2)];
753 double v3 = v[fdim*((k0 + 4*k) + 3)];
758 pts2[iR * dim + k] +=
759 fabs(v0 + v1 - 2*val0 - ratio * (v2 + v3 - 2*val0));
761 #ifdef O2SCL_NEVER_DEFINED 766 for (k = 0; k < numRR0_0fs(dim); ++k) {
767 sum4 += v[fdim*(k0 + k)];
771 for (k = 0; k < numR_Rfs(dim); ++k) {
772 sum5 += v[fdim*(k0 + k)];
776 result = R[iR].h.vol * (r->
weight1 * val0 + weight2 * sum2 +
777 r->
weight3 * sum3 + weight4 * sum4 +
779 res5th = R[iR].h.vol * (r->
weightE1 * val0 + weightE2 * sum2 +
780 r->
weightE3 * sum3 + weightE4 * sum4);
782 R[iR].ee[j].val = result;
783 R[iR].ee[j].err = fabs(res5th - result);
790 for (iR = 0; iR < nR; ++iR) {
792 size_t dimDiffMax = 0;
794 for (i = 0; i < dim; ++i) {
795 if (pts2[iR*dim + i] > maxdiff) {
796 maxdiff = pts2[iR*dim + i];
800 R[iR].splitDim = dimDiffMax;
805 #ifdef O2SCL_NEVER_DEFINED 815 O2SCL_ERR(
"this rule does not support 1d integrals",
824 if (dim >=
sizeof(
size_t) * 8) {
825 O2SCL_ERR(
"this rule does not support large dims",
830 make_rule(dim, fdim,num0_0(dim) + 2 * numR0_0fs(dim)
831 + numRR0_0fs(dim) + numR_Rfs(dim),r);
833 r.
weight1=(12824.0-9120.0*dim+400.0*dim*dim)/19683.0;
834 r.
weight3=(1820.0-400.0*dim)/19683.0;
835 r.
weight5=6859.0/19683.0/((double)(1U << dim));
836 r.
weightE1=(729.0-950.0*dim+50.0*dim*dim)/729.0;
837 r.
weightE3=(265.0-100.0*dim)/1458.0;
847 int rule15gauss_evalError
848 (
rule &r,
size_t fdim, func_t &f,
size_t nR, std::vector<region> &R) {
850 static const double cub_dbl_min=std::numeric_limits<double>::min();
851 static const double cub_dbl_eps=std::numeric_limits<double>::epsilon();
857 const double xgk[8] = {
858 0.991455371120812639206854697526329,
859 0.949107912342758524526189684047851,
860 0.864864423359769072789712788640926,
861 0.741531185599394439863864773280788,
862 0.586087235467691130294144838258730,
863 0.405845151377397166906606412076961,
864 0.207784955007898467600689403773245,
865 0.000000000000000000000000000000000
869 static const double wg[4] = {
870 0.129484966168869693270611432679082,
871 0.279705391489276667901467771423780,
872 0.381830050505118944950369775488975,
873 0.417959183673469387755102040816327
875 static const double wgk[8] = {
876 0.022935322010529224963732008058970,
877 0.063092092629978553290700663189204,
878 0.104790010322250183839876322541518,
879 0.140653259715525918745189590510238,
880 0.169004726639267902826583426598550,
881 0.190350578064785409913256402421014,
882 0.204432940075298892414161999234649,
883 0.209482141084727828012999174891714
889 alloc_rule_pts(r, nR);
893 for (iR = 0; iR < nR; ++iR) {
894 const double center = R[iR].h.data[0];
895 const double halfwidth = R[iR].h.data[1];
897 pts[npts++] = center;
899 for (j = 0; j < (n - 1) / 2; ++j) {
901 double w = halfwidth * xgk[j2];
902 pts[npts++] = center - w;
903 pts[npts++] = center + w;
905 for (j = 0; j < n/2; ++j) {
907 double w = halfwidth * xgk[j2];
908 pts[npts++] = center - w;
909 pts[npts++] = center + w;
915 if (f(1, npts, pts, fdim, vals)) {
919 for (k = 0; k < fdim; ++k) {
920 const double *vk = vals + k;
921 for (iR = 0; iR < nR; ++iR) {
922 const double halfwidth = R[iR].h.data[1];
923 double result_gauss = vk[0] * wg[n/2 - 1];
924 double result_kronrod = vk[0] * wgk[n - 1];
925 double result_abs = fabs(result_kronrod);
926 double result_asc, mean, err;
930 for (j = 0; j < (n - 1) / 2; ++j) {
932 double v = vk[fdim*npts] + vk[fdim*npts+fdim];
933 result_gauss += wg[j] * v;
934 result_kronrod += wgk[j2] * v;
935 result_abs += wgk[j2] * (fabs(vk[fdim*npts])
936 + fabs(vk[fdim*npts+fdim]));
939 for (j = 0; j < n/2; ++j) {
941 result_kronrod += wgk[j2] * (vk[fdim*npts]
942 + vk[fdim*npts+fdim]);
943 result_abs += wgk[j2] * (fabs(vk[fdim*npts])
944 + fabs(vk[fdim*npts+fdim]));
949 R[iR].ee[k].val = result_kronrod * halfwidth;
955 mean = result_kronrod * 0.5;
956 result_asc = wgk[n - 1] * fabs(vk[0] - mean);
958 for (j = 0; j < (n - 1) / 2; ++j) {
960 result_asc += wgk[j2] * (fabs(vk[fdim*npts]-mean)
961 + fabs(vk[fdim*npts+fdim]-mean));
964 for (j = 0; j < n/2; ++j) {
966 result_asc += wgk[j2] * (fabs(vk[fdim*npts]-mean)
967 + fabs(vk[fdim*npts+fdim]-mean));
970 err = fabs(result_kronrod - result_gauss) * halfwidth;
971 result_abs *= halfwidth;
972 result_asc *= halfwidth;
973 if (result_asc != 0 && err != 0) {
974 double scale = pow((200 * err / result_asc), 1.5);
975 err = (scale < 1) ? result_asc * scale : result_asc;
977 if (result_abs > cub_dbl_min / (50 * cub_dbl_eps)) {
978 double min_err = 50 * cub_dbl_eps * result_abs;
979 if (min_err > err) err = min_err;
981 R[iR].ee[k].err = err;
998 return make_rule(dim,fdim,15,r);
1055 h.
items.resize(nalloc);
1069 for (
size_t i = 0; i < fdim; ++i) {
1073 heap_resize(h, nalloc);
1093 size_t fdim = h.
fdim;
1095 for (
size_t i = 0; i < fdim; ++i) {
1096 h.
ee[i].val += hi.
ee[i].val;
1097 h.
ee[i].err += hi.
ee[i].err;
1101 heap_resize(h, h.
n * 2);
1105 int parent = (insert - 1) / 2;
1112 h.
items[insert] = hi;
1119 for (
size_t i = 0; i < ni; ++i) {
1133 O2SCL_ERR(
"Attempted to pop an empty heap in cubature.",
1140 while ((child = i * 2 + 1) < n) {
1145 if (h.
items[child].errmax <= h.
items[i].errmax) {
1151 if (++child < n && h.
items[largest].errmax <
1152 h.
items[child].errmax) {
1160 h.
items[i = largest] = swap;
1164 size_t i, fdim = h.
fdim;
1165 for (i = 0; i < fdim; ++i) {
1166 h.
ee[i].val -= ret.
ee[i].val;
1167 h.
ee[i].err -= ret.
ee[i].err;
1177 double reqAbsError,
double reqRelError,
1186 for (j = 0; j < fdim; ++j) {
1187 if (ee[j].err > reqAbsError && ee[j].err >
1188 fabs(ee[j].val)*reqRelError) {
1196 for (j = 0; j+1 < fdim; j += 2) {
1197 double maxerr, serr, err, maxval, sval, val;
1199 maxerr = ee[j].err > ee[j+1].err ? ee[j].err : ee[j+1].err;
1200 maxval = ee[j].val > ee[j+1].val ? ee[j].val : ee[j+1].val;
1201 serr = maxerr > 0 ? 1/maxerr : 1;
1202 sval = maxval > 0 ? 1/maxval : 1;
1203 err=std::hypot(ee[j].err*serr,ee[j+1].err*serr)*maxerr;
1204 val=std::hypot(ee[j].val*sval,ee[j+1].val*sval)*maxval;
1205 if (err > reqAbsError && err > val*reqRelError) {
1212 if (ee[j].err > reqAbsError && ee[j].err >
1213 fabs(ee[j].val)*reqRelError) {
1222 double err = 0, val = 0;
1223 for (j = 0; j < fdim; ++j) {
1225 val += fabs(ee[j].val);
1227 return err <= reqAbsError || err <= val*reqRelError;
1232 double err = 0, val = 0;
1233 for (j = 0; j < fdim; ++j) {
1234 double absval = fabs(ee[j].val);
1235 if (ee[j].err > err) err = ee[j].err;
1236 if (absval > val) val = absval;
1238 return err <= reqAbsError || err <= val*reqRelError;
1243 double maxerr = 0, maxval = 0, serr, sval, err = 0, val = 0;
1245 for (j = 0; j < fdim; ++j) {
1246 double absval = fabs(ee[j].val);
1247 if (ee[j].err > maxerr) maxerr = ee[j].err;
1248 if (absval > maxval) maxval = absval;
1250 serr = maxerr > 0 ? 1/maxerr : 1;
1251 sval = maxval > 0 ? 1/maxval : 1;
1252 for (j = 0; j < fdim; ++j) {
1253 err += (ee[j].err * serr)*(ee[j].err * serr);
1254 val += fabs(ee[j].val) * sval*fabs(ee[j].val) * sval;
1256 err = sqrt(err) * maxerr;
1257 val = sqrt(val) * maxval;
1258 return err <= reqAbsError || err <= val*reqRelError;
1267 template<
class vec_t>
1270 double reqAbsError,
double reqRelError,
1272 vec_t &err,
int parallel) {
1278 std::vector<region> R;
1279 size_t nR_alloc = 0;
1280 std::vector<esterr> ee(fdim);
1287 regions = heap_alloc(1, fdim);
1291 make_region(h, fdim, R[0]);
1292 if (eval_regions(1, R, f, r) || heap_push(regions, R[0])) {
1299 while (numEval < maxEval || !maxEval) {
1301 if (converged(fdim, regions.
ee, reqAbsError, reqRelError, norm)) {
1334 for (j = 0; j < fdim; ++j) ee[j] = regions.
ee[j];
1337 if (nR + 2 > nR_alloc) {
1338 nR_alloc = (nR + 2) * 2;
1341 R[nR] = heap_pop(regions);
1342 for (j = 0; j < fdim; ++j) ee[j].err -= R[nR].ee[j].err;
1345 if (converged(fdim, ee, reqAbsError, reqRelError, norm)) {
1349 }
while (regions.
n > 0 && (numEval < maxEval || !maxEval));
1351 if (eval_regions(nR, R, f, r)
1352 || heap_push_many(regions, nR, &(R[0]))) {
1364 R[0] = heap_pop(regions);
1365 if (cut_region(R[0], R[1]) || eval_regions(2, R, f, r)
1366 || heap_push_many(regions, 2, &(R[0]))) {
1376 for (j = 0; j < fdim; ++j) {
1380 for (i = 0; i < regions.n; ++i) {
1381 for (j = 0; j < fdim; ++j) {
1382 val[j] += regions.items[i].ee[j].val;
1383 err[j] += regions.items[i].ee[j].err;
1386 destroy_hypercube(regions.items[i].h);
1387 regions.items[i].ee.clear();
1399 template<
class vec_t>
1400 int cubature(
size_t fdim, func_t &f,
size_t dim,
const vec_t &xmin,
1401 const vec_t &xmax,
size_t maxEval,
double reqAbsError,
1402 double reqRelError,
error_norm norm, vec_t &val,
1403 vec_t &err,
int parallel) {
1413 if (f(0, 1, &(xmin[0]), fdim, &(val[0]))) {
1416 for (
size_t i = 0; i < fdim; ++i) err[i] = 0;
1422 make_rule15gauss(dim,fdim,r);
1423 make_hypercube_range(dim,xmin,xmax,h);
1424 status = rulecubature(r,fdim,f,h,maxEval,reqAbsError,
1425 reqRelError,norm,val,err,parallel);
1428 make_rule75genzmalik(dim,fdim,r);
1429 make_hypercube_range(dim,xmin,xmax,h);
1430 status = rulecubature(r,fdim,f,h,maxEval,reqAbsError,
1431 reqRelError,norm,val,err,parallel);
1433 destroy_hypercube(h);
1449 template<
class vec_t>
1450 int integ(
size_t fdim, func_t &f,
size_t dim,
1451 const vec_t &xmin,
const vec_t &xmax,
size_t maxEval,
1452 double reqAbsError,
double reqRelError,
error_norm norm,
1453 vec_t &val, vec_t &err) {
1459 return cubature(fdim,f,dim,xmin,xmax,maxEval,reqAbsError,
1460 reqRelError,norm,val,err,use_parallel);
1466 #ifdef O2SCL_NEVER_DEFINED 1485 template<
class func_t,
class vec_t,
class vec_crange_t,
class vec_range_t>
1492 static const size_t MAXDIM=20;
1527 std::vector<size_t>
m;
1541 size_t &vali,
size_t fdim, func_t &f,
size_t dim,
1542 size_t id, std::vector<double> &p,
const vec_t &xmin,
1543 const vec_t &xmax, vec_t &buf,
size_t nbuf,
1548 for(
size_t k=0;k<dim;k++) {
1549 buf[k+ibuf*dim]=p[k];
1554 (((
const vec_t &)buf),0,buf.size());
1557 if (f(dim, nbuf, &(buf[0]), fdim, &(val[0]) + vali)) {
1560 vali += ibuf * fdim;
1566 double c = (xmin[id] + xmax[id]) * 0.5;
1567 double r = (xmax[id] - xmin[id]) * 0.5;
1568 const double *x = clencurt_x
1569 + ((
id == mi) ? (m[
id] ? (1 << (m[
id] - 1)) : 0) : 0);
1571 size_t nx = (
id == mi ? (m[id] ? (1 << (m[id] - 1)) : 1)
1575 if (compute_cacheval(m, mi, val, vali, fdim, f,
1577 xmin, xmax, buf, nbuf, ibuf)) {
1581 for (i = 0; i < nx; ++i) {
1582 p[id] = c + r * x[i];
1583 if (compute_cacheval(m, mi, val, vali, fdim, f,
1585 xmin, xmax, buf, nbuf, ibuf)) {
1588 p[id] = c - r * x[i];
1589 if (compute_cacheval(m, mi, val, vali, fdim, f,
1591 xmin, xmax, buf, nbuf, ibuf)) {
1605 for (
size_t i = 0; i < dim; ++i) {
1607 if (m[i+i_shift]==0) {
1610 nval*= (1 << (m[i+i_shift]));
1613 nval *= (1 << (m[i+i_shift] + 1)) + 1;
1622 const std::vector<size_t> &m,
1623 size_t mi,
size_t fdim, func_t &f,
size_t dim,
1624 const vec_t &xmin,
const vec_t &xmax, vec_t &buf,
1627 size_t ic = vc.size();
1628 size_t nval, vali = 0, ibuf = 0;
1629 std::vector<double> p(MAXDIM);
1635 for(
size_t j=0;j<dim;j++) {
1638 nval = fdim * num_cacheval(m,mi,dim,0);
1639 vc[ic].val.resize(nval);
1641 if (compute_cacheval(m,mi,vc[ic].val,vali,fdim,f,
1642 dim,0,p,xmin,xmax,buf,nbuf,ibuf)) {
1649 (((
const vec_t &)buf),0,buf.size());
1651 return f(dim, ibuf, &(buf[0]), fdim, &((vc[ic].val)[vali]));
1665 size_t eval(
const std::vector<size_t> &cm,
size_t cmi,
1666 vec_t &cval,
const std::vector<size_t> &m,
size_t md,
1667 size_t fdim,
size_t dim,
size_t id,
1668 double weight, vec_t &val,
size_t voff2) {
1673 for (i = 0; i < fdim; ++i) {
1674 val[i] += cval[i+voff2] * weight;
1678 }
else if (m[
id] == 0 &&
id == md) {
1681 voff = eval(cm, cmi, cval, m, md, fdim, dim,
id+1, weight*2, val,voff2);
1682 voff += fdim * (1 << cm[id]) * 2
1683 * num_cacheval(cm, cmi - (
id+1), dim - (
id+1),
id+1);
1689 size_t mid = m[id] - (
id == md);
1690 const double *w = clencurt_w + mid + (1 << mid) - 1
1691 + (
id == cmi ? (cm[
id] ? 1 + (1 << (cm[
id]-1)) : 1) : 0);
1692 size_t cnx = (
id == cmi ? (cm[id] ? (1 << (cm[id]-1)) : 1)
1694 size_t nx = cm[id] <= mid ? cnx : (1 << mid);
1697 voff = eval(cm, cmi, cval, m, md, fdim, dim,
id + 1,
1698 weight * w[0], val,voff2);
1701 for (i = 0; i < nx; ++i) {
1702 voff += eval(cm, cmi, cval, m, md, fdim, dim,
id + 1,
1703 weight * w[i], val,voff+voff2);
1704 voff += eval(cm, cmi, cval, m, md, fdim, dim,
id + 1,
1705 weight * w[i], val,voff+voff2);
1708 voff += (cnx - nx) * fdim * 2
1709 * num_cacheval(cm, cmi - (
id+1), dim - (
id+1),
id+1);
1719 void evals(std::vector<cache> &vc,
const std::vector<size_t> &m,
1720 size_t md,
size_t fdim,
size_t dim,
double V, vec_t &val) {
1722 for(
size_t k=0;k<fdim;k++) {
1725 for (
size_t i = 0; i < vc.size(); ++i) {
1726 if (vc[i].mi >= dim ||
1727 vc[i].m[vc[i].mi] + (vc[i].mi == md) <= m[vc[i].mi]) {
1728 eval(vc[i].m, vc[i].mi, vc[i].val,
1729 m, md, fdim, dim, 0, V, val,0);
1743 const std::vector<size_t> &m,
1744 size_t fdim,
size_t dim,
double V,
1745 size_t &mi, vec_t &val, vec_t &err, vec_t &val1) {
1750 evals(vc, m, dim, fdim, dim, V, val);
1755 for(
size_t j=0;j<fdim;j++) {
1759 for (i = 0; i < dim; ++i) {
1761 evals(vc, m, i, fdim, dim, V, val1);
1762 for (j = 0; j < fdim; ++j) {
1763 double e = fabs(val[j] - val1[j]);
1764 if (e > emax) emax = e;
1765 if (e > err[j]) err[j] = e;
1767 if (emax > maxerr) {
1777 template<
class vec2_t>
1778 int converged(
size_t fdim,
const vec2_t &vals,
const vec2_t &errs,
1779 double reqAbsError,
double reqRelError,
error_norm norm) {
1785 for (
size_t j = 0; j < fdim; ++j) {
1786 if (errs[j] > reqAbsError && errs[j] > fabs(vals[j])*reqRelError) {
1797 for (j = 0; j+1 < fdim; j += 2) {
1798 double maxerr, serr, err, maxval, sval, val;
1800 maxerr = errs[j] > errs[j+1] ? errs[j] : errs[j+1];
1801 maxval = vals[j] > vals[j+1] ? vals[j] : vals[j+1];
1802 serr = maxerr > 0 ? 1/maxerr : 1;
1803 sval = maxval > 0 ? 1/maxval : 1;
1804 err=std::hypot(errs[j]*serr,errs[j+1]*serr)*maxerr;
1805 val=std::hypot(vals[j]*sval,vals[j+1]*sval)*maxval;
1806 if (err > reqAbsError && err > val*reqRelError) {
1812 if (errs[j] > reqAbsError && errs[j] > fabs(vals[j])*reqRelError) {
1822 double err = 0, val = 0;
1823 for (
size_t j = 0; j < fdim; ++j) {
1825 val += fabs(vals[j]);
1827 return err <= reqAbsError || err <= val*reqRelError;
1833 double err = 0, val = 0;
1834 for (
size_t j = 0; j < fdim; ++j) {
1835 double absval = fabs(vals[j]);
1836 if (errs[j] > err) err = errs[j];
1837 if (absval > val) val = absval;
1839 return err <= reqAbsError || err <= val*reqRelError;
1845 double maxerr = 0, maxval = 0, serr, sval, err = 0, val = 0;
1847 for (
size_t j = 0; j < fdim; ++j) {
1848 double absval = fabs(vals[j]);
1849 if (errs[j] > maxerr) maxerr = errs[j];
1850 if (absval > maxval) maxval = absval;
1852 serr = maxerr > 0 ? 1/maxerr : 1;
1853 sval = maxval > 0 ? 1/maxval : 1;
1854 for (
size_t j = 0; j < fdim; ++j) {
1855 err += (errs[j] * serr)*(errs[j] * serr);
1856 val += (fabs(vals[j]) * sval)*(fabs(vals[j]) * sval);
1858 err = sqrt(err) * maxerr;
1859 val = sqrt(val) * maxval;
1860 return err <= reqAbsError || err <= val*reqRelError;
1864 O2SCL_ERR(
"Improper value of 'norm' in cubature::converged().",
1884 int integ_v_buf(
size_t fdim, func_t &f,
size_t dim,
const vec_t &xmin,
1885 const vec_t &xmax,
size_t maxEval,
double reqAbsError,
1886 double reqRelError,
error_norm norm, std::vector<size_t> &m,
1887 vec_t &buf,
size_t &nbuf,
size_t max_nbuf, vec_t &val,
1892 size_t numEval = 0, new_nbuf;
1894 std::vector<cache> vc;
1912 for (i = 0; i < fdim; ++i) err[i] = 0;
1916 for (i = 0; i < fdim; ++i) {
1921 for (i = 0; i < dim; ++i) {
1923 V *= (xmax[i] - xmin[i]) * 0.5;
1926 new_nbuf = num_cacheval(m,dim,dim,0);
1928 if (max_nbuf < 1) max_nbuf = 1;
1929 if (new_nbuf > max_nbuf) new_nbuf = max_nbuf;
1930 if (nbuf < new_nbuf) {
1931 buf.resize((nbuf=new_nbuf)*dim);
1935 if (add_cacheval(vc,m, dim, fdim, f, dim, xmin, xmax,
1942 eval_integral(vc,m,fdim,dim,V,mi,val,err,val1);
1944 if (converged(fdim,val,err,reqAbsError, reqRelError, norm) ||
1945 (numEval > maxEval && maxEval)) {
1950 if (m[mi] > clencurt_M) {
1955 new_nbuf = num_cacheval(m, mi, dim,0);
1956 if (new_nbuf > nbuf && nbuf < max_nbuf) {
1958 if (nbuf > max_nbuf) nbuf = max_nbuf;
1959 buf.resize(nbuf*dim);
1962 if (add_cacheval(vc,m, mi, fdim, f,
1967 numEval += new_nbuf;
1975 static const size_t DEFAULT_MAX_NBUF=(1U << 20);
1979 int integ(
size_t fdim, func_t &f,
size_t dim,
1980 const vec_t &xmin,
const vec_t &xmax,
size_t maxEval,
1981 double reqAbsError,
double reqRelError,
error_norm norm,
1982 vec_t &val, vec_t &err) {
1986 std::vector<size_t> m(dim);
1990 ret = integ_v_buf(fdim,f,dim,xmin,xmax,
1991 maxEval,reqAbsError,reqRelError,norm,
1992 m,buf,nbuf,16,val,err);
size_t num_points
The number of evaluation points.
size_t eval(const std::vector< size_t > &cm, size_t cmi, vec_t &cval, const std::vector< size_t > &m, size_t md, size_t fdim, size_t dim, size_t id, double weight, vec_t &val, size_t voff2)
Desc.
Base class for integration routines from the Cubature library.
size_t num0_0(size_t dim)
Desc.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
double vol
cache volume = product of widths
int rulecubature(rule &r, size_t fdim, func_t &f, const hypercube &h, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err, int parallel)
Desc.
int heap_push(heap &h, heap_item &hi)
Desc.
int cubature(size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err, int parallel)
Desc.
sanity check failed - shouldn't happen
std::vector< double > p
Desc.
invalid argument supplied by user
void make_rule15gauss(size_t dim, size_t fdim, rule &r)
Desc.
void evals(std::vector< cache > &vc, const std::vector< size_t > &m, size_t md, size_t fdim, size_t dim, double V, vec_t &val)
Desc.
void make_rule(size_t dim, size_t fdim, size_t num_points, rule &r)
Desc.
size_t num_regions
The max number of regions evaluated at once.
size_t numR0_0fs(size_t dim)
Desc.
std::vector< heap_item > items
Desc.
int converged(size_t fdim, const std::vector< esterr > &ee, double reqAbsError, double reqRelError, error_norm norm)
Desc.
size_t vals_ix
num_regions * num_points * fdim
void eval_integral(std::vector< cache > &vc, const std::vector< size_t > &m, size_t fdim, size_t dim, double V, size_t &mi, vec_t &val, vec_t &err, vec_t &val1)
Desc.
size_t fdim
The number of functions.
int integ(size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err)
Desc.
int heap_push_many(heap &h, size_t ni, heap_item *hi)
Desc.
double weight1
dimension-dependent constants
void make_region(const hypercube &h, size_t fdim, region &R)
Desc.
void heap_free(heap &h)
Note that heap_free does not deallocate anything referenced by the items.
int integ_v_buf(size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, std::vector< size_t > &m, vec_t &buf, size_t &nbuf, size_t max_nbuf, vec_t &val, vec_t &err)
Desc.
ubvector pts
points to eval: num_regions * num_points * dim
individual relerr criteria in each component
int compute_cacheval(const std::vector< size_t > &m, size_t mi, vec_t &val, size_t &vali, size_t fdim, func_t &f, size_t dim, size_t id, std::vector< double > &p, const vec_t &xmin, const vec_t &xmax, vec_t &buf, size_t nbuf, size_t &ibuf)
Desc.
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
void destroy_hypercube(hypercube &h)
Desc.
error_norm
Different ways of measuring the absolute and relative error.
void make_hypercube2(size_t dim, const std::vector< double > &dat, hypercube &h)
Desc.
std::vector< size_t > m
Desc.
Cache of the values for the m[dim] grid.
void evalRR0_0fs2(ubvector &pts, size_t pts_ix, size_t dim, std::vector< double > &p, size_t p_ix, const std::vector< double > &c, size_t c_ix, const std::vector< double > &r, size_t r_ix)
Desc.
int eval_regions(size_t nR, std::vector< region > &R, func_t &f, rule &r)
Desc.
size_t numR_Rfs(size_t dim)
Desc.
void evalR_Rfs2(ubvector &pts, size_t pts_ix, size_t dim, std::vector< double > &p, size_t p_ix, const std::vector< double > &c, size_t c_ix, const std::vector< double > &r, size_t r_ix)
Evaluate the integration points for all points (+/-r,...+/-r)
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
abserr is L_2 norm |e|, and relerr is |e|/|v|
abserr is L_1 norm |e|, and relerr is |e|/|v|
Adaptive multidimensional integration on hyper-rectangles using cubature rules from the Cubature libr...
void make_hypercube(size_t dim, const vec_t ¢er, const vec_t &halfwidth, hypercube &h)
Desc.
double errMax(const std::vector< esterr > &ee)
Return the maximum error from ee.
heap heap_alloc(size_t nalloc, size_t fdim)
Desc.
int integ(size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err)
Desc.
Integration by p-adaptive cubature from the Cubature library.
double errmax
max ee[k].err
abserr is norm |e|, and relerr is |e|/|v|
heap_item heap_pop(heap &h)
Desc.
std::vector< double > data
length 2*dim = center followed by half-widths
std::vector< esterr > ee
array of length fdim
void alloc_rule_pts(rule &r, size_t num_regions)
Desc.
size_t fdim
dimensionality of vector integrand
size_t numRR0_0fs(size_t dim)
Desc.
size_t num_cacheval(const std::vector< size_t > &m, size_t mi, size_t dim, size_t i_shift)
Desc.
size_t ls0(size_t n)
Functions to loop over points in a hypercube.
int add_cacheval(std::vector< cache > &vc, const std::vector< size_t > &m, size_t mi, size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, vec_t &buf, size_t nbuf)
Desc.
int converged(size_t fdim, const vec2_t &vals, const vec2_t &errs, double reqAbsError, double reqRelError, error_norm norm)
Desc.
double compute_vol(const hypercube &h)
Desc.
void make_rule75genzmalik(size_t dim, size_t fdim, rule75genzmalik &r)
Desc.
const dat_t * const_vector_range(const dat_t *v, size_t start, size_t last)
Vector range function for const pointers.
void heap_resize(heap &h, size_t nalloc)
Desc.
paired L2 norms of errors in each component, mainly for integrating vectors of complex numbers ...
size_t dim
The dimensionality.
int cut_region(region &R, region &R2)
Desc.