include/adapt02.hpp
00001 /* adapt02.hpp
00002  * 
00003  * Structures for the adaptive quadrature routines.
00004  */
00005 
00006 #ifndef MULTIQUAD_ADAPT02_HPP
00007 #define MULTIQUAD_ADAPT02_HPP
00008 
00016 template <typename FP>
00017 struct integral 
00018 {
00019     FP *a, *b, *result, *rawerr, *resabs;
00020     unsigned int dim;
00021 };
00022 
00026 template<typename FP>
00027 struct less_abs 
00028 {
00029     bool operator()(const integral<FP> e1, const integral<FP> e2)
00030     {   
00031         return *(e1.rawerr) < *(e2.rawerr);
00032     }
00033 };
00034 
00038 template<typename FP>
00039 struct less_rel 
00040 {
00041     bool operator()(const integral<FP> e1, const integral<FP> e2)
00042     {   
00043         return (*(e1.rawerr) / *(e1.resabs)) < (*(e2.rawerr) / *(e2.resabs));
00044     }
00045 };
00046 
00053 template <typename FP>
00054 class integralblock 
00055 {
00056     unsigned int dim_, sub1d_, N_;
00057     FP *a_, *b_;
00058     integral<FP>* integrals_;
00059     
00060 public:
00062     void init(unsigned int sub1d, unsigned int dim, const FP* a0, const FP* b0);
00063     unsigned int size() {return N_;}
00065     const FP* a() const {return a_;} 
00067     const FP* b() const {return b_;} 
00068     FP* result;
00069     FP* rawerr;
00070     FP* resabs;
00072     integral<FP>& region(unsigned int j)  {return integrals_[j];} 
00073 
00074     integralblock() : N_(0), a_(NULL), b_(NULL), result(NULL), rawerr(NULL), resabs(NULL) {}
00075     ~integralblock();
00076 };
00077 
00078 template <typename FP>
00079 void integralblock<FP>::init(unsigned int sub1d, unsigned int dim, const FP* a0, const FP* b0)
00080 {
00081     dim_ = dim; sub1d_ = sub1d; N_ = 1;
00082     for (int k = 0; k < dim_; ++k) 
00083         N_ *= sub1d_;
00084     
00085     a_ = new FP[N_ * dim_]; 
00086     b_ = new FP[N_ * dim_];
00087     result = new FP[N_]; 
00088     rawerr = new FP[N_]; 
00089     resabs = new FP[N_];
00090     integrals_ = new integral<FP>[N_];
00091     
00092     for (unsigned int j = 0; j < N_; ++j) {
00093         integrals_[j].a = &a_[j * dim_];
00094         integrals_[j].b = &b_[j * dim_];
00095         integrals_[j].result = &result[j];
00096         integrals_[j].rawerr = &rawerr[j];
00097         integrals_[j].resabs = &resabs[j];
00098 
00099         /* The main region is the product of intervals (a0[i], b0[i]).
00100          * The "multi-indices" k=(k0,...,k_{dim-1}) enumerate the subintervals 
00101          * for respective axes. The latter are used to locate the boundaries of 
00102          * each subregion.
00103          */
00104         unsigned int m = 1; 
00105         for (unsigned int i = 0; i < dim_; ++i) {
00106             unsigned int ki = (j / m) % sub1d_;
00107             FP len = (b0[i] - a0[i]) / sub1d_;
00108             // "left" endpoint of ki-th subinterval:
00109             a_[j * dim_ + i] = ki * len + a0[i]; 
00110             b_[j * dim_ + i] = (ki + 1) * len + a0[i];
00111             m *= sub1d_;
00112         }
00113         
00114         result[j] = FP(0);
00115         rawerr[j] = FP(0);
00116     }   
00117 }
00118 
00119 template <typename FP>
00120 integralblock<FP>::~integralblock()
00121 {
00122     if (N_) {
00123         delete[] a_; 
00124         delete[] b_; 
00125         delete[] result; 
00126         delete[] rawerr; 
00127         delete[] resabs;
00128         delete[] integrals_;
00129     }
00130 }
00131 
00148 #endif // MULTIQUAD_ADAPT02_HPP
 All Classes Functions Variables