Skip to content

GitLab

Commit 7fb217ca authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

kdtreelib: added node pool.

parent 8ab4d09a
Hide whitespace changes
Inline Side-by-side
Showing with 201 additions and 96 deletions
src/Afterburner.cc
View file @ 7fb217ca
......@@ -2310,7 +2310,7 @@ int afterburner(int argc, char *argv[])
omp_set_num_threads(numThreads);
if ( omp_get_max_threads() > omp_get_num_procs() )
fprintf(stdout, " Number of threads is greater than number of Cores=%d!\n", omp_get_num_procs());
fprintf(stdout, " OpenMP: num_procs = %d max_threads = %d\n", omp_get_num_procs(), omp_get_max_threads());
fprintf(stdout, " OpenMP: num_procs=%d max_threads=%d\n", omp_get_num_procs(), omp_get_max_threads());
#else
if ( numThreads > 0 )
{
......
src/cdo.cc
View file @ 7fb217ca
......@@ -1557,10 +1557,10 @@ int main(int argc, char *argv[])
fprintf(stderr, "Warning: omp_get_max_threads() returns %d!\n", ompNumThreads);
if ( cdoVerbose )
{
fprintf(stderr, " OpenMP: num_procs = %d max_threads = %d", omp_get_num_procs(), omp_get_max_threads());
fprintf(stderr, " OpenMP: num_procs=%d max_threads=%d", omp_get_num_procs(), omp_get_max_threads());
#if defined(HAVE_OPENMP4)
#if !defined(__ICC)
fprintf(stderr, " num_devices = %d", omp_get_num_devices());
fprintf(stderr, " num_devices=%d", omp_get_num_devices());
#endif
#endif
fprintf(stderr, "\n");
......
src/grid_search.cc
View file @ 7fb217ca
......@@ -135,8 +135,8 @@ struct gridsearch *gridsearch_create_reg2d(bool is_cyclic, size_t dims[2], const
return gs;
}
struct kdNode *gs_create_kdtree(size_t n, const double *restrict lons, const double *restrict lats)
static
kdTree_t *gs_create_kdtree(size_t n, const double *restrict lons, const double *restrict lats)
{
struct kd_point *pointlist = (struct kd_point *) Malloc(n*sizeof(struct kd_point));
// see example_cartesian.c
......@@ -152,7 +152,7 @@ struct kdNode *gs_create_kdtree(size_t n, const double *restrict lons, const dou
{
point = pointlist[i].point;
LLtoXYZ_kd(lons[i], lats[i], point);
for ( size_t j = 0; j < 3; ++j )
for ( unsigned j = 0; j < 3; ++j )
{
min[j] = point[j] < min[j] ? point[j] : min[j];
max[j] = point[j] > max[j] ? point[j] : max[j];
......@@ -160,14 +160,15 @@ struct kdNode *gs_create_kdtree(size_t n, const double *restrict lons, const dou
pointlist[i].index = i;
}
struct kdNode *kdt = kd_buildTree(pointlist, n, min, max, 3, ompNumThreads);
kdTree_t *kdt = kd_buildTree(pointlist, n, min, max, 3, ompNumThreads);
if ( pointlist ) Free(pointlist);
if ( kdt == NULL ) cdoAbort("kd_buildTree failed!");
return kdt;
}
struct kdNode *gs_create_kdsph(size_t n, const double *restrict lons, const double *restrict lats)
static
kdTree_t *gs_create_kdsph(size_t n, const double *restrict lons, const double *restrict lats)
{
struct kd_point *pointlist = (struct kd_point *) Malloc(n*sizeof(struct kd_point)); // kd_point contains 3d point
// see example_cartesian.c
......@@ -185,7 +186,7 @@ struct kdNode *gs_create_kdsph(size_t n, const double *restrict lons, const doub
point[0] = lons[i];
point[1] = lats[i];
point[2] = 0; // dummy
for ( size_t j = 0; j < 2; ++j )
for ( unsigned j = 0; j < 2; ++j )
{
min[j] = point[j] < min[j] ? point[j] : min[j];
max[j] = point[j] > max[j] ? point[j] : max[j];
......@@ -193,7 +194,7 @@ struct kdNode *gs_create_kdsph(size_t n, const double *restrict lons, const doub
pointlist[i].index = i;
}
struct kdNode *kdt = kd_sph_buildTree(pointlist, n, min, max, ompNumThreads);
kdTree_t *kdt = kd_sph_buildTree(pointlist, n, min, max, ompNumThreads);
if ( pointlist ) Free(pointlist);
return kdt;
......@@ -443,7 +444,7 @@ struct gridsearch *gridsearch_create(size_t n, const double *restrict lons, cons
gs->n = n;
if ( n == 0 ) return gs;
gs->kdt = gs_create_kdtree(n, lons, lats);
gs->kdt = gs_create_kdtree(n, lons, lats);
gs->search_radius = cdo_default_search_radius();
......@@ -511,7 +512,7 @@ double gs_set_range(double *prange)
}
static
size_t gs_nearest_kdtree(kdNode *kdt, double lon, double lat, double *prange)
size_t gs_nearest_kdtree(kdTree_t *kdt, double lon, double lat, double *prange)
{
size_t index = GS_NOT_FOUND;
if ( kdt == NULL ) return index;
......@@ -522,7 +523,7 @@ size_t gs_nearest_kdtree(kdNode *kdt, double lon, double lat, double *prange)
kdata_t point[3];
LLtoXYZ_kd(lon, lat, point);
kdNode *node = kd_nearest(kdt, point, &range, 3);
kdNode *node = kd_nearest(kdt->node, point, &range, 3);
float frange = KDATA_INVSCALE(range);
if ( !(frange < range0) ) node = NULL;
......@@ -534,7 +535,7 @@ size_t gs_nearest_kdtree(kdNode *kdt, double lon, double lat, double *prange)
}
static
size_t gs_nearest_kdsph(kdNode *kdt, double lon, double lat, double *prange)
size_t gs_nearest_kdsph(kdTree_t *kdt, double lon, double lat, double *prange)
{
size_t index = GS_NOT_FOUND;
if ( kdt == NULL ) return index;
......@@ -546,7 +547,7 @@ size_t gs_nearest_kdsph(kdNode *kdt, double lon, double lat, double *prange)
point[0] = lon;
point[1] = lat;
kdNode *node = kd_nearest(kdt, point, &range, 3);
kdNode *node = kd_nearest(kdt->node, point, &range, 3);
float frange = KDATA_INVSCALE(range);
if ( !(frange < range0) ) node = NULL;
......@@ -667,7 +668,7 @@ struct pqueue *gridsearch_qnearest(struct gridsearch *gs, double lon, double lat
if ( gs )
{
result = kd_qnearest(gs->kdt, point, &range, nnn, 3);
result = kd_qnearest(gs->kdt->node, point, &range, nnn, 3);
// printf("range %g %g %g %p\n", lon, lat, range, node);
float frange = KDATA_INVSCALE(range);
......
src/grid_search.h
View file @ 7fb217ca
......@@ -34,7 +34,7 @@ struct gridsearch {
size_t dims[2];
struct gsNear *near;
struct kdNode *kdt;
struct kdTree *kdt;
struct gsFull *full;
double search_radius;
......
src/kdtreelib/kdtree.h
View file @ 7fb217ca
......@@ -12,6 +12,7 @@
_compPoints: compare index if points[axis] are equal
20171102: renamed kd_buildArg() to kd_initArg(), changed interface and memory handling
changed data pointer to size_t
20171112: added node pool
*/
#ifndef KDTREE_H_
#define KDTREE_H_
......@@ -22,6 +23,9 @@
#include <stdint.h>
#include <stdio.h>
#include <string.h>
#include <atomic>
//#define MEMPOOL 1
#define KD_FLOAT 1
#define KD_INT 2
......@@ -81,13 +85,31 @@ int qcmp(struct kd_point *a, struct kd_point *b, int axis)
typedef struct kdNode {
struct kdNode *left; /*!<the left child of the tree node */
struct kdNode *right; /*!<the right child of the tree node */
kdata_t location[KD_MAX_DIM]; /*!<vector to the node's location */
kdata_t min[KD_MAX_DIM]; /*!<vector to the min coordinates of the hyperrectangle */
kdata_t max[KD_MAX_DIM]; /*!<vector to the max coordinates of the hyperrectangle */
kdata_t location[KD_MAX_DIM]; /*!<vector to the node's location */
int split; /*!<axis along which the tree bifurcates */
size_t index; /*!<optional index value */
} kdNode;
typedef struct kdNodePool {
size_t size;
std::atomic<size_t> imp;
kdNode *pool;
} kdNodePool_t;
typedef struct kdTree {
kdNode *node;
kdNodePool_t *nodepool;
} kdTree_t;
kdNodePool_t *kd_nodepool_new(size_t n);
void kd_nodepool_delete(kdNodePool_t *pool);
/*!
* \struct resItem
* \brief result items, member of a priority queue
......@@ -113,13 +135,14 @@ typedef struct pqueue {
* \brief arguments passed to the threaded kd-Tree construction
*/
typedef struct kd_thread_data {
struct kd_point *points;
kdata_t min[KD_MAX_DIM];
kdata_t max[KD_MAX_DIM];
size_t nPoints;
int max_threads;
int depth;
int dim;
kdNodePool_t *nodepool;
struct kd_point *points;
kdata_t min[KD_MAX_DIM];
kdata_t max[KD_MAX_DIM];
size_t nPoints;
int max_threads;
int depth;
int dim;
} kd_thread_data;
......@@ -153,16 +176,19 @@ void kd_printNode(struct kdNode *node);
void kd_printTree(struct kdNode *node);
/* Functions for building and destroying trees */
struct kdNode *kd_allocNodeP(kdNodePool *nodepool, struct kd_point *points, size_t pivot,
kdata_t *min, kdata_t *max, int dim, int axis);
struct kdNode *kd_allocNode(struct kd_point *points, size_t pivot,
kdata_t *min, kdata_t *max, int dim, int axis);
void kd_destroyTree(struct kdNode *node);
void kd_initArg(struct kd_thread_data *d, struct kd_point *points, size_t nPoints,
void kd_doDestroyTree(struct kdNode *node);
void kd_destroyTree(kdTree_t *tree);
void kd_initArg(struct kd_thread_data *d, kdNodePool_t *nodepool, struct kd_point *points, size_t nPoints,
kdata_t *min, kdata_t *max, int depth, int max_threads, int dim);
struct kdNode *kd_buildTree(struct kd_point *points, size_t nPoints,
kdTree_t *kd_buildTree(struct kd_point *points, size_t nPoints,
kdata_t *min, kdata_t *max, int dim, int max_threads);
void *kd_doBuildTree(void *threadarg);
struct kdNode *kd_sph_buildTree(struct kd_point *points, size_t nPoints,
kdata_t *min, kdata_t *max, int max_threads);
kdTree_t *kd_sph_buildTree(struct kd_point *points, size_t nPoints,
kdata_t *min, kdata_t *max, int max_threads);
void *kd_sph_doBuildTree(void *threadarg);
/* Functions for range searches
......@@ -208,4 +234,4 @@ int kd_sph_doRange(struct kdNode *node, kdata_t *p, kdata_t *max_dist_sq,
int kd_insertResTree(struct kdNode *node, struct pqueue *res);
#endif /* KDTREE_H_ */
#endif /* KDTREE_H_ */
src/kdtreelib/kdtree_cartesian.cc
View file @ 7fb217ca
......@@ -69,13 +69,23 @@ kdata_t kd_min(kdata_t x, kdata_t y)
*
* \return root node of the tree
*/
struct kdNode *
kdTree_t *
kd_buildTree(struct kd_point *points, size_t nPoints,
kdata_t *min, kdata_t *max, int dim, int max_threads)
{
kdTree_t *tree = (kdTree_t *) calloc(1, sizeof(kdTree_t));
#ifdef MEMPOOL
tree->nodepool = kd_nodepool_new(nPoints*2);
#endif
struct kd_thread_data my_data;
kd_initArg(&my_data, points, nPoints, min, max, 0, max_threads, dim);
struct kdNode *tree = (kdNode *)kd_doBuildTree(&my_data);
kd_initArg(&my_data, tree->nodepool, points, nPoints, min, max, 0, max_threads, dim);
tree->node = (kdNode *)kd_doBuildTree(&my_data);
if ( tree->node == NULL )
{
free(tree);
tree = NULL;
}
return tree;
}
......
src/kdtreelib/kdtree_common.cc
View file @ 7fb217ca
......@@ -79,9 +79,10 @@ kd_printTree(struct kdNode *node)
/* End helper functions */
void kd_initArg(struct kd_thread_data *d, struct kd_point *points, size_t nPoints,
void kd_initArg(struct kd_thread_data *d, kdNodePool *nodepool, struct kd_point *points, size_t nPoints,
kdata_t *min, kdata_t *max, int depth, int max_threads, int dim)
{
d->nodepool = nodepool;
d->points = points;
d->nPoints = nPoints;
memcpy(d->min, min, dim*sizeof(kdata_t));
......@@ -97,6 +98,111 @@ void kd_initArg(struct kd_thread_data *d, struct kd_point *points, size_t nPoint
******************************************************************** */
/*!
* \brief free the kd-tree data structure,
*
* \param node the root node of the tree to be destroyed
*
* \param *destr a pointer to the destructor function for the data container.
*
* \return This function does not return a value
*/
static
void kd_freeNode(kdNode *node)
{
if ( node ) free(node);
}
kdNodePool_t *kd_nodepool_new(size_t n)
{
kdNodePool_t *nodepool = (kdNodePool_t *) malloc(sizeof(kdNodePool_t));
nodepool->size = n;
nodepool->imp = 0;
nodepool->pool = (kdNode *) malloc(n*sizeof(kdNode));
printf("kd_nodepool_new: size=%zu\n", n);
return nodepool;
}
void kd_nodepool_delete(kdNodePool_t *nodepool)
{
if ( nodepool == NULL ) return;
if ( nodepool->pool ) free(nodepool->pool);
free(nodepool);
}
struct kdNode *
kd_allocNode(struct kd_point *points, size_t pivot,
kdata_t *min, kdata_t *max, int axis, int dim)
{
struct kdNode *node;
if ((node = (kdNode *)kd_malloc(sizeof(kdNode), "kd_allocNode (node): ")) == NULL)
return NULL;
node->split = axis;
memcpy(node->location, points[pivot].point, dim * sizeof(kdata_t));
memcpy(node->min, min, dim * sizeof(kdata_t));
memcpy(node->max, max, dim * sizeof(kdata_t));
node->left = node->right = NULL;
node->index = 0;
return node;
}
struct kdNode *
kd_allocNodeP(kdNodePool *nodepool, struct kd_point *points, size_t pivot,
kdata_t *min, kdata_t *max, int axis, int dim)
{
struct kdNode *node;
if ( nodepool )
{
if ( nodepool->imp >= nodepool->size ) return NULL;
node = &nodepool->pool[nodepool->imp++];
}
else
{
if ((node = (kdNode *)kd_malloc(sizeof(kdNode), "kd_allocNode (node): ")) == NULL)
return NULL;
}
node->split = axis;
memcpy(node->location, points[pivot].point, dim * sizeof(kdata_t));
memcpy(node->min, min, dim * sizeof(kdata_t));
memcpy(node->max, max, dim * sizeof(kdata_t));
node->left = node->right = NULL;
node->index = 0;
return node;
}
void
kd_doDestroyTree(struct kdNode *node)
{
if ( node == NULL ) return;
kd_doDestroyTree(node->left);
kd_doDestroyTree(node->right);
kd_freeNode(node);
}
void
kd_destroyTree(kdTree_t *tree)
{
if ( tree == NULL ) return;
if ( tree->nodepool )
kd_nodepool_delete(tree->nodepool);
else
kd_doDestroyTree(tree->node);
free(tree);
}
void *kd_doBuildTree(void *threadarg)
{
kdata_t tmpMaxLeft[KD_MAX_DIM], tmpMinRight[KD_MAX_DIM];
......@@ -106,6 +212,7 @@ void *kd_doBuildTree(void *threadarg)
struct kd_thread_data argleft, argright;
struct kd_thread_data *my_data = (struct kd_thread_data *) threadarg;
kdNodePool_t *nodepool = my_data->nodepool;
struct kd_point *points = my_data->points;
size_t nPoints = my_data->nPoints;
kdata_t *min = my_data->min;
......@@ -120,7 +227,7 @@ void *kd_doBuildTree(void *threadarg)
int sortaxis = depth % dim;
if (nPoints == 1) {
if ((node = kd_allocNode(points, 0, min, max, sortaxis, dim)) == NULL)
if ((node = kd_allocNodeP(nodepool, points, 0, min, max, sortaxis, dim)) == NULL)
return NULL;
node->index = points[0].index;
return node;
......@@ -133,12 +240,12 @@ void *kd_doBuildTree(void *threadarg)
pmergesort(points, nPoints, sortaxis, max_threads);
size_t pivot = nPoints / 2;
if ((node = kd_allocNode(points, pivot, min, max, sortaxis, dim)) == NULL)
if ((node = kd_allocNodeP(nodepool, points, pivot, min, max, sortaxis, dim)) == NULL)
return NULL;
memcpy(tmpMaxLeft, max, dim * sizeof(kdata_t));
tmpMaxLeft[sortaxis] = node->location[sortaxis];
kd_initArg(&argleft, points, pivot, min, tmpMaxLeft,
kd_initArg(&argleft, nodepool, points, pivot, min, tmpMaxLeft,
depth + 1, max_threads / 2, dim);
if (max_threads > 1) {
......@@ -146,14 +253,14 @@ void *kd_doBuildTree(void *threadarg)
} else {
node->left = (kdNode *)kd_doBuildTree((void *) &argleft);
if (!node->left) {
kd_destroyTree(node);
kd_doDestroyTree(node);
return NULL;
}
}
memcpy(tmpMinRight, min, dim * sizeof(kdata_t));
tmpMinRight[sortaxis] = node->location[sortaxis];
kd_initArg(&argright, &points[pivot], nPoints - pivot, tmpMinRight, max,
kd_initArg(&argright, nodepool, &points[pivot], nPoints - pivot, tmpMinRight, max,
depth + 1, max_threads / 2, dim);
if (max_threads > 1) {
......@@ -161,7 +268,7 @@ void *kd_doBuildTree(void *threadarg)
} else {
node->right = (kdNode *)kd_doBuildTree((void *) &argright);
if (!node->right) {
kd_destroyTree(node);
kd_doDestroyTree(node);
return NULL;
}
}
......@@ -170,7 +277,7 @@ void *kd_doBuildTree(void *threadarg)
pthread_join(threads[0], (void **) (&node->left));
pthread_join(threads[1], (void **) (&node->right));
if (!node->left || !node->right) {
kd_destroyTree(node);
kd_doDestroyTree(node);
return NULL;
}
}
......@@ -178,51 +285,6 @@ void *kd_doBuildTree(void *threadarg)
return (void *) node;
}
static
void kd_freeNode(kdNode *node)
{
if ( node ) free(node);
}
struct kdNode *
kd_allocNode(struct kd_point *points, size_t pivot,
kdata_t *min, kdata_t *max, int axis, int dim)
{
struct kdNode *node;
if ((node = (kdNode *)kd_malloc(sizeof(kdNode), "kd_allocNode (node): ")) == NULL)
return NULL;
node->split = axis;
memcpy(node->location, points[pivot].point, dim * sizeof(kdata_t));
memcpy(node->min, min, dim * sizeof(kdata_t));
memcpy(node->max, max, dim * sizeof(kdata_t));
node->left = node->right = NULL;
node->index = 0;
return node;
}
/*!
* \brief free the kd-tree data structure,
*
* \param node the root node of the tree to be destroyed
*
* \param *destr a pointer to the destructor function for the data container.
*
* \return This function does not return a value
*/
void
kd_destroyTree(struct kdNode *node)
{
if ( node == NULL ) return;
kd_destroyTree(node->left);
kd_destroyTree(node->right);
kd_freeNode(node);
}
/* end of tree construction and destruction */
/* Functions dealing with result heaps */
......
src/kdtreelib/kdtree_spherical.cc
View file @ 7fb217ca
......@@ -103,13 +103,19 @@ kd_sph_orth_dist(kdata_t *p1, kdata_t *p2, int split)
*
* \return root node of the tree
*/
struct kdNode *
kdTree_t *
kd_sph_buildTree(struct kd_point *points, size_t nPoints,
kdata_t *min, kdata_t *max, int max_threads)
{
kdTree_t *tree = (kdTree_t *) malloc(sizeof(kdTree_t));
#ifdef MEMPOOL
tree->nodepool = kd_nodepool_new(nPoints*2);
#endif
struct kd_thread_data my_data;
kd_initArg(&my_data, points, nPoints, min, max, 0,max_threads, 2);
return (kdNode *)kd_doBuildTree(&my_data);
kd_initArg(&my_data, tree->nodepool, points, nPoints, min, max, 0,max_threads, 2);
tree->node = (kdNode *)kd_doBuildTree(&my_data);
return tree;
}
......@@ -176,7 +182,7 @@ kd_sph_isRectInRect(struct kdNode *node, kdata_t *min, kdata_t *max)
point->point[1] = node->min[1];
tmpNode = kd_allocNode(point, 0, point->point, point->point, 0, 2);
if (!kd_sph_isPointInRect(tmpNode, min, max)) {
kd_destroyTree(tmpNode);
kd_doDestroyTree(tmpNode);
free(point);
return 0;
}
......@@ -184,7 +190,7 @@ kd_sph_isRectInRect(struct kdNode *node, kdata_t *min, kdata_t *max)
point->point[1] = node->max[1];
tmpNode = kd_allocNode(point, 0, point->point, point->point, 0, 2);
if (!kd_sph_isPointInRect(tmpNode, min, max)) {
kd_destroyTree(tmpNode);
kd_doDestroyTree(tmpNode);
free(point);
return 0;
}
......@@ -192,7 +198,7 @@ kd_sph_isRectInRect(struct kdNode *node, kdata_t *min, kdata_t *max)
point->point[1] = node->min[1];
tmpNode = kd_allocNode(point, 0, point->point, point->point, 0, 2);
if (!kd_sph_isPointInRect(tmpNode, min, max)) {
kd_destroyTree(tmpNode);
kd_doDestroyTree(tmpNode);
free(point);
return 0;
}
......@@ -200,11 +206,11 @@ kd_sph_isRectInRect(struct kdNode *node, kdata_t *min, kdata_t *max)
point->point[1] = node->max[1];
tmpNode = kd_allocNode(point, 0, point->point, point->point, 0, 2);
if (!kd_sph_isPointInRect(tmpNode, min, max)) {
kd_destroyTree(tmpNode);
kd_doDestroyTree(tmpNode);
free(point);
return 0;
}
kd_destroyTree(tmpNode);
kd_doDestroyTree(tmpNode);
free(point);
return 1;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment