GKIntegration< F, GKRule > Class Template Reference

Collaboration diagram for GKIntegration< F, GKRule >:

Classes
class	IntervalResult

Public Member Functions
	GKIntegration (F &f_)
void	SetAbsoluteErrorMode ()
void	SetRelativeErrorMode ()
double	Integrate (const double a, const double b, const double acc)
double	Integrate (const double a, const double b, const double accAbs, const double accRel, const bool andFlag)

Private Member Functions
void	GK (IntervalResult &r)
void	GKGeneral (const int n, const double xgk[], const double wg[], const double wgk[], IntervalResult &r)
void	PrintList ()
void	PrintMax ()
void	CheckList ()
void	CheckError (const double err)
double	RecomputeError ()
void	Insert (const IntervalResult &r)
double	Integrate (const double a, const double b, const double absError, const bool absErrorFlag, const double relError, const bool relErrorFlag, const bool andFlag)

Private Attributes
std::list< IntervalResult >	ir
F &	f
bool	relativeErrors

Detailed Description

template<class F, class GKRule = GK31>
class GKIntegration< F, GKRule >

Definition at line 90 of file GKIntegration.h.

Constructor & Destructor Documentation

GKIntegration()

GKIntegration ( F &  f_ )

inline

Definition at line 115 of file GKIntegration.h.

: f(f_), relativeErrors(false) {}

Member Function Documentation

CheckError()

void CheckError ( const double  err )

inlineprivate

Definition at line 215 of file GKIntegration.h.

References qmcplusplus::abs(), BMWrite4, qmcplusplus::Units::charge::e, error(), and GKIntegration< F, GKRule >::ir.

  {
    double errorSum = 0.0;

    if (ir.size() > 0)
    {
      for (typename std::list<IntervalResult>::iterator p = ir.begin(); p != ir.end(); p++)
      {
        errorSum += (*p).err;
      }
    }

    if (errorSum == 0.0)
    {
      if (err != 0.0)
        error("CheckError", errorSum, err);
    }
    else
    {
      if (err / errorSum - 1.0 > 1e-8 && std::abs(err - errorSum) > 1e-14)
        error("CheckError", errorSum, err, errorSum - err);
    }

    BMWrite4("PassedErrorCheck", errorSum, err, errorSum - err);
  }

CheckList()

void CheckList ( )

inlineprivate

Definition at line 196 of file GKIntegration.h.

References BMWrite2, error(), GKIntegration< F, GKRule >::ir, and GKIntegration< F, GKRule >::PrintList().

  {
    if (ir.size() == 0)
      return;

    for (typename std::list<IntervalResult>::iterator p = ir.begin(); p != ir.end(); p++)
    {
      typename std::list<IntervalResult>::iterator pn = p;
      pn++;
      if (pn != ir.end())
        if (((*p).err) < ((*pn).err))
        {
          PrintList();
          BMWrite2(*pn, *p);
          ::error("Ordering problem in list");
        }
    }
  }

GK()

void GK ( IntervalResult &  r )

inlineprivate

Definition at line 122 of file GKIntegration.h.

References GKIntegration< F, GKRule >::GKGeneral(), and qmcplusplus::n.

Referenced by GKIntegration< F, GKRule >::Integrate().

{ GKGeneral(GKRule::n, GKRule::xgk, GKRule::wg, GKRule::wgk, r); }

GKGeneral()

void GKGeneral ( const int  n, const double  xgk[], const double  wg[], const double  wgk[], IntervalResult &  r  )

inlineprivate

Definition at line 125 of file GKIntegration.h.

References GKIntegration< F, GKRule >::IntervalResult::a, qmcplusplus::abs(), GKIntegration< F, GKRule >::IntervalResult::b, GKIntegration< F, GKRule >::IntervalResult::delta, GKIntegration< F, GKRule >::IntervalResult::err, GKIntegration< F, GKRule >::f, qmcplusplus::n, and GKIntegration< F, GKRule >::IntervalResult::result.

Referenced by GKIntegration< F, GKRule >::GK().

  {
    const double center     = 0.5 * (r.a + r.b);
    const double halfLength = 0.5 * r.delta;
    const double fCenter    = f(center);

    double resultGauss   = 0;
    double resultKronrod = fCenter * wgk[n - 1];

    if (n % 2 == 0)
    {
      resultGauss = fCenter * wg[n / 2 - 1];
    }

    for (int j = 0; j < (n - 1) / 2; j++)
    {
      const int jtw      = j * 2 + 1; // j=1,2,3 jtw=2,4,6
      const double xx    = halfLength * xgk[jtw];
      const double fval1 = f(center - xx);
      const double fval2 = f(center + xx);
      const double fsum  = fval1 + fval2;
      resultGauss   += wg[j] * fsum;
      resultKronrod += wgk[jtw] * fsum;
    }

    for (int j = 0; j < n / 2; j++)
    {
      int jtwm1          = j * 2;
      const double xx    = halfLength * xgk[jtwm1];
      const double fval1 = f(center - xx);
      const double fval2 = f(center + xx);
      resultKronrod += wgk[jtwm1] * (fval1 + fval2);
    };

    /* scale by the width of the integration region */
    resultGauss   *= halfLength;
    resultKronrod *= halfLength;

    r.result = resultKronrod;
    r.err    = std::abs(resultKronrod - resultGauss);
    //r.err = pow(200.0 * std::abs(resultKronrod - resultGauss), 1.5);
    //    BMWrite(r);
  }

Insert()

void Insert ( const IntervalResult &  r )

inlineprivate

Definition at line 256 of file GKIntegration.h.

References GKIntegration< F, GKRule >::IntervalResult::err, and GKIntegration< F, GKRule >::ir.

Referenced by GKIntegration< F, GKRule >::Integrate().

  {
    //    std::cout << "Inserting.." << std::endl;
    //    PrintList();

    if (ir.empty())
    {
      ir.push_back(r);
      return;
    }

    if (ir.back().err >= r.err)
    {
      ir.push_back(r);
      return;
    }

    if (r.err >= ir.front().err)
    {
      ir.push_front(r);
      return;
    }

    // size must be >=2
    typename std::list<IntervalResult>::iterator p = ir.end();
    p--;

    // p cannot become less the begin() because of check above
    while (r.err > (*p).err)
      p--;

    // go one down because insert put the element before p
    p++;
    ir.insert(p, r);
    //    CheckList();
  }

Integrate() [1/3]

double Integrate ( const double  a, const double  b, const double  absError, const bool  absErrorFlag, const double  relError, const bool  relErrorFlag, const bool  andFlag  )

inlineprivate

Definition at line 293 of file GKIntegration.h.

References qmcplusplus::abs(), GKIntegration< F, GKRule >::IntervalResult::delta, qmcplusplus::Units::charge::e, GKIntegration< F, GKRule >::IntervalResult::err, GKIntegration< F, GKRule >::GK(), GKIntegration< F, GKRule >::Insert(), GKIntegration< F, GKRule >::ir, qmcplusplus::ldexp(), GKIntegration< F, GKRule >::PrintMax(), GKIntegration< F, GKRule >::RecomputeError(), GKIntegration< F, GKRule >::IntervalResult::result, and warning().

Referenced by GKIntegration< F, GKRule >::Integrate().

  {
    ir.clear();

    // #define PRINT_IT
#ifdef PRINT_IT
    std::cout << "Beginning integration" << std::endl;
#endif

    double errorUnresolved = 0.0;
    const int iterationMax = 30;
    // double lengthMin = (b-a)*pow(0.5,iterationMax);
    double lengthMin = ldexp(b - a, -iterationMax);

    IntervalResult r0(a, b, b - a);
    GK(r0);
    double result  = r0.result;
    double err     = r0.err;
    double errLast = err;

    ir.push_back(r0);

    bool quitFlag;
    do
    {
      // PrintList();

      // Test if the interval with the biggest error has already been subdivided
      // the maximum number of times. If this is the case throw it out and print add
      // this contribution to the 'unresolved' errors to be printed at the end
      while (ir.size() > 0)
      {
        IntervalResult& rTest(ir.front());
        double lengthTest = rTest.delta;
        if (lengthTest < lengthMin)
        {
          warning("KC:Interval was divided too many times", iterationMax, rTest.a, rTest.b, rTest.err, ir.size());
          //    warning("KC:current result=",result,"error=",err);
          if (absErrorFlag)
            warning("Absolute accuracy = ", absError);
          if (relErrorFlag)
            warning("Relative accuracy = ", relError, "->absolute accuracy=", relError * std::abs(result));
          // this means there is a problem with the integrand->you could exit here
          //    exit(1);
          errorUnresolved += rTest.err;
          //    PrintList();
          ir.pop_front();
        }
        else
          break;
      }
      // do you want to exit with a warning after the first unresolved sub-interval occurred?
      if (ir.size() == 0 || errorUnresolved > 0.0)
        break;
      // or print as many as occur
      //      if (ir.size()==0) break;

      IntervalResult& r(ir.front());

      double center = 0.5 * (r.a + r.b);
      IntervalResult r1(r.a, center, 0.5 * r.delta);
      IntervalResult r2(center, r.b, 0.5 * r.delta);

      GK(r1);
      GK(r2);

      // must not use r after popping
      result += r1.result + r2.result - r.result;
      err    += r1.err + r2.err - r.err;

#ifdef PRINT_IT
      std::cout.setf(std::ios::scientific);
      std::cout << "Refined [ " << r.a << " " << r.b << " ] err/L=" << (r1.err + r2.err) / (r.b - r.a)
                << " error=" << err << std::endl;
#endif

      // must remove old element first because new ones could move to top
      ir.pop_front();

      Insert(r1);
      Insert(r2);

      // In rare events, the error decreases over may (>10) orders of magnitude
      // during the refinement. Rounding errors from the beginning can prevent
      // err from becoming small enough. Recompute err after a substantial decrease.
      if (err < 1e-6 * errLast)
      {
        err     = RecomputeError();
        errLast = err;
      }

      //      CheckError(err);
      //       PrintList();

      const bool relOk = (err < relError * std::abs(result) || result == 0.0);
      const bool absOk = (err < absError);

      if (absErrorFlag && relErrorFlag)
      {
        quitFlag = andFlag ? (relOk && absOk) : (relOk || absOk);
      }
      else
      {
        quitFlag = absErrorFlag ? absOk : relOk;
      }

    } while (!quitFlag);

#ifdef PRINT_IT
    PrintMax();
#endif

    if (errorUnresolved > 0.0)
    {
      warning("KC:Unresolved error sum=", errorUnresolved, "for integration interval", a, b);
      warning("KC:--> Result=", result, "total error=", err,
              "rel. error=", ((result != 0.0) ? err / std::abs(result) : 0.0));
      //      if (absErrorFlag) warning("Absolute accuracy = ",absError);
      //      if (relErrorFlag) warning("Relative accuracy = ",relError,
      //        "->absolute accuracy=",relError*std::abs(result));
    }

    //    CheckList();
#ifdef PRINT_IT
    std::cout << "End integration" << std::endl;
#endif
    double sum       = 0.0;
    int numIntervals = 0;
    for (typename std::list<IntervalResult>::iterator p = ir.begin(); p != ir.end(); p++)
    {
      sum += p->result;
      numIntervals++;
    }

    //     if (numIntervals > 2000)
    //       std::cerr << "Number of intervals = " << numIntervals << std::endl;

    double badSum = std::abs((result - sum) / sum);
    if ((badSum > 1.0e-7) && (std::abs(result - sum) > absError))
    {
      std::cerr << "absError tolerance = " << absError << std::endl;
      std::cerr << "Percent error = " << badSum * 100.0 << std::endl;
      std::cerr << "Number of intervals = " << numIntervals << std::endl;
      std::cerr << "result = " << result << " sum = " << sum << std::endl;
    }


    return result;
  }

Integrate() [2/3]

double Integrate ( const double  a, const double  b, const double  acc  )

inline

Definition at line 450 of file GKIntegration.h.

References GKIntegration< F, GKRule >::Integrate(), and GKIntegration< F, GKRule >::relativeErrors.

  {
    return Integrate(a, b, acc, !relativeErrors, acc, relativeErrors, false);
  }

Integrate() [3/3]

double Integrate ( const double  a, const double  b, const double  accAbs, const double  accRel, const bool  andFlag  )

inline

Definition at line 454 of file GKIntegration.h.

References GKIntegration< F, GKRule >::Integrate().

  {
    return Integrate(a, b, accAbs, true, accRel, true, andFlag);
  }

PrintList()

void PrintList ( )

inlineprivate

Definition at line 169 of file GKIntegration.h.

References BMWrite2, and GKIntegration< F, GKRule >::ir.

Referenced by GKIntegration< F, GKRule >::CheckList().

  {
    std::cout << "/------------------------------------------\\" << std::endl;
    int i = 0;
    for (typename std::list<IntervalResult>::iterator p = ir.begin(); p != ir.end(); p++)
    {
      BMWrite2(i, *p);
      i++;
    }
    std::cout << "\\------------------------------------------/" << std::endl;
  }

PrintMax()

void PrintMax ( )

inlineprivate

Definition at line 183 of file GKIntegration.h.

References BMWrite, and GKIntegration< F, GKRule >::ir.

Referenced by GKIntegration< F, GKRule >::Integrate().

  {
    typename std::list<IntervalResult>::iterator rMax(ir.begin());

    for (typename std::list<IntervalResult>::iterator r = ir.begin()++; r != ir.end(); r++)
    {
      if ((*r).ErrorL() > (*rMax).ErrorL())
        rMax = r;
    }

    BMWrite(*rMax);
  }

RecomputeError()

double RecomputeError ( )

inlineprivate

Definition at line 241 of file GKIntegration.h.

References GKIntegration< F, GKRule >::ir.

Referenced by GKIntegration< F, GKRule >::Integrate().

  {
    double errorSum = 0.0;

    if (ir.size() > 0)
    {
      for (typename std::list<IntervalResult>::iterator p = ir.begin(); p != ir.end(); p++)
      {
        errorSum += (*p).err;
      }
    }

    return errorSum;
  }

SetAbsoluteErrorMode()

void SetAbsoluteErrorMode ( )

inline

Definition at line 117 of file GKIntegration.h.

References GKIntegration< F, GKRule >::relativeErrors.

{ relativeErrors = false; }

SetRelativeErrorMode()

void SetRelativeErrorMode ( )

inline

Definition at line 118 of file GKIntegration.h.

References GKIntegration< F, GKRule >::relativeErrors.

{ relativeErrors = true; }

Member Data Documentation

f

F& f

private

Definition at line 111 of file GKIntegration.h.

Referenced by GKIntegration< F, GKRule >::GKGeneral().

ir

std::list<IntervalResult> ir

private

Definition at line 110 of file GKIntegration.h.

Referenced by GKIntegration< F, GKRule >::CheckError(), GKIntegration< F, GKRule >::CheckList(), GKIntegration< F, GKRule >::Insert(), GKIntegration< F, GKRule >::Integrate(), GKIntegration< F, GKRule >::PrintList(), GKIntegration< F, GKRule >::PrintMax(), and GKIntegration< F, GKRule >::RecomputeError().

relativeErrors

bool relativeErrors

private

Definition at line 112 of file GKIntegration.h.

Referenced by GKIntegration< F, GKRule >::Integrate(), GKIntegration< F, GKRule >::SetAbsoluteErrorMode(), and GKIntegration< F, GKRule >::SetRelativeErrorMode().

---

The documentation for this class was generated from the following file:

• /home/pk7/projects/qmc/for_cron_doxygen/qmcpack/src/QMCTools/ppconvert/src/common/GKIntegration.h