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