.. _program_listing_file_src_util_integrate.cpp:

Program Listing for File integrate.cpp
======================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_util_integrate.cpp>` (``src_util/integrate.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #include <chrono>
   #include "global_parameter.hpp"
   #include "integrate.hpp"
   #include "vector3D.hpp"
   #include "QCM.hpp"
   #include "qcm_ED.hpp"
   #include "cubature.h"
   #ifdef _OPENMP
   #include <omp.h>
   #endif
   
   //------------------------------------------------------------------------------
   // Context structs and cubature callbacks (local to this file)
   
   namespace {
   
   struct WkContext {
       function<void(Complex, vector3D<double>&, const int*, double[])>& f;
       double w_domain;
       double w_offset;  
       double iw_cutoff; 
       bool high;        
       bool verb;
       int& count;
   };
   
   int wk_cb(unsigned ndim, size_t npts, const double *x, void *fdata, unsigned fdim, double *fval)
   {
       auto& ctx = *static_cast<WkContext*>(fdata);
       const int nv = (int)fdim;
   
       if(!ctx.high) {
           for(size_t i = 0; i < npts; i++) {
               const double *xi = x + i*ndim;
               Complex w(0, xi[0]*ctx.w_domain + ctx.w_offset);
               vector3D<double> k(ndim>1 ? xi[1] : 0.0, ndim>2 ? xi[2] : 0.0, ndim>3 ? xi[3] : 0.0);
               ctx.f(w, k, &nv, &fval[fdim*i]);
           }
       } else {
           for(size_t i = 0; i < npts; i++) {
               const double *xi = x + i*ndim;
               double iw = xi[0]*ctx.w_domain;
               double *fv = &fval[fdim*i];
               if(iw < ctx.iw_cutoff) {
                   for(unsigned j = 0; j < fdim; ++j) fv[j] = 0.0;
               } else {
                   Complex w(0, 1.0/iw);
                   vector3D<double> k(ndim>1 ? xi[1] : 0.0, ndim>2 ? xi[2] : 0.0, ndim>3 ? xi[3] : 0.0);
                   ctx.f(w, k, &nv, fv);
                   double iw2 = iw*iw;
                   for(unsigned j = 0; j < fdim; ++j) fv[j] /= iw2;
               }
           }
       }
       ctx.count += (int)npts;
       if(ctx.count % 1000 == 0 && ctx.verb) cout << ctx.count << " points" << endl;
       return 0;
   }
   
   
   struct KContext {
       function<void(vector3D<double>&, const int*, double[])>& f;
   };
   
   int k_cb(unsigned ndim, size_t npts, const double *x, void *fdata, unsigned fdim, double *fval)
   {
       auto& ctx = *static_cast<KContext*>(fdata);
       const int nv = (int)fdim;
   
       for(size_t i = 0; i < npts; i++) {
           const double *xi = x + i*ndim;
           vector3D<double> k(xi[0], ndim>1 ? xi[1] : 0.0, ndim>2 ? xi[2] : 0.0);
           ctx.f(k, &nv, &fval[fdim*i]);
       }
       return 0;
   }
   
   } // anonymous namespace
   
   
   
   
   void QCM::wk_integral(int dim, function<void (Complex w, vector3D<double> &k, const int *nv, double I[])> f, vector<double> &Iv, const double accuracy, bool verb)
   {
       int ndim = dim+1;
       int maxpoints;
   
       auto cubature = global_bool("use_pcubature") ? pcubature_v : hcubature_v;
       if(verb) cout << "Cubature integration (" << (global_bool("use_pcubature") ? "p" : "h") << "-adaptive)" << endl;
   
       int ncomp = (int)Iv.size();
   
       if(ndim==2)      maxpoints = 500000;
       else if(ndim==3) maxpoints = 10000000;
       else if(ndim==4) maxpoints = 20000000;
       else             maxpoints = 10000;
   
       int fail;
       const double xmin[] = {0,0,0,0};
       const double xmax[] = {1,1,1,1};
   
       const double small_scale = global_double("small_scale");
       const double large_scale = global_double("large_scale");
       const double iw_cutoff = 1.0/global_double("cutoff_scale");
       int count = 0;
   
       //------------------------------------------------------------------------------
       // first region : frequencies below small_scale
       {
           double w_domain = small_scale;
           vector<double> value(ncomp, 0.0), err(ncomp, 0.0);
           WkContext ctx{f, w_domain, 0.0, iw_cutoff, false, verb, count};
           auto t1 = std::chrono::high_resolution_clock::now();
           fail = cubature(ncomp, wk_cb, &ctx, ndim, xmin, xmax, (size_t)maxpoints, accuracy*M_PI/w_domain, 1e-10, ERROR_INDIVIDUAL, value.data(), err.data());
           if(fail) qcm_throw("error in Cubature integral : fail = "+to_string<int>(fail));
           if(verb) cout << "region 1 : " << (double)std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now()-t1).count()/1000 << " seconds" << endl;
           double fac = w_domain*M_1_PI;
           for(int i=0; i<ncomp; ++i) Iv[i] += value[i]*fac;
       }
   
       //------------------------------------------------------------------------------
       // second region : frequencies between small_scale and large_scale
       {
           double w_domain = large_scale - small_scale;
           vector<double> value(ncomp, 0.0), err(ncomp, 0.0);
           WkContext ctx{f, w_domain, small_scale, iw_cutoff, false, verb, count};
           auto t1 = std::chrono::high_resolution_clock::now();
           fail = cubature(ncomp, wk_cb, &ctx, ndim, xmin, xmax, (size_t)maxpoints, accuracy*M_PI/w_domain, 1e-10, ERROR_INDIVIDUAL, value.data(), err.data());
           if(fail) qcm_throw("error in Cubature integral : fail = "+to_string<int>(fail));
           if(verb) cout << "region 2 : " << (double)std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now()-t1).count()/1000 << " seconds" << endl;
           double fac = w_domain*M_1_PI;
           for(int i=0; i<ncomp; ++i) Iv[i] += value[i]*fac;
       }
   
       //------------------------------------------------------------------------------
       // third region : inverse frequencies below 1/large_scale
       {
           double w_domain = 1.0/large_scale;
           vector<double> value(ncomp, 0.0), err(ncomp, 0.0);
           WkContext ctx{f, w_domain, 0.0, iw_cutoff, true, verb, count};
           auto t1 = std::chrono::high_resolution_clock::now();
           fail = cubature(ncomp, wk_cb, &ctx, ndim, xmin, xmax, (size_t)maxpoints, accuracy*M_PI/w_domain, 1e-10, ERROR_INDIVIDUAL, value.data(), err.data());
           if(fail) qcm_throw("error in Cubature integral : fail = "+to_string<int>(fail));
           if(verb) cout << "region 3 : " << (double)std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now()-t1).count()/1000 << " seconds" << endl;
           double fac = w_domain*M_1_PI;
           for(int i=0; i<ncomp; ++i) Iv[i] += value[i]*fac;
       }
   }
   
   
   
   
   void QCM::k_integral_grid(int dim, int nkx, int nky, int nkz, function<void (vector3D<double> &k, const int *nv, double I[])> f, vector<double> &Iv)
   {
       int nv = (int)Iv.size();
       vector<double> I(nv, 0.0);
       double ikx = 1.0/nkx;
       double iky = 1.0/nky;
       double ikz = 1.0/nkz;
       if(dim==1){
           #pragma omp parallel
           {
               vector<double> Ik(nv);
               vector<double> Ip(nv, 0.0);
               #pragma omp for nowait
               for(int i = 0; i < nkx; i++){
                   vector3D<double> k(i*ikx, 0.0, 0.0);
                   f(k, &nv, Ik.data());
                   Ip += Ik;
               }
               #pragma omp critical
               I += Ip;
           }
           I *= ikx;
       }
       else if (dim==2){
           #pragma omp parallel
           {
               vector<double> Ik(nv);
               vector<double> Ip(nv, 0.0);
               #pragma omp for collapse(2) nowait
               for(int i = 0; i < nkx; i++){
                   for(int j = 0; j < nky; j++){
                       vector3D<double> k(i*ikx, j*iky, 0.0);
                       f(k, &nv, Ik.data());
                       Ip += Ik;
                   }
               }
               #pragma omp critical
               I += Ip;
           }
           I *= ikx*iky;
       }
       else if (dim==3){
           #pragma omp parallel
           {
               vector<double> Ik(nv);
               vector<double> Ip(nv, 0.0);
               #pragma omp for collapse(3) nowait
               for(int i = 0; i < nkx; i++){
                   for(int j = 0; j < nky; j++){
                       for(int l = 0; l < nkz; l++){
                           vector3D<double> k(i*ikx, j*iky, l*ikz);
                           f(k, &nv, Ik.data());
                           Ip += Ik;
                       }
                   }
               }
               #pragma omp critical
               I += Ip;
           }
           I *= ikx*iky*ikz;
       }
       Iv += I;
   }
   
   
   
   
   void QCM::k_integral(int dim, function<void (vector3D<double> &k, const int *nv, double I[])> f, vector<double> &Iv, const double accuracy, bool verb)
   {
       int maxpoints;
       if(dim==1)       maxpoints = 100000;
       else if(dim==2)  maxpoints = 300000;
       else             maxpoints = 4000000;
   
       int ncomp = (int)Iv.size();
       vector<double> value(ncomp, 0.0), err(ncomp, 0.0);
       int fail;
   
       auto cubature = global_bool("use_pcubature") ? pcubature_v : hcubature_v;
       const double xmin[] = {0,0,0,0};
       const double xmax[] = {1,1,1,1};
   
       KContext ctx{f};
       auto t1 = std::chrono::high_resolution_clock::now();
       fail = cubature(ncomp, k_cb, &ctx, dim, xmin, xmax, (size_t)maxpoints, accuracy, 1e-10, ERROR_INDIVIDUAL, value.data(), err.data());
       if(verb) cout << "Cubature : done in " << (double)std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now()-t1).count()/1000 << " seconds" << endl;
       for(int i=0; i<ncomp; ++i) Iv[i] += value[i];
   }