Allow data_grid to use uBLAS vectors as DATA_TYPE

[campreilly]

There a cases where it can be useful to interpolate more than one ordinate using the same abscissa. However, There are several instances were we assume that the abscissa is the same date type as the ordinate. It may be much faster to interpolate these items in parallel than it is to create a separate data_grid for each ordinate.

This is pretty easy to do for the linear version of interpolation

Before the change:

        DATA_TYPE linear(int dim, const size_t* index, const double* location,
                DATA_TYPE& deriv, DATA_TYPE* deriv_vec) const
        {
            DATA_TYPE result, da, db;

            // build interpolation coefficients

            const DATA_TYPE a = interp(dim - 1, index, location, da, deriv_vec);
            size_t next[NUM_DIMS];
            memcpy(next, index, NUM_DIMS * sizeof(size_t));
            ++next[dim];
            const DATA_TYPE b = interp(dim - 1, next, location, db, deriv_vec);
            const size_t k = index[dim];
            seq_vector* ax = _axis[dim];

            // compute field value in this dimension

            const DATA_TYPE h = (DATA_TYPE) ax->increment(k) ;
            const DATA_TYPE u = (location[dim] - (*ax)(k)) / h ;
            result = a * (1.0 - u) + b * u;

            // compute derivative in this dimension and prior dimension

            if (deriv_vec) {
                deriv = (b - a) / h ;
                deriv_vec[dim] = deriv;
                if (dim > 0) {
                    deriv_vec[dim - 1] = da * (1.0 - u) + db * u;
                }
            }

            // use results for dim+1 iteration

            return result;
        }

After the change:

        DATA_TYPE linear(int dim, const size_t* index, const double* location,
                DATA_TYPE& deriv, DATA_TYPE* deriv_vec) const
        {
            DATA_TYPE result, da, db;

            // build interpolation coefficients

            const DATA_TYPE a = interp(dim - 1, index, location, da, deriv_vec);
            size_t next[NUM_DIMS];
            memcpy(next, index, NUM_DIMS * sizeof(size_t));
            ++next[dim];
            const DATA_TYPE b = interp(dim - 1, next, location, db, deriv_vec);
            const size_t k = index[dim];
            seq_vector* ax = _axis[dim];

            // compute field value in this dimension

            const double h = (double) ax->increment(k) ;
            const double u = (location[dim] - (*ax)(k)) / h ;
            result = a * (1.0 - u) + b * u;

            // compute derivative in this dimension and prior dimension

            if (deriv_vec) {
                deriv = (b - a) / h ;
                deriv_vec[dim] = deriv;
                if (dim > 0) {
                    deriv_vec[dim - 1] = da * (1.0 - u) + db * u;
                }
            }

            // use results for dim+1 iteration

            return result;
        }

But, is much harder for pchip, basically because of the cases where we check the sing of the slope. I'm also unclear on how it should effect the ability to write to a netCDF file.

So we'll put this off for another day.

[Tibonium]

Adding a note here only, as Sean and I had talked about this can be accomplished with functors. I have a few concepts that could address this concern, but will not begin development at this time.

[Tibonium]

I was thinking that I'd steered you in the wrong direction slighting. The change to this line:

initialize<DATA_TYPE>::zero(dy0, dy1, dy2, dy3) ;

should be

initialize<DATA_TYPE>::zero(dy0, dy1, dy2, dy3,*_data) ;

which will properly address the known issue with VS2010 and provide a sizing element required by the functor.

[campreilly]

Had some valgrind errors.

[campreilly]

Fixed valgrid errors.