| /* --------------------------------------------------------------------- |
| Author: H. Carter Edwards |
| hcedwar@sandia.gov |
| |
| Copyright: Copyright (C) 1997 H. Carter Edwards |
| Graduate Student |
| University of Texas |
| |
| Re-release: Copyright (C) 2011-2012 H. Carter Edwards |
| |
| Purpose: Domain paritioning based upon Hilbert Space-Filling Curve |
| ordering. |
| |
| License: Re-release under the less-restrictive CLAMR software terms. |
| Permitted by email with H. Carter Edwards on 9/13/2011 |
| |
| Disclaimer: |
| |
| These routines comes with ABSOLUTELY NO WARRANTY; |
| This is free software, and you are welcome to redistribute it |
| under certain conditions. See License terms in file 'LICENSE'. |
| --------------------------------------------------------------------- */ |
| |
| #include <limits.h> |
| #include <stdlib.h> |
| |
| #include "hsfc.h" |
| |
| /*--------------------------------------------------------------------*/ |
| /* Make it callable from FORTRAN: |
| * Interface types: INTEGER and REAL*8 |
| */ |
| |
| void hsfc2sort( |
| const int N , /* IN: Number of points */ |
| const double * X , /* IN: array of X-Coordinates */ |
| const double * Y , /* IN: array of Y-Coordinates */ |
| const int ibase, /* 0 for C and 1 for Fortran */ |
| int * Info , /* OUT: (1 <= LDInfo) [ HSFC ordering ] |
| (2 <= LDInfo) [ HSFC index, #1 ] |
| (3 <= LDInfo) [ HSFC index, #2 ] */ |
| int LDInfo /* IN: Leading dimension of Info */ |
| ); |
| |
| /*--------------------------------------------------------------------*/ |
| |
| #define MaxBits ( sizeof(unsigned) * CHAR_BIT ) |
| |
| #define NBITC (32) /* 32 Bits per coordinate, resolve data at 2^31 */ |
| #define NKEY(ND) ((NBITC * ND + MaxBits - 1) / MaxBits) |
| |
| /*--------------------------------------------------------------------*/ |
| |
| static int ui1comp( const void * const I1 , const void * const I2 ) |
| { |
| return ( |
| ( ((const unsigned *)I1)[0] != ((const unsigned *)I2)[0] ) ? ( |
| ( ((const unsigned *)I1)[0] < ((const unsigned *)I2)[0] ) ? -1 : 1 ) : ( |
| 0 )); |
| } |
| |
| static int ui2comp( const void * const I1 , const void * const I2 ) |
| { |
| return ( |
| ( ((const unsigned *)I1)[0] != ((const unsigned *)I2)[0] ) ? ( |
| ( ((const unsigned *)I1)[0] < ((const unsigned *)I2)[0] ) ? -1 : 1 ) : ( |
| ( ((const unsigned *)I1)[1] != ((const unsigned *)I2)[1] ) ? ( |
| ( ((const unsigned *)I1)[1] < ((const unsigned *)I2)[1] ) ? -1 : 1 ) : ( |
| 0 ))); |
| } |
| |
| /*--------------------------------------------------------------------*/ |
| |
| static int ui3comp( const void * const I1 , const void * const I2 ) |
| { |
| return ( |
| ( ((const unsigned *)I1)[0] != ((const unsigned *)I2)[0] ) ? ( |
| ( ((const unsigned *)I1)[0] < ((const unsigned *)I2)[0] ) ? -1 : 1 ) : ( |
| ( ((const unsigned *)I1)[1] != ((const unsigned *)I2)[1] ) ? ( |
| ( ((const unsigned *)I1)[1] < ((const unsigned *)I2)[1] ) ? -1 : 1 ) : ( |
| ( ((const unsigned *)I1)[2] != ((const unsigned *)I2)[2] ) ? ( |
| ( ((const unsigned *)I1)[2] < ((const unsigned *)I2)[2] ) ? -1 : 1 ) : ( |
| 0 )))); |
| } |
| |
| static int N_uiNcomp = 0 ; |
| |
| static int uiNcomp( const void * const I1 , const void * const I2 ) |
| { |
| const int N = N_uiNcomp ; |
| register int i ; |
| |
| for ( i = 0 ; i < N && |
| ((const unsigned *)I1)[i] != ((const unsigned *)I2)[i] ; ++i ); |
| |
| return ( i < N ) ? ( |
| ( ((const unsigned *)I1)[i] < ((const unsigned *)I2)[i] ) ? -1 : 1 ) : 0 ; |
| } |
| |
| /*--------------------------------------------------------------------*/ |
| |
| void hsfc2sort( |
| const int N , /* IN: Number of points */ |
| const double * X , /* IN: array of X-Coordinates */ |
| const double * Y , /* IN: array of Y-Coordinates */ |
| const int ibase, /* 0 for C and 1 for Fortran */ |
| int * Info , /* OUT: (1 <= LDInfo) [ HSFC ordering ] |
| (2 <= LDInfo) [ HSFC index, #1 ] |
| (3 <= LDInfo) [ HSFC index, #2 ] */ |
| int LDInfo )/* IN: Leading dimension of Info */ |
| { |
| /*------------------------------------------------------------------*/ |
| |
| const double imax = ((double) ~(0u)) ; |
| |
| const unsigned ldinfo = LDInfo ; |
| const unsigned long long npt = N ; |
| const unsigned nkey = NKEY(2) ; |
| const unsigned ldT = nkey + 1 ; |
| |
| unsigned * const T = (unsigned *) malloc( sizeof(unsigned) * ldT * npt ); |
| |
| int i , ix , iy , ii , it ; |
| |
| /* Fill SFC table */ |
| |
| for ( i = it = ix = iy = 0 ; (unsigned long long)i < npt ; |
| ++i , ix++ , iy++ , it += ldT ) { |
| double xy[2] ; |
| unsigned coord[2] ; |
| |
| xy[0] = X[ix] ; |
| xy[1] = Y[iy] ; |
| |
| coord[0] = xy[0] * imax ; |
| coord[1] = xy[1] * imax ; |
| |
| hsfc2d( coord , nkey , T + it ); |
| T[it+nkey] = i ; |
| } |
| |
| /* SFC Key output */ |
| |
| if ( 2 < ldinfo && 1 < nkey ) { |
| for ( ii = 1, it = 0, i = 0 ; (unsigned long long)i < npt ; ++i, ii += ldinfo, it += ldT ) { |
| Info[ii] = T[it]; |
| Info[ii+1] = T[it+1]; |
| } |
| } |
| else if ( 1 < ldinfo ) { |
| for ( ii = 1, it = 0 ,i = 0 ; (unsigned long long)i < npt ; ++i, ii += ldinfo, it += ldT ) { |
| Info[ii] = T[it] ; |
| } |
| } |
| |
| /* Sort */ |
| |
| switch ( nkey ) { |
| case 0: break ; |
| case 1: qsort( T , npt , sizeof(unsigned) * ldT , ui1comp ); break ; |
| case 2: qsort( T , npt , sizeof(unsigned) * ldT , ui2comp ); break ; |
| case 3: qsort( T , npt , sizeof(unsigned) * ldT , ui3comp ); break ; |
| default: |
| N_uiNcomp = nkey ; |
| qsort( T , npt , sizeof(unsigned) * ldT , uiNcomp ); |
| N_uiNcomp = 0 ; |
| break ; |
| } |
| |
| for (ii = 0, i = 0, it = nkey ; (unsigned long long)i < npt ; ++i, ii += ldinfo, it += ldT) { |
| Info[ii] = T[it] + ibase; /* 1 -- FORTRAN convention, 0 -- C */ |
| } |
| |
| free( (void *) T ); |
| |
| return ; |
| } |
| |
| /*--------------------------------------------------------------------*/ |
| |
| void hsfc3sort( |
| const int N , /* IN: Number of points */ |
| const double * X , /* IN: array of X-Coordinates */ |
| const double * Y , /* IN: array of Y-Coordinates */ |
| const double * Z , /* IN: array of Y-Coordinates */ |
| const int ibase , /* IN: Stride for Y array */ |
| int * Info , /* OUT: (1 <= LDInfo) [ HSFC ordering ] |
| (2 <= LDInfo) [ HSFC index, #1 ] |
| (3 <= LDInfo) [ HSFC index, #2 ] |
| (4 <= LDInfo) [ HSFC index, #3 ] */ |
| int LDInfo )/* IN: Leading dimension of Info */ |
| { |
| /*------------------------------------------------------------------*/ |
| |
| const double imax = ((double) ~(0u)) ; |
| |
| const unsigned ldinfo = LDInfo ; |
| const unsigned long long npt = N ; |
| const unsigned nkey = NKEY(3) ; |
| const unsigned ldT = nkey + 1 ; |
| |
| unsigned * const T = (unsigned *) malloc( sizeof(unsigned) * ldT * npt ); |
| |
| int i , ix , iy , iz , ii , it ; |
| |
| /* Fill SFC table */ |
| |
| for ( i = it = ix = iy = iz = 0 ; (unsigned long long)i < npt ; |
| ++i , ix++ , iy++ , iz++ , it += ldT ) { |
| double xyz[3] ; |
| unsigned coord[3] ; |
| |
| xyz[0] = X[ix] ; |
| xyz[1] = Y[iy] ; |
| xyz[2] = Z[iz] ; |
| |
| coord[0] = xyz[0] * imax ; |
| coord[1] = xyz[1] * imax ; |
| coord[2] = xyz[2] * imax ; |
| |
| hsfc3d( coord , nkey , T + it ); |
| T[it+nkey] = i ; |
| } |
| |
| /* SFC Key output */ |
| |
| if ( 3 < ldinfo && 2 < nkey ) { |
| for ( ii = 1, it = 0, i = 0 ; (unsigned long long)i < npt ; ++i, ii += ldinfo, it += ldT ) { |
| Info[ii] = T[it]; |
| Info[ii+1] = T[it+1]; |
| Info[ii+2] = T[it+2]; |
| } |
| } |
| else if ( 2 < ldinfo && 1 < nkey ) { |
| for ( ii = 1, it = 0, i = 0 ; (unsigned long long)i < npt ; ++i, ii += ldinfo, it += ldT ) { |
| Info[ii] = T[it]; |
| Info[ii+1] = T[it+1]; |
| } |
| } |
| else if ( 1 < ldinfo ) { |
| for ( ii = 1, it = 0 ,i = 0 ; (unsigned long long)i < npt ; ++i, ii += ldinfo, it += ldT ) { |
| Info[ii] = T[it] ; |
| } |
| } |
| |
| /* Sort */ |
| |
| switch ( nkey ) { |
| case 0: break ; |
| case 1: qsort( T , npt , sizeof(unsigned) * ldT , ui1comp ); break ; |
| case 2: qsort( T , npt , sizeof(unsigned) * ldT , ui2comp ); break ; |
| case 3: qsort( T , npt , sizeof(unsigned) * ldT , ui3comp ); break ; |
| default: |
| N_uiNcomp = nkey ; |
| qsort( T , npt , sizeof(unsigned) * ldT , uiNcomp ); |
| N_uiNcomp = 0 ; |
| break ; |
| } |
| |
| for (ii = 0, i = 0, it = nkey ; (unsigned long long)i < npt ; ++i, ii += ldinfo, it += ldT) { |
| Info[ii] = T[it] + ibase ; /* FORTRAN convention */ |
| } |
| |
| free( (void *) T ); |
| |
| return ; |
| } |
| |