rule75genzmalik_evalError()

static int rule75genzmalik_evalError ( rule *  r_, unsigned  fdim, integrand_v  f, void *  fdata, unsigned  nR, region *  R  )

static

{
     /* lambda2 = sqrt(9/70), lambda4 = sqrt(9/10), lambda5 = sqrt(9/19) */
     const double lambda2 = 0.3585685828003180919906451539079374954541;
     const double lambda4 = 0.9486832980505137995996680633298155601160;
     const double lambda5 = 0.6882472016116852977216287342936235251269;
     const double weight2 = 980. / 6561.;
     const double weight4 = 200. / 19683.;
     const double weightE2 = 245. / 486.;
     const double weightE4 = 25. / 729.;
     const double ratio = (lambda2 * lambda2) / (lambda4 * lambda4);

     rule75genzmalik *r = (rule75genzmalik *) r_;
     unsigned i, j, iR, dim = r_->dim, npts = 0;
     double *diff, *pts, *vals;

     if (alloc_rule_pts(r_, nR)) return FAILURE;
     pts = r_->pts; vals = r_->vals;

     for (iR = 0; iR < nR; ++iR) {
          const double *center = R[iR].h.data;
          const double *halfwidth = R[iR].h.data + dim;
          
          for (i = 0; i < dim; ++i)
               r->p[i] = center[i];
          
          for (i = 0; i < dim; ++i)
               r->widthLambda2[i] = halfwidth[i] * lambda2;
          for (i = 0; i < dim; ++i)
               r->widthLambda[i] = halfwidth[i] * lambda4;

          /* Evaluate points in the center, in (lambda2,0,...,0) and
             (lambda3=lambda4, 0,...,0).  */
          evalR0_0fs4d(pts + npts*dim, dim, r->p, center, 
                       r->widthLambda2, r->widthLambda);
          npts += num0_0(dim) + 2 * numR0_0fs(dim);

          /* Calculate points for (lambda4, lambda4, 0, ...,0) */
          evalRR0_0fs(pts + npts*dim, dim, r->p, center, r->widthLambda);
          npts += numRR0_0fs(dim);

          /* Calculate points for (lambda5, lambda5, ..., lambda5) */
          for (i = 0; i < dim; ++i)
               r->widthLambda[i] = halfwidth[i] * lambda5;
          evalR_Rfs(pts + npts*dim, dim, r->p, center, r->widthLambda);
          npts += numR_Rfs(dim);
     }

     /* Evaluate the integrand function(s) at all the points */
     f(dim, npts, pts, fdata, fdim, vals);

     /* we are done with the points, and so we can re-use the pts
        array to store the maximum difference diff[i] in each dimension 
        for each hypercube */
     diff = pts;
     for (i = 0; i < dim * nR; ++i) diff[i] = 0;

     for (j = 0; j < fdim; ++j) {
          for (iR = 0; iR < nR; ++iR) {
               double result, res5th;
               double val0, sum2=0, sum3=0, sum4=0, sum5=0;
               unsigned k, k0 = 0;
               
               /* accumulate j-th function values into j-th integrals
                     NOTE: this relies on the ordering of the eval functions
                     above, as well as on the internal structure of
                     the evalR0_0fs4d function */

               val0 = vals[0]; /* central point */
               k0 += 1;

               for (k = 0; k < dim; ++k) {
                    double v0 = vals[k0 + 4*k];
                    double v1 = vals[(k0 + 4*k) + 1];
                    double v2 = vals[(k0 + 4*k) + 2];
                    double v3 = vals[(k0 + 4*k) + 3];
                    
                    sum2 += v0 + v1;
                    sum3 += v2 + v3;
                    
                    diff[iR * dim + k] += 
                         fabs(v0 + v1 - 2*val0 - ratio * (v2 + v3 - 2*val0));
               }
               k0 += 4*k;

               for (k = 0; k < numRR0_0fs(dim); ++k)
                    sum4 += vals[k0 + k];
               k0 += k;
               
               for (k = 0; k < numR_Rfs(dim); ++k)
                    sum5 += vals[k0 + k];
               
               /* Calculate fifth and seventh order results */
               result = R[iR].h.vol * (r->weight1 * val0 + weight2 * sum2 + r->weight3 * sum3 + weight4 * sum4 + r->weight5 * sum5);
               res5th = R[iR].h.vol * (r->weightE1 * val0 + weightE2 * sum2 + r->weightE3 * sum3 + weightE4 * sum4);
               
               R[iR].ee[j].val = result;
               R[iR].ee[j].err = fabs(res5th - result);
               
               vals += r_->num_points;
          }
     }

     /* figure out dimension to split: */
     for (iR = 0; iR < nR; ++iR) {
          double maxdiff = 0;
          unsigned dimDiffMax = 0;

          for (i = 0; i < dim; ++i)
               if (diff[iR*dim + i] > maxdiff) {
                    maxdiff = diff[iR*dim + i];
                    dimDiffMax = i;
               }
          R[iR].splitDim = dimDiffMax;
     }
     return SUCCESS;
}