1 /*
2   This file is part of MADNESS.
3 
4   Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6   This program is free software; you can redistribute it and/or modify
7   it under the terms of the GNU General Public License as published by
8   the Free Software Foundation; either version 2 of the License, or
9   (at your option) any later version.
10 
11   This program is distributed in the hope that it will be useful,
12   but WITHOUT ANY WARRANTY; without even the implied warranty of
13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14   GNU General Public License for more details.
15 
16   You should have received a copy of the GNU General Public License
17   along with this program; if not, write to the Free Software
18   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20   For more information please contact:
21 
22   Robert J. Harrison
23   Oak Ridge National Laboratory
24   One Bethel Valley Road
25   P.O. Box 2008, MS-6367
26 
27   email: harrisonrj@ornl.gov
28   tel:   865-241-3937
29   fax:   865-572-0680
30 
31   $Id$
32 */
33 #ifndef MADNESS_MRA_ADQUAD_H__INCLUDED
34 #define MADNESS_MRA_ADQUAD_H__INCLUDED
35 
36 #include <madness/mra/legendre.h>
37 
38 namespace madness {
39 
40     namespace detail {
41         template <typename T>
norm(const T & t)42         double norm(const T& t) {
43             return std::abs(t);
44         }
45 
46         template <typename T>
norm(const Tensor<T> & t)47         double norm(const Tensor<T>& t) {
48             return t.normf();
49         }
50     }
51 
52     template <typename funcT>
do_adq(double lo,double hi,const funcT & func,int n,const double * x,const double * w)53     typename funcT::returnT do_adq(double lo, double hi, const funcT& func,
54                                    int n, const double* x, const double* w) {
55         // x and w provide the Gauss-Legendre quadrature rule of order n on [0,1]
56         double range = (hi-lo);
57         typename funcT::returnT sum = func(lo + range*x[0])*w[0];
58         for (int i=1; i<n; ++i) sum += func(lo + range*x[i])*w[i];
59         return sum*range;
60     }
61 
62 
63     template <typename funcT>
adq1(double lo,double hi,const funcT & func,double thresh,int n,const double * x,const double * w,int level)64     typename funcT::returnT adq1(double lo, double hi, const funcT& func, double thresh,
65                                  int n, const double* x, const double* w, int level) {
66         static const int MAX_LEVEL=14;
67         double d = (hi-lo)/2;
68         // Twoscale by any other name would smell just as sweet.
69         typename funcT::returnT full = do_adq(lo, hi, func, n, x, w);
70         typename funcT::returnT half = do_adq(lo, lo+d, func, n, x, w) + do_adq(lo+d, hi, func, n, x, w);
71 
72         double abserr = madness::detail::norm(full-half);
73         double norm = madness::detail::norm(half);
74         double relerr = (norm==0.0) ? 0.0 : abserr/norm;
75         //for (int i=0; i<level; ++i) std::cout << "! ";
76         //std::cout << norm << " " << abserr << " " << relerr << " " << thresh << std::endl;
77 
78         bool converged = (relerr < 1e-14) || (abserr<thresh && relerr<0.01);
79         if (converged) {
80             return half;
81         }
82         else {
83             if (level == MAX_LEVEL) {
84                 //throw "Adaptive quadrature: failed : runaway refinement?";
85                 return half;
86             }
87             else {
88                 return adq1(lo, lo+d, func, thresh*0.5, n, x, w, level+1) +
89                        adq1(lo+d, hi, func, thresh*0.5, n, x, w, level+1);
90             }
91         }
92     }
93 
94     template <typename funcT>
adq(double lo,double hi,const funcT & func,double thresh)95     typename funcT::returnT adq(double lo, double hi, const funcT& func, double thresh) {
96         const int n = 20;
97         double x[n], y[n];
98         gauss_legendre(n, 0.0, 1.0, x, y);
99 
100         return adq1(lo, hi, func, thresh, n, x, y, 0);
101     }
102 
103     namespace detail {
104         struct adqtest {
105             typedef double returnT;
operatoradqtest106             double operator()(double x) const {
107                 // int(exp(-x^2)*cos(a*x),x=-inf..inf) = sqrt(pi)*exp(-a^2/4)
108                 return exp(-x*x)*cos(3*x);
109             }
110 
exactadqtest111             static double exact() {
112                 const double pi = 3.1415926535897932384;
113                 return sqrt(pi)*exp(-9.0/4.0);
114             }
115 
runtestadqtest116             static bool runtest() {
117                 double test = madness::adq(-6.0,6.0,adqtest(),1e-14);
118                 return std::abs(test-exact())<1e-14;
119             }
120 
121         };
122     }
123 }
124 
125 #endif // MADNESS_MRA_ADQUAD_H__INCLUDED
126