00001 # ifndef CPPAD_ROMBERG_MUL_INCLUDED
00002 # define CPPAD_ROMBERG_MUL_INCLUDED
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172 # include <cppad/romberg_one.hpp>
00173 # include <cppad/check_numeric_type.hpp>
00174 # include <cppad/check_simple_vector.hpp>
00175
00176 namespace CppAD {
00177
00178 template <class Fun, class FloatVector>
00179 class SliceLast {
00180 typedef typename FloatVector::value_type Float;
00181 private:
00182 Fun *F;
00183 size_t last;
00184 FloatVector x;
00185 public:
00186 SliceLast( Fun *F_, size_t last_, const FloatVector &x_ )
00187 : F(F_) , last(last_), x(last + 1)
00188 { size_t i;
00189 for(i = 0; i < last; i++)
00190 x[i] = x_[i];
00191 }
00192 double operator()(const Float &xlast)
00193 { x[last] = xlast;
00194 return (*F)(x);
00195 }
00196 };
00197
00198 template <class Fun, class SizeVector, class FloatVector, class Float>
00199 class IntegrateLast {
00200 private:
00201 Fun *F;
00202 const size_t last;
00203 const FloatVector a;
00204 const FloatVector b;
00205 const SizeVector n;
00206 const SizeVector p;
00207 Float esum;
00208 size_t ecount;
00209
00210 public:
00211 IntegrateLast(
00212 Fun *F_ ,
00213 size_t last_ ,
00214 const FloatVector &a_ ,
00215 const FloatVector &b_ ,
00216 const SizeVector &n_ ,
00217 const SizeVector &p_ )
00218 : F(F_) , last(last_), a(a_) , b(b_) , n(n_) , p(p_)
00219 { }
00220 Float operator()(const FloatVector &x)
00221 { Float r, e;
00222 SliceLast<Fun, FloatVector > S(F, last, x);
00223 r = CppAD::RombergOne(
00224 S, a[last], b[last], n[last], p[last], e
00225 );
00226 esum = esum + e;
00227 ecount++;
00228 return r;
00229 }
00230 void ClearEsum(void)
00231 { esum = 0.; }
00232 Float GetEsum(void)
00233 { return esum; }
00234
00235 void ClearEcount(void)
00236 { ecount = 0; }
00237 size_t GetEcount(void)
00238 { return ecount; }
00239 };
00240
00241 template <class Fun, class SizeVector, class FloatVector, size_t m>
00242 class RombergMul {
00243 typedef typename FloatVector::value_type Float;
00244 public:
00245 RombergMul(void)
00246 { }
00247 Float operator() (
00248 Fun &F ,
00249 const FloatVector &a ,
00250 const FloatVector &b ,
00251 const SizeVector &n ,
00252 const SizeVector &p ,
00253 Float &e )
00254 { Float r;
00255
00256 typedef IntegrateLast<
00257 Fun ,
00258 SizeVector ,
00259 FloatVector ,
00260 Float > IntegrateOne;
00261
00262 IntegrateOne Fm1(&F, m-1, a, b, n, p);
00263 RombergMul<
00264 IntegrateOne,
00265 SizeVector ,
00266 FloatVector ,
00267 m-1 > RombergMulM1;
00268
00269 Fm1.ClearEsum();
00270 Fm1.ClearEcount();
00271
00272 r = RombergMulM1(Fm1, a, b, n, p, e);
00273
00274 size_t i, j;
00275 Float prod = 1;
00276 size_t pow2 = 1;
00277 for(i = 0; i < m-1; i++)
00278 { prod *= (b[i] - a[i]);
00279 for(j = 0; j < (n[i] - 1); j++)
00280 pow2 *= 2;
00281 }
00282 assert( Fm1.GetEcount() == (pow2+1) );
00283
00284 e = e + Fm1.GetEsum() * prod / Fm1.GetEcount();
00285
00286 return r;
00287 }
00288 };
00289
00290 template <class Fun, class SizeVector, class FloatVector>
00291 class RombergMul <Fun, SizeVector, FloatVector, 1> {
00292 typedef typename FloatVector::value_type Float;
00293 public:
00294 Float operator() (
00295 Fun &F ,
00296 const FloatVector &a ,
00297 const FloatVector &b ,
00298 const SizeVector &n ,
00299 const SizeVector &p ,
00300 Float &e )
00301 { Float r;
00302 typedef IntegrateLast<
00303 Fun ,
00304 SizeVector ,
00305 FloatVector ,
00306 Float > IntegrateOne;
00307
00308
00309 CheckSimpleVector<Float, FloatVector>();
00310
00311
00312 CheckNumericType<Float>();
00313
00314 IntegrateOne F0(&F, 0, a, b, n, p);
00315
00316 F0.ClearEsum();
00317 F0.ClearEcount();
00318
00319 r = F0(a);
00320
00321 assert( F0.GetEcount() == 1 );
00322 e = F0.GetEsum();
00323
00324 return r;
00325 }
00326 };
00327
00328 }
00329
00330 # endif