Commit 144b550d authored by Célestin Marot's avatar Célestin Marot
Browse files

initial commit

parents
cmake_minimum_required(VERSION 3.9)
project(HXT C)
include(CheckCCompilerFlag)
set(EXEC_NAME hxtDelaunay_sequential)
add_executable(${EXEC_NAME}
./src/hxt_vertices.c
./src/hxt_tetrahedra.c
./src/predicates.c
./src/main.c
)
# add link-time optimization -flto
# set_property(TARGET ${EXEC_NAME} PROPERTY INTERPROCEDURAL_OPTIMIZATION True)
# set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -s -fvisibility=hidden") #to avoid any function names in the executable
if (CMAKE_C_COMPILER_ID STREQUAL "GNU" OR CMAKE_C_COMPILER_ID STREQUAL "Clang")
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Ofast -march=native -ffast-math")
set_source_files_properties(./src/predicates.c PROPERTIES COMPILE_FLAGS "-fno-unsafe-math-optimizations -ffp-contract=off")
endif()
if (CMAKE_C_COMPILER_ID STREQUAL "Intel")
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Ofast -xHost -std=c99 -ansi-alias -restrict -fp-model fast")
# set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -D_Float128=__float128") # a bug in my icc
set_source_files_properties(./src/predicates.c PROPERTIES COMPILE_FLAGS "-fp-model strict")
endif()
#debug
set(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -Wall -g -DDEBUG")
This diff is collapsed.
HXTDelaunay sequential version
==============================
This folder contain a sequential version of our Delaunay algorithm.
A better parallel version is provided in [Gmsh][1], in the `contrib/hxt` folder.
However, the parallel version on a unique thread is slower than this version.
With this source code, we provide a way for people to verify that the sequential performances presented in
[our paper][2] are not faked.
Building
--------
mkdir build; cd build; cmake ..; make -j4
---
see www.hextreme.eu for more information
[1]: http://gmsh.info
[2]: https://arxiv.org/abs/1805.08831
\ No newline at end of file
This diff is collapsed.
/* This file is part of hxt_SeqDel. *
*
hxt_SeqDel is free software: you can redistribute it and/or modify *
it under the terms of the GNU General Public License as published by *
the Free Software Foundation, either version 3 of the License, or *
(at your option) any later version. *
*
hxt_SeqDel is distributed in the hope that it will be useful, *
but WITHOUT ANY WARRANTY; without even the implied warranty of *
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
GNU General Public License for more details. *
*
You should have received a copy of the GNU General Public License *
along with hxt_SeqDel. If not, see <http://www.gnu.org/licenses/>. *
*
See the COPYING file for the GNU General Public License . *
*
Author: Célestin Marot (celestin.marot@uclouvain.be) */
#ifndef _HXT_TETRAHEDRA_
#define _HXT_TETRAHEDRA_
#include "hxt_tools.h"
#define HXT_GHOST_VERTEX UINT32_MAX
typedef struct vertex_t_struct{
double coord[3];
uint64_t dist;
} vertex_t;
typedef struct bbox_t_struct{
double hxt_declare_aligned min[3];
double hxt_declare_aligned max[3];
} bbox_t;
typedef struct HXT_mesh_struct {
// vertices
uint32_t num_vertices;
uint32_t size_vertices;
vertex_t* vertices;
// bounding box
bbox_t bbox;
// tetrahedra
struct {
uint32_t* node;
uint64_t* neigh;
double* subdet;
uint64_t num;
uint64_t size;
uint32_t num_vertices;
} tetrahedra;
} mesh_t;
status_t HXT_mesh_create ( mesh_t** mesh);
status_t HXT_mesh_delete ( mesh_t** mesh);
// add vertices from mesh->vertices[from] to mesh->vertices[to-1] to the triangulation
status_t HXT_tetrahedra_compute(mesh_t* mesh);
#endif
\ No newline at end of file
/* This file is part of hxt_SeqDel. *
*
hxt_SeqDel is free software: you can redistribute it and/or modify *
it under the terms of the GNU General Public License as published by *
the Free Software Foundation, either version 3 of the License, or *
(at your option) any later version. *
*
hxt_SeqDel is distributed in the hope that it will be useful, *
but WITHOUT ANY WARRANTY; without even the implied warranty of *
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
GNU General Public License for more details. *
*
You should have received a copy of the GNU General Public License *
along with hxt_SeqDel. If not, see <http://www.gnu.org/licenses/>. *
*
See the COPYING file for the GNU General Public License . *
*
Author: Célestin Marot (celestin.marot@uclouvain.be) */
#ifndef HXT_TOOLS_H
#define HXT_TOOLS_H
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <stdarg.h>
/* HEXTREME FUNCTIONS ONLY RETURN A STATUS (except HXT_get_message_string)*/
typedef enum status_enum status_t;
enum status_enum
{
HXT_STATUS_OK = 0,
HXT_STATUS_FAILED = 1,
HXT_STATUS_OUT_OF_MEMORY = 2,
HXT_STATUS_ERROR = 3,
HXT_STATUS_FILE_CANNOT_BE_OPENED = 4,
HXT_STATUS_ASSERTION_FAILED = 5,
HXT_STATUS_POINTER_ERROR = 6
};
#define HXT_INFO(fmt,...) printf("Info : " fmt "\n", ## __VA_ARGS__ );
#define HXT_WARNING(fmt,...) fprintf(stderr,"/!\\ Warning : " fmt "\n", ## __VA_ARGS__ );
#define HXT_ERROR(status) HXT_ERROR_MSG(status,"")
#ifdef DEBUG
/* print an error message corresponding to a status */
#define HXT_ERROR_MSG(status, fmt,...) HXT_Message_Error(status, __func__, __FILE__, __LINE__, fmt, ## __VA_ARGS__ )
/* give trace message if status is not OK, but does not return */
#define HXT_TRACE(status) HXT_Message_Trace_Error(status, __func__, __FILE__, __LINE__)
/* give trace message and return the error status if not ok */
#define HXT_CHECK(status) do{ status_t s = status; if(s!=HXT_STATUS_OK) {HXT_Message_Trace_Error(s, __func__, __FILE__, __LINE__); return s; }}while(0)
/* use to check some expression inside of function, throw error if exp is not true */
#define HXT_ASSERT_MSG(exp, fmt, ...) do{ if(!(exp)) { return HXT_ERROR_MSG(HXT_STATUS_ASSERTION_FAILED, fmt, ## __VA_ARGS__ ); }} while(0)
#define HXT_ASSERT(exp) HXT_ASSERT_MSG(exp, "assertion (" #exp ") Failed")
static inline status_t HXT_Message_Error(status_t status, const char* func, const char* file, const int line, const char* fmt, ...){
if(status==HXT_STATUS_OK)
return status;
fprintf(stderr, "= X = Error (in %s -> %s:%d) :\n", func, file, line);
va_list args;
va_start(args, fmt);
vfprintf(stderr, fmt, args);
va_end (args);
return status;
}
static inline status_t HXT_Message_Trace_Error(status_t status, const char* func, const char* file, const int line){
if(status==HXT_STATUS_OK)
return status;
fprintf(stderr, " - trace - (in %s -> %s:%d) :\n", func, file, line);
return status;
}
#else
/* print an error message corresponding to a status */
#define HXT_ERROR_MSG(status, fmt,...) HXT_Message_Error(status, fmt, ## __VA_ARGS__ )
/* give trace message and return the error status if not ok */
#define HXT_CHECK(status) do{ status_t _s = status; if(_s!=HXT_STATUS_OK){ return _s; }}while(0)
/* use to check some expression inside of function, throw error if exp is not true */
#define HXT_ASSERT_MSG(exp, fmt, ...)
#define HXT_ASSERT(exp)
static inline status_t HXT_Message_Error(status_t status, const char* fmt, ...){
fprintf(stderr, "= X = Error :\n");
va_list args;
va_start(args, fmt);
vfprintf(stderr, fmt, args);
va_end (args);
return status;
}
#endif
// declare alignement of pointer allocated on the stack or in a struct
#ifdef HXT_MICROSOFT
#define hxt_declare_aligned __declspec(align(64))
#ifndef __restrict__
#define __restrict__ __restrict
#endif
// no attribute in MSVC
//#ifndef __attribute__
//#define __attribute__(x)
//#endif
#else
#define hxt_declare_aligned __attribute__((aligned(64)))
#endif
/*********************************************************
* Hextreme malloc implementation
*********************************************************/
static inline status_t HXT_malloc(void* ptr_to_ptr, size_t size)
{
void** p = (void**)ptr_to_ptr;
*p = malloc(size);
if (*p==NULL)
return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
return HXT_STATUS_OK;
}
static inline void HXT_free(void* ptr_to_ptr)
{
void** p = (void**)ptr_to_ptr;
free(*p);
*p=NULL;
}
static inline status_t HXT_realloc(void* ptr_to_ptr, size_t size)
{
void** p = (void**)ptr_to_ptr;
void* newptr = realloc(*p, size);
if (newptr==NULL && *p!=NULL && size!=0)
return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
*p = newptr;
return HXT_STATUS_OK;
}
#endif
/* This file is part of hxt_SeqDel. *
*
hxt_SeqDel is free software: you can redistribute it and/or modify *
it under the terms of the GNU General Public License as published by *
the Free Software Foundation, either version 3 of the License, or *
(at your option) any later version. *
*
hxt_SeqDel is distributed in the hope that it will be useful, *
but WITHOUT ANY WARRANTY; without even the implied warranty of *
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
GNU General Public License for more details. *
*
You should have received a copy of the GNU General Public License *
along with hxt_SeqDel. If not, see <http://www.gnu.org/licenses/>. *
*
See the COPYING file for the GNU General Public License . *
*
Author: Célestin Marot (celestin.marot@uclouvain.be) */
#include <string.h>
#include "hxt_vertices.h"
/* create a nextbefore macro
* for a strictly positive value, round to previous double value */
#if (defined (__STD_VERSION__) && (__STD_VERSION__ >= 199901L)) // nextafter only from C99
#include <math.h>
#define nextbefore(x) nextafter(x,0.0);
#else
#include <float.h>
#define nextbefore(x) (x*(1.0-DBL_EPSILON))
#endif
#define MIN(x,y) ((x)<(y)?(x):(y))
#define SWAP(x,y) do{uint32_t tmp=x; x=y; y=tmp;}while(0)
#define INVE(x) x=s-1-x
uint32_t HXT_hilbert_advised_bits(uint32_t n)
{
uint32_t nlog2 = (32-__builtin_clz(n))*3/2; // if we have n points, we want approximately n^1.5 possible hilbert coordinates
return MIN(63, nlog2>15?(nlog2+10)/11*11:nlog2); // we prefer multiple of 11
}
/*================================= compute the hilberts distance ===========================================
*
* This part is for computing the hilbert coordinates of the vertices (X,Y,Z)
*/
static status_t HXT_vertices_hilbert_dist(bbox_t* bbox, vertex_t* vertices, uint32_t n, uint32_t* nbits)
{
uint32_t level;
if(*nbits==0){
*nbits = HXT_hilbert_advised_bits(n);
}
else if(*nbits>63){
*nbits = 63;
}
level = (*nbits+2)/3;
*nbits = level*3;
uint32_t nmax = 1U<<level;
double divX = nextbefore(nmax/(bbox->max[0]-bbox->min[0]));
double divY = nextbefore(nmax/(bbox->max[1]-bbox->min[1]));
double divZ = nextbefore(nmax/(bbox->max[2]-bbox->min[2]));
for (uint32_t i=0; i<n; i++)
{
uint32_t x = (vertices[i].coord[0]-bbox->min[0])*divX;
uint32_t y = (vertices[i].coord[1]-bbox->min[1])*divY;
uint32_t z = (vertices[i].coord[2]-bbox->min[2])*divZ;
uint64_t bits = 0;
// #ifdef DEBUG
// if(x<nmax && y<nmax && z<nmax,"coordinate out of bbox")
// i=n+1;
// #endif
uint32_t s;
for (s=nmax>>1; s; s>>=1) {
uint32_t rx = (x & s) != 0;
uint32_t ry = (y & s) != 0;
uint32_t rz = (z & s) != 0;
uint32_t gc = (7*rx)^((3*ry)^rz);
// __assume(gc<=7);
bits += (uint64_t) s*s*s* gc;
if(gc & 4){
if(gc & 2){
if(gc & 1){ // 7
SWAP(x,z);
INVE(x);
INVE(z);
}
else{ // 6
SWAP(x,y);
INVE(x);
INVE(y);
}
}
else{
if(!(gc & 1)){ // 4
SWAP(x,z);
SWAP(y,z);
}
}
}
else{
if(gc & 2){
if(gc & 1){ // 3
SWAP(x,z);
SWAP(y,z);
INVE(x);
INVE(y);
}
}
else{
if(gc & 1) // 1
SWAP(x,y);
else // 0
SWAP(x,z);
}
}
}
vertices[i].dist = bits;
}
return HXT_STATUS_OK;
}
/************************************** sort functions ***********************************************/
// quickest sort for less than 32 entries
static inline void HXT_insertion_sort(vertex_t* array, uint32_t from, uint32_t to){
uint32_t r;
for (r=from+1; r<to; r++)
{
vertex_t tmp = array[r];
uint64_t a = tmp.dist;
uint32_t l;
for (l=r; l>from && array[l-1].dist>a; l--){
array[l] = array[l-1];
}
array[l] = tmp;
}
}
// quickest sort for less than 10000 entries
static void HXT_radix2_sort_in_place(vertex_t* array, uint32_t from, uint32_t to, uint64_t mask)
{
uint32_t l = from, r;
while(mask && to > from + 32){
r = to - 1;
while (1) {
/* find left most with bit, and right most without bit, swap */
while (l < r && (array[l].dist & mask)==0) l++;
while (l < r && (array[r].dist & mask)) r--;
if (l >= r) break;
// swap
vertex_t tmp = array[l];
array[l] = array[r];
array[r] = tmp;
}
if (!(mask & array[l].dist) && l < to) l++;
mask >>= 1;
HXT_radix2_sort_in_place(array, from, l, mask);
from=l;
}
HXT_insertion_sort(array, from, to);
}
/********************************************/
// exclusive scan prefix sum
static inline void scan2048(uint32_t* h)
{
unsigned i;
uint32_t sum = 0;
for (i=0; i<2048; i++) {
uint32_t tsum = h[i] + sum;
h[i] = sum; // sum of the preceding, but not with itself !!
sum = tsum;
}
}
static inline void rad11_lsb(vertex_t* __restrict__ array, vertex_t* __restrict__ buffer, const uint32_t n, const unsigned npass)
{
unsigned pass;
for (pass=0; pass<npass; pass++)
{
uint32_t hxt_declare_aligned h[2048] = {0};
for (uint32_t i=0; i<n; i++)
h[(array[i].dist >> (11*pass)) & 2047]++;
scan2048(h);
for (uint32_t i=0; i<n; i++)
{
vertex_t ai = array[i];
buffer[h[(ai.dist >> (11*pass)) & 2047]++] = ai;
}
vertex_t* tmp = array;
array = buffer;
buffer = tmp;
}
}
// quickest sort for over 10000 entries
static status_t HXT_vertices_sort_no_copy(vertex_t* const __restrict__ array, vertex_t* const __restrict__ buffer, const uint32_t n, uint32_t nbits, int* to_copy)
{
uint32_t hxt_declare_aligned hist[2048] = {0};
const unsigned bits_rem = (nbits<=11)?0:nbits-11; // bits remaining after first pass
const unsigned pass_rem = (bits_rem+10)/11;
*to_copy = (pass_rem&1)==0;
/***** we first make a MSB radix11 sort pass to ensure locality of data ***/
for (uint32_t i=0; i<n; i++)
hist[array[i].dist>>bits_rem & 2047]++;
// calculate total sums
scan2048(hist);
// read/write histogram, copy
for (uint32_t i=0; i < n; i++) {
vertex_t ai = array[i];
buffer[hist[ai.dist>>bits_rem & 2047]++] = ai;
}
/***** then we make a LSB radix sort pass with each of the 2048 pieces of array ***/
uint32_t start = 0;
for (uint32_t i=0; i<2048; i++)
{
uint32_t end = hist[i];
rad11_lsb(buffer + start, array + start, end-start, pass_rem);
start = end;
}
return HXT_STATUS_OK;
}
/* wrapper over HXT_vertices_sort_no_copy that copy the result in the initial array if needed
+ use HXT_radix_sort2_in_place if the array is small */
static status_t HXT_vertices_sort(vertex_t* const __restrict__ array, const uint32_t n, uint32_t nbits)
{
if(nbits==0){
return HXT_STATUS_OK;
}
if(n<=16384){
if(nbits>64){
nbits = 64;
}
uint64_t one = 1;
HXT_radix2_sort_in_place(array, 0, n, one<<(nbits-1));
return HXT_STATUS_OK;
}
else{
vertex_t* buffer = NULL;
HXT_CHECK(
HXT_malloc(&buffer, n*sizeof(vertex_t)) );
int to_copy = 0;
HXT_CHECK(
HXT_vertices_sort_no_copy(array, buffer, n, nbits, &to_copy) );
if(to_copy){
memcpy(array, buffer, n*sizeof(vertex_t));
}
HXT_free(&buffer);
return HXT_STATUS_OK;
}
}
/*********************************** shuffle functions ***********************************************/
static inline uint32_t fast_hash(uint32_t x) {
x = ((x >> 16) ^ x) * 0x45d9f3b;
x = ((x >> 16) ^ x) * 0x45d9f3b;
x = (x >> 16) ^ x;
return x;
}
/****************************************** BRIO *****************************************************/
/* automatic biased ranomized insertion order */
status_t HXT_vertices_BRIO(bbox_t* bbox, vertex_t* vertices, const uint32_t n){
for (uint32_t i=0; i<n; i++){
vertices[i].dist = fast_hash(i);
}
HXT_CHECK(
HXT_vertices_sort(vertices, n, 22) ); // just make two pass of the sort (no copy needed)
uint32_t nbits=0;
HXT_CHECK(
HXT_vertices_hilbert_dist(bbox, vertices, n, &nbits) );
// sort the 32768 first number, than, because there will be approximately 7* more tetrahedron, make it times 7 every time
uint32_t end = n;
// int to_copy = 1;
while(end > 1500) {
uint32_t start = end/7.5;
HXT_CHECK(
HXT_vertices_sort(vertices+start, end-start, nbits) );
end = start;
}
for (uint32_t i=0; i<n; i++){
vertices[i].dist = UINT64_C(0x0000000000000000);
}
return HXT_STATUS_OK;
}
/* This file is part of hxt_SeqDel. *
*
hxt_SeqDel is free software: you can redistribute it and/or modify *
it under the terms of the GNU General Public License as published by *
the Free Software Foundation, either version 3 of the License, or *
(at your option) any later version. *
*
hxt_SeqDel is distributed in the hope that it will be useful, *
but WITHOUT ANY WARRANTY; without even the implied warranty of *
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
GNU General Public License for more details. *
*
You should have received a copy of the GNU General Public License *
along with hxt_SeqDel. If not, see <http://www.gnu.org/licenses/>. *
*
See the COPYING file for the GNU General Public License . *
*
Author: Célestin Marot (celestin.marot@uclouvain.be) */
#ifndef _HEXTREME_VERTICES_