diff --git a/SparseGrids/tsgGridSequence.cpp b/SparseGrids/tsgGridSequence.cpp index 12e9ca783..5ebee339c 100644 --- a/SparseGrids/tsgGridSequence.cpp +++ b/SparseGrids/tsgGridSequence.cpp @@ -351,9 +351,10 @@ std::vector GridSequence::getMultiIndex(const double x[]){ return p; } void GridSequence::loadConstructedPoint(const double x[], const std::vector &y){ - auto p = getMultiIndex(x); + std::vector p = getMultiIndex(x); + std::vector scratch(num_dimensions); - if (MultiIndexManipulations::isLowerComplete(p, points)){ + if (MultiIndexManipulations::isLowerComplete(p, points, scratch)){ std::vector approx_value(num_outputs), surplus(num_outputs);; if (!points.empty()){ evaluate(x, approx_value.data()); diff --git a/SparseGrids/tsgIndexManipulator.hpp b/SparseGrids/tsgIndexManipulator.hpp index f5d824138..c09527346 100644 --- a/SparseGrids/tsgIndexManipulator.hpp +++ b/SparseGrids/tsgIndexManipulator.hpp @@ -440,14 +440,19 @@ MultiIndexSet createPolynomialSpace(const MultiIndexSet &tensors, std::function< * \ingroup TasmanianMultiIndexManipulations * \brief Assuming that \b mset is lower complete, return \b true if adding the \b point will preserve completeness. * + * \param point is the new point to add to an existing set + * \param mset is an existing lower complete set + * \param scratch is scratch space, must have size equal to point.size() and will be used for temporary storage + * the scratch reduces allocations + * * \endinternal */ -inline bool isLowerComplete(std::vector const &point, MultiIndexSet const &mset){ - auto dad = point; - for(auto &d : dad){ +inline bool isLowerComplete(std::vector const &point, MultiIndexSet const &mset, std::vector &scratch){ + std::copy(point.begin(), point.end(), scratch.begin()); + for(int &d : scratch){ if (d > 0){ d--; - if (mset.missing(dad)) return false; + if (mset.missing(scratch)) return false; d++; } } @@ -473,6 +478,7 @@ inline MultiIndexSet getLargestCompletion(MultiIndexSet const ¤t, MultiInd } } + std::vector scratch(num_dimensions); bool loopon = true; while(loopon){ Data2D update(num_dimensions, 0); @@ -481,9 +487,9 @@ inline MultiIndexSet getLargestCompletion(MultiIndexSet const ¤t, MultiInd if (!result.empty()) total += result; for(int i=0; i kid(total.getIndex(i), total.getIndex(i) + num_dimensions); - for(auto &k : kid){ + for(int &k : kid){ k++; // construct the kid in the new direction - if (!candidates.missing(kid) && result.missing(kid) && isLowerComplete(kid, total)) + if (!candidates.missing(kid) && result.missing(kid) && isLowerComplete(kid, total, scratch)) update.appendStrip(kid); k--; } @@ -508,6 +514,7 @@ template MultiIndexSet addExclusiveChildren(const MultiIndexSet &tensors, const MultiIndexSet &exclude, const std::vector level_limits){ int num_dimensions = (int) tensors.getNumDimensions(); Data2D tens(num_dimensions, 0); + std::vector scratch(num_dimensions); for(int i=0; i kid(t, t + num_dimensions); @@ -515,7 +522,7 @@ MultiIndexSet addExclusiveChildren(const MultiIndexSet &tensors, const MultiInde for(auto &k : kid){ k++; if (exclude.missing(kid) && tensors.missing(kid)){ // if the kid is not to be excluded and if not included in the current set - if (isLowerComplete(kid, tensors)){ + if (isLowerComplete(kid, tensors, scratch)){ if (limited){ if ((*ilimit == -1) || (k <= *ilimit)) tens.appendStrip(kid);