From 954e605416cec5fbca92f74504ac5ca52c444789 Mon Sep 17 00:00:00 2001 From: jthorn Date: Thu, 28 Jun 2001 16:20:49 +0000 Subject: design changes to support Cactus-style optimizations, start implementation of interpolate() git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@117 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/jtutil/interpolate.cc | 78 +++++++++++++++++++++++++++++++++----- src/jtutil/interpolate.hh | 96 +++++++++++++++++++++++++++++++++-------------- 2 files changed, 136 insertions(+), 38 deletions(-) (limited to 'src/jtutil') diff --git a/src/jtutil/interpolate.cc b/src/jtutil/interpolate.cc index b1aef2f..661f81d 100644 --- a/src/jtutil/interpolate.cc +++ b/src/jtutil/interpolate.cc @@ -4,6 +4,7 @@ // interpolation::Lagrange_uniform_1d::Lagrange_uniform // interpolation::Lagrange_uniform_1d::setup_data_x // interpolation::Lagrange_uniform_1d::setup_data_y +// interpolation::Lagrange_uniform_1d::setup_interp_x // interpolation::Lagrange_uniform_1d::bdod // interpolation::Lagrange_uniform_1d::interpolate // interpolation::Lagrange_uniform_1d::Jacobian @@ -30,8 +31,26 @@ namespace interpolation template Lagrange_uniform_1d::Lagrange_uniform_1d(int order_in) : order(order_in), N_interp(order_in+1), + x_interp_valid(false), x_interp(0.0), end_action('?'), x_map_ptr(NULL), y_array(NULL) { +switch (order) + { +case 1: +case 2: +case 3: +case 4: +case 5: +case 6: + // yes, we support this ==> no-op here + break; +default: + error_exit(ERROR_EXIT, +"***** Lagrange_uniform_1d::Lagrange_uniform_1d():\n" +" order=%d isn't supported!\n", +, + order); /*NOTREACHED*/ + } } } @@ -47,6 +66,7 @@ template Lagrange_uniform_1d::setup_data_x(const linear_map& x_map_in) { x_map_ptr = &x_map_in; + if (x_map_ptr->N_points() < N_interp) then error_exit(ERROR_EXIT, "***** Lagrange_uniform_1d::setup_data_x():\n" @@ -83,18 +103,58 @@ y_array = y_array_in; //***************************************************************************** // -// This function does the actual interpolation. +// This function is used by client code to set up +// * the x coordinate at which interpolation is to be done +// * what sort of off-the-end-of-the-data handling is desired // namespace interpolation { template - fpt Lagrange_uniform_1d::interpolate(fpt x, char end_action = 'o') - const + void Lagrange_uniform_1d::setup_interp_x(fpt x_interp_in, + char end_action_in = 'o') +{ +x_interp = x_interp_in, +end_action = end_action_in; +x_interp_valid = true; +} + } + +//***************************************************************************** + +// +// This function is the user-level driver to do the actual interpolation. +// +namespace interpolation + { +template + fpt Lagrange_uniform_1d::interpolate() const { assert(x_map_ptr != NULL); assert(y_array != NULL); +assert(x_interp_valid); -// FIXME FIXME FIXME +int i_lo, i_hi; +const int i_center = bdod(i_lo, i_hi); + +const fpt x_relative = x - x_map_ptr->fp_of_int(i_center); +const fpt x_scaled = x_map_ptr->inverse_delta_fp() * x_relative; + +switch (order) + { +case 1: return interpolate_local_order1(x_scaled, & y_array[i_center]); +case 2: return interpolate_local_order2(x_scaled, & y_array[i_center]); +case 3: return interpolate_local_order3(x_scaled, & y_array[i_center]); +case 4: return interpolate_local_order4(x_scaled, & y_array[i_center]); +case 5: return interpolate_local_order5(x_scaled, & y_array[i_center]); +case 6: return interpolate_local_order6(x_scaled, & y_array[i_center]); +default: + error_exit(PANIC_EXIT, +"***** Lagrange_uniform_1d::intgerpolate():\n" +" order=%d isn't supported!\n", +" (this should never happen!)\n" +, + order); /*NOTREACHED*/ + } } } @@ -110,9 +170,7 @@ assert(y_array != NULL); namespace interpolation { template - void Lagrange_uniform_1d::bdod(fpt x, char end_action = 'o', - int& i_lo, int& i_hi) - const + int Lagrange_uniform_1d::bdod(int& i_lo, int& i_hi) const { // // Suppose we have a data array indexed by [i], from which we wish to @@ -140,6 +198,7 @@ template // assert(x_map_ptr != NULL); +assert(x_valid_flag); const enum linear_map::noninteger_action floor_or_round = jtutil::is_even(N_interp) ? linear_map::floor @@ -172,6 +231,7 @@ case 'e': double(x_map_ptr->max_fp()), i_center, i_lo, i_hi); /*NOTREACHED*/ case 'o': +case '*': if (i_lo < x_map_ptr->min_int()) then { const int shift = x_map_ptr->min_int() - i_lo; @@ -189,11 +249,11 @@ case 'o': break; } break; - -case '*': + default: error_exit(ERROR_EXIT, "***** Lagrange_uniform::bdod(): invalid end_action='%c'='0x%02x'!\n", +" (this should never happen!)\n" int(end_action), unsigned int(end_action)); /*NOTREACHED*/ } diff --git a/src/jtutil/interpolate.hh b/src/jtutil/interpolate.hh index c636fa5..715ac65 100644 --- a/src/jtutil/interpolate.hh +++ b/src/jtutil/interpolate.hh @@ -68,7 +68,7 @@ // interpolation position lies (fuzzily) outside the data range. // '*' ==> Like 'o', but also allow the interpolation position to lie // (fuzzily) outside the data range. In other words, this case -// case accepts *any* interpolation position. +// accepts *any* interpolation position. // // For Lagrange interpolation all three possibilities are valid; // for spline interpolation the "off-centering as needed" is intrinsic @@ -89,15 +89,11 @@ // 4. Set up the x coordinates where interpolated data is desired. // 5. Do the actual interpolation. // -// The routines here all allow 2-5, 3-5, and/or 4-5 to be repeated as needed. // The idea is that earlier steps can preallocate storage arrays and/or // precompute coefficients to be used by later stages, so repeating the // later stages may have lower overheads than redoing the entire process. -// -// One could also imagine doing these steps in the order 12435, i.e. -// interpolating many different sets of date to the same position(s). -// The Cactus interpolator interface supports this nicely, but alas -// the interface here doesn't support this at all. :( :( +// The interface here requires that 1 be done first, then allows 2-4 +// to be done in any order before 5. // // @@ -109,15 +105,17 @@ //****************************************************************************** +namespace interpolation + { + +//************************************** + // // This class provides 1D Lagrange global polynomial interpolation // for uniformly spaced data. // // ... const qualifiers refer to results of interpolate() // - -namespace interpolation - { template class Lagrange_uniform_1d { @@ -127,12 +125,21 @@ public: void setup_data_x(const linear_map& x_map_in); // this function sets up the y data values + // n.b. indexing off this pointer is [i], i.e. it need not be 0-origin void setup_data_y(const fpt y_array_in[]); + // this function sets up + // * the x coordinate where interpolation is to be done + // * what sort of off-the-end-of-the-data handling is desired + void setup_interp_x(fpt x_interp_in, char end_action_in = 'o'); + // the actual interpolation - // ... needs setup_data_x() and setup_data_y() - // to have already been called - fpt interpolate(fpt x, char end_action = 'o') const; + // ... requires that + // setup_data_x() + // setup_data_y() + // setup_interp_x() + // have already been called + fpt interpolate(); // this function computes the backwards domain of dependence // of the interpolation, i.e. the interval of [i] such that @@ -141,20 +148,24 @@ public: // for interpolate() // // ... returns domain-of-dependence interval in [i_lo,i_hi] - // ... returns nominal center point (i_center in source code) - // as function value - // ... needs setup_data_x() to have already been called, - // but does *not* need setup_data_y() to have already been called - int bdod(fpt x, char end_action = 'o', int& i_lo, int& i_hi) const; - - // Jacobian d interpolate(x,end_action) / d y[i] - // ... needs setup_data_x() to have already been called, - // but does *not* need setup_data_y() to have already been called - fpt Jacobian(fpt x, char end_action = 'o', int i) const; + // ... returns nominal center point of underlying local interpolant + // (i_center in source code) as function value + // ... requires that + // setup_data_x() + // setup_interp_x() + // have already been called + int bdod(int& i_lo, int& i_hi) const; + + // Jacobian d (value returned by interpolate()) / d y[i] + // ... requires that + // setup_data_x() + // setup_interp_x() + // have already been called + fpt Jacobian(int i) const; // constructor, destructor + // ... compiler-generated no-op destructor is ok Lagrange_uniform_1d(int order_in); - ~Lagrange_uniform_1d() { } private: // we forbid copying and passing by value @@ -166,13 +177,40 @@ private: private: const int order; - const int N_interp; // number of points in local - // interpolator, i.e. - // linear = 2 points, - // quadratic = 3 points, etc + const int N_interp; // number of points in local + // interpolator, i.e. + // linear = 2 points, + // quadratic = 3 points, etc + + bool x_interp_valid; // valid flag for x_interp, end_action + fpt x_interp; + char end_action; // n.b. non-const pointers to const data! const linear_map *x_map_ptr; // --> client-owned data const fpt *y_array; // --> client-owned data }; - }; + +//************************************** + +// +// ***** underlying local interpolants of varying orders ***** +// + +// ... x here is in units of the grid spacing, +// relative to the center point +// ... these make no checks for data being in range +template + fpt interpolate_local_order1(const fpt x, const fpt *y_center); +template + fpt interpolate_local_order2(const fpt x, const fpt *y_center); +template + fpt interpolate_local_order3(const fpt x, const fpt *y_center); +template + fpt interpolate_local_order4(const fpt x, const fpt *y_center); +fpt interpolate_local_order5(const fpt x, const fpt *y_center); +fpt interpolate_local_order6(const fpt x, const fpt *y_center); + +//************************************** + + }; // close namespace interpolation:: -- cgit v1.2.3