/* ND.c */
#include "../misc.h"
/*--------------------------------------------------------------------*/
/*
---------------------------------------------------------------
this file contains procedures to create a nested dissection
permutation in one, two or three dimensions.
---------------------------------------------------------------
*/
#define DEBUG 0
/*--------------------------------------------------------------------*/
/*
-----------------
static procedures
-----------------
*/
static void SplitX ( int n1, int n2, int n3, int newToOld[],
int west, int east, int south, int north,
int bottom, int top ) ;
static void SplitY ( int n1, int n2, int n3, int newToOld[],
int west, int east, int south, int north,
int bottom, int top ) ;
static void SplitZ ( int n1, int n2, int n3, int newToOld[],
int west, int east, int south, int north,
int bottom, int top ) ;
static int WhichCut ( int west, int east, int south, int north,
int bottom, int top ) ;
static int MidPoint ( int left, int right, int glob_center ) ;
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------------
purpose -- main procedure. if the input region is 1 x 1 x 1,
the node is put into the permutation vector.
otherwise the region is split into three pieces,
two subregions and a separator, and recursive
calls are made to order the subregions
input --
n1 -- number of points in the first direction
n2 -- number of points in the second direction
n3 -- number of points in the third direction
newToOld -- pointer to the permutation vector
west -- west coordinate
east -- east coordinate
south -- south coordinate
north -- north coordinate
bottom -- bottom coordinate
top -- top coordinate
created -- 95nov15, cca
------------------------------------------------------------
*/
void
mkNDperm (
int n1,
int n2,
int n3,
int newToOld[],
int west,
int east,
int south,
int north,
int bottom,
int top
) {
if ( n1 <= 0 || n2 <= 0 || n3 <= 0 || newToOld == NULL
|| west < 0 || east >= n1
|| south < 0 || north >= n2
|| bottom < 0 || top >= n3 ) {
fprintf(stderr,
"\n fatal error in mkNDperm(%d,%d,%d,%p,%d,%d,%d,%d,%d,%d)"
"\n bad input data\n",
n1, n2, n3, newToOld,
west, east, south, north, bottom, top) ;
exit(-1) ;
}
if ( west == east && south == north && bottom == top ) {
# if DEBUG > 0
fprintf(stdout, "\n ND : ordering %d",
west + south * n1 + bottom * n1 * n2) ;
# endif
newToOld[0] = west + south * n1 + bottom * n1 * n2 ; }
else {
switch ( WhichCut(west, east, south, north, bottom, top) ) {
case 1 :
# if DEBUG > 0
fprintf(stdout, "\n ND : calling Split9X") ;
# endif
SplitX(n1, n2, n3, newToOld,
west, east, south, north, bottom, top) ;
break ;
case 2 :
# if DEBUG > 0
fprintf(stdout, "\n ND : calling Split9Y") ;
# endif
SplitY(n1, n2, n3, newToOld,
west, east, south, north, bottom, top) ;
break ;
case 3 :
# if DEBUG > 0
fprintf(stdout, "\n ND : calling Split9Z") ;
# endif
SplitZ(n1, n2, n3, newToOld,
west, east, south, north, bottom, top) ;
break ; } }
return ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------------------------------
purpose -- to decide which way to partition a subregion
created -- 95nov16, cca
-------------------------------------------------------
*/
static int
WhichCut (
int west,
int east,
int south,
int north,
int bottom,
int top
) {
int d1, d2, d3 ;
d1 = east - west + 1 ;
d2 = north - south + 1 ;
d3 = top - bottom + 1 ;
#if DEBUG > 0
fprintf(stdout, "\n WhichCut : (d1,d2,d3) = (%d,%d,%d)", d1,d2,d3) ;
#endif
if ( d1 > d2 && d1 > d3 ) {
return(1) ; }
else if ( d2 > d1 && d2 > d3 ) {
return(2) ; }
else if ( d3 > d1 && d3 > d2 ) {
return(3) ; }
else if ( d1 == d2 && d1 > d3 ) {
return(1) ; }
else if ( d1 == d3 && d1 > d2 ) {
return(1) ; }
else if ( d2 == d3 && d2 > d1 ) {
return(2) ; }
else {
return(3) ; } }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------
purpose -- to split a subregion with a Y-Z plane
------------------------------------------------
*/
static void
SplitX (
int n1,
int n2,
int n3,
int newToOld[],
int west,
int east,
int south,
int north,
int bottom,
int top
) {
int east1, east2, j, k, k1, k2, k3, m, m1, west1, west2 ;
m1 = MidPoint(west, east, n1/2) ;
west1 = west ;
east1 = m1 - 1 ;
west2 = m1 + 1 ;
east2 = east ;
k1 = 0 ;
k2 = k1 + (m1 - west) * (north - south + 1) * (top - bottom + 1) ;
k3 = k2 + (east - m1) * (north - south + 1) * (top - bottom + 1) ;
if ( m1 > west ) {
# if DEBUG > 0
fprintf(stdout, "\n SplitX : 1. calling ND") ;
# endif
mkNDperm(n1, n2, n3, newToOld + k1,
west1, east1, south, north, bottom, top) ; }
if ( east > m1 ) {
# if DEBUG > 0
fprintf(stdout, "\n SplitX : 2. calling ND") ;
# endif
mkNDperm(n1, n2, n3, newToOld + k2,
west2, east2, south, north, bottom, top) ; }
# if DEBUG > 0
fprintf(stdout, "\n SplitX : 3. calling ND") ;
fprintf(stdout, "\n m1 = %d, south = %d, north = %d",
m1, south, north) ;
# endif
m = k3 ;
for ( k = bottom ; k <= top ; k++ ) {
for ( j = south ; j <= north ; j++ ) {
newToOld[m++] = m1 + j * n1 + k * n1 * n2 ; } }
return ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------
purpose -- to split a subregion with a X-Z plane
created -- 95nov16, cca
------------------------------------------------
*/
static void
SplitY (
int n1,
int n2,
int n3,
int newToOld[],
int west,
int east,
int south,
int north,
int bottom,
int top
) {
int i, k, k1, k2, k3, m, m1, north1, north2, south1, south2 ;
m1 = MidPoint(south, north, n2/2) ;
south1 = south ;
north1 = m1 - 1 ;
south2 = m1 + 1 ;
north2 = north ;
k1 = 0 ;
k2 = k1 + (m1 - south) * (east - west + 1) * (top - bottom + 1) ;
k3 = k2 + (north - m1) * (east - west + 1) * (top - bottom + 1) ;
if ( m1 > south ) {
# if DEBUG > 0
fprintf(stdout, "\n SplitY : calling ND(%d:%d,%d:%d,%d:%d)",
west, east, south1, north1, bottom, top) ;
# endif
mkNDperm(n1, n2, n3, newToOld + k1,
west, east, south1, north1, bottom, top) ; }
if ( north > m1 ) {
# if DEBUG > 0
fprintf(stdout, "\n SplitY : calling ND(%d:%d,%d:%d,%d:%d)",
west, east, south2, north2, bottom, top) ;
# endif
mkNDperm(n1, n2, n3, newToOld + k2,
west, east, south2, north2, bottom, top) ; }
# if DEBUG > 0
fprintf(stdout, "\n SplitY : ordering (%d:%d,%d:%d,%d:%d)",
west, east, m1, m1, bottom, top) ;
# endif
m = k3 ;
for ( k = bottom ; k <= top ; k++ ) {
for ( i = west ; i <= east ; i++ ) {
# if DEBUG > 0
fprintf(stdout, "\n SplitY : newToOld[%d] = %d",
m, i + m1 * n1 + k * n1 * n2) ;
# endif
newToOld[m++] = i + m1 * n1 + k * n1 * n2 ; } }
return ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------
purpose -- to split a subregion with a X-Y plane
created -- 95nov16, cca
------------------------------------------------
*/
static void
SplitZ (
int n1,
int n2,
int n3,
int newToOld[],
int west,
int east,
int south,
int north,
int bottom,
int top
) {
int bottom1, bottom2, i, j, k1, k2, k3, m, m1, top1, top2 ;
m1 = MidPoint(bottom, top, n2/2) ;
bottom1 = bottom ;
top1 = m1 - 1 ;
bottom2 = m1 + 1 ;
top2 = top ;
k1 = 0 ;
k2 = k1 + (m1 - bottom) * (east - west + 1) * (north - south + 1) ;
k3 = k2 + (top - m1) * (east - west + 1) * (north - south + 1) ;
if ( m1 > bottom ) {
# if DEBUG > 0
fprintf(stdout, "\n SplitZ : 1. calling ND") ;
# endif
mkNDperm(n1, n2, n3, newToOld + k1,
west, east, south, north, bottom1, top1) ; }
if ( top > m1 ) {
# if DEBUG > 0
fprintf(stdout, "\n SplitZ : 2. calling ND") ;
# endif
mkNDperm(n1, n2, n3, newToOld + k2,
west, east, south, north, bottom2, top2) ; }
# if DEBUG > 0
fprintf(stdout, "\n SplitZ : 3. calling ND") ;
# endif
m = k3 ;
for ( j = south ; j <= north ; j++ ) {
for ( i = west ; i <= east ; i++ ) {
newToOld[m++] = i + j * n1 + m1 * n1 * n2 ; } }
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the midpoint of a region
input --
left -- left coordinate
right -- right coordinate
glob_center -- glob_center coordinate
----------------------------------------------
*/
static int
MidPoint (
int left,
int right,
int glob_center
) {
int mid ;
mid = (left + right)/2 ;
if ( (left + right) % 2 == 0 ) {
return(mid) ; }
else if ( mid >= glob_center ) {
return(mid) ; }
else {
return(mid+1) ; } }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1