1 /* ND.c */
2
3 #include "../misc.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 ---------------------------------------------------------------
8 this file contains procedures to create a nested dissection
9 permutation in one, two or three dimensions.
10 ---------------------------------------------------------------
11 */
12
13 #define DEBUG 0
14
15 /*--------------------------------------------------------------------*/
16 /*
17 -----------------
18 static procedures
19 -----------------
20 */
21 static void SplitX ( int n1, int n2, int n3, int newToOld[],
22 int west, int east, int south, int north,
23 int bottom, int top ) ;
24 static void SplitY ( int n1, int n2, int n3, int newToOld[],
25 int west, int east, int south, int north,
26 int bottom, int top ) ;
27 static void SplitZ ( int n1, int n2, int n3, int newToOld[],
28 int west, int east, int south, int north,
29 int bottom, int top ) ;
30 static int WhichCut ( int west, int east, int south, int north,
31 int bottom, int top ) ;
32 static int MidPoint ( int left, int right, int glob_center ) ;
33
34 /*--------------------------------------------------------------------*/
35 /*
36 ------------------------------------------------------------
37 purpose -- main procedure. if the input region is 1 x 1 x 1,
38 the node is put into the permutation vector.
39 otherwise the region is split into three pieces,
40 two subregions and a separator, and recursive
41 calls are made to order the subregions
42
43 input --
44
45 n1 -- number of points in the first direction
46 n2 -- number of points in the second direction
47 n3 -- number of points in the third direction
48 newToOld -- pointer to the permutation vector
49 west -- west coordinate
50 east -- east coordinate
51 south -- south coordinate
52 north -- north coordinate
53 bottom -- bottom coordinate
54 top -- top coordinate
55
56 created -- 95nov15, cca
57 ------------------------------------------------------------
58 */
59 void
mkNDperm(int n1,int n2,int n3,int newToOld[],int west,int east,int south,int north,int bottom,int top)60 mkNDperm (
61 int n1,
62 int n2,
63 int n3,
64 int newToOld[],
65 int west,
66 int east,
67 int south,
68 int north,
69 int bottom,
70 int top
71 ) {
72 if ( n1 <= 0 || n2 <= 0 || n3 <= 0 || newToOld == NULL
73 || west < 0 || east >= n1
74 || south < 0 || north >= n2
75 || bottom < 0 || top >= n3 ) {
76 fprintf(stderr,
77 "\n fatal error in mkNDperm(%d,%d,%d,%p,%d,%d,%d,%d,%d,%d)"
78 "\n bad input data\n",
79 n1, n2, n3, newToOld,
80 west, east, south, north, bottom, top) ;
81 exit(-1) ;
82 }
83 if ( west == east && south == north && bottom == top ) {
84 # if DEBUG > 0
85 fprintf(stdout, "\n ND : ordering %d",
86 west + south * n1 + bottom * n1 * n2) ;
87 # endif
88 newToOld[0] = west + south * n1 + bottom * n1 * n2 ; }
89 else {
90 switch ( WhichCut(west, east, south, north, bottom, top) ) {
91 case 1 :
92 # if DEBUG > 0
93 fprintf(stdout, "\n ND : calling Split9X") ;
94 # endif
95 SplitX(n1, n2, n3, newToOld,
96 west, east, south, north, bottom, top) ;
97 break ;
98 case 2 :
99 # if DEBUG > 0
100 fprintf(stdout, "\n ND : calling Split9Y") ;
101 # endif
102 SplitY(n1, n2, n3, newToOld,
103 west, east, south, north, bottom, top) ;
104 break ;
105 case 3 :
106 # if DEBUG > 0
107 fprintf(stdout, "\n ND : calling Split9Z") ;
108 # endif
109 SplitZ(n1, n2, n3, newToOld,
110 west, east, south, north, bottom, top) ;
111 break ; } }
112
113 return ; }
114
115 /*--------------------------------------------------------------------*/
116 /*
117 -------------------------------------------------------
118 purpose -- to decide which way to partition a subregion
119
120 created -- 95nov16, cca
121 -------------------------------------------------------
122 */
123 static int
WhichCut(int west,int east,int south,int north,int bottom,int top)124 WhichCut (
125 int west,
126 int east,
127 int south,
128 int north,
129 int bottom,
130 int top
131 ) {
132 int d1, d2, d3 ;
133
134 d1 = east - west + 1 ;
135 d2 = north - south + 1 ;
136 d3 = top - bottom + 1 ;
137 #if DEBUG > 0
138 fprintf(stdout, "\n WhichCut : (d1,d2,d3) = (%d,%d,%d)", d1,d2,d3) ;
139 #endif
140 if ( d1 > d2 && d1 > d3 ) {
141 return(1) ; }
142 else if ( d2 > d1 && d2 > d3 ) {
143 return(2) ; }
144 else if ( d3 > d1 && d3 > d2 ) {
145 return(3) ; }
146 else if ( d1 == d2 && d1 > d3 ) {
147 return(1) ; }
148 else if ( d1 == d3 && d1 > d2 ) {
149 return(1) ; }
150 else if ( d2 == d3 && d2 > d1 ) {
151 return(2) ; }
152 else {
153 return(3) ; } }
154
155 /*--------------------------------------------------------------------*/
156 /*
157 ------------------------------------------------
158 purpose -- to split a subregion with a Y-Z plane
159 ------------------------------------------------
160 */
161 static void
SplitX(int n1,int n2,int n3,int newToOld[],int west,int east,int south,int north,int bottom,int top)162 SplitX (
163 int n1,
164 int n2,
165 int n3,
166 int newToOld[],
167 int west,
168 int east,
169 int south,
170 int north,
171 int bottom,
172 int top
173 ) {
174 int east1, east2, j, k, k1, k2, k3, m, m1, west1, west2 ;
175
176 m1 = MidPoint(west, east, n1/2) ;
177 west1 = west ;
178 east1 = m1 - 1 ;
179 west2 = m1 + 1 ;
180 east2 = east ;
181 k1 = 0 ;
182 k2 = k1 + (m1 - west) * (north - south + 1) * (top - bottom + 1) ;
183 k3 = k2 + (east - m1) * (north - south + 1) * (top - bottom + 1) ;
184 if ( m1 > west ) {
185 # if DEBUG > 0
186 fprintf(stdout, "\n SplitX : 1. calling ND") ;
187 # endif
188 mkNDperm(n1, n2, n3, newToOld + k1,
189 west1, east1, south, north, bottom, top) ; }
190 if ( east > m1 ) {
191 # if DEBUG > 0
192 fprintf(stdout, "\n SplitX : 2. calling ND") ;
193 # endif
194 mkNDperm(n1, n2, n3, newToOld + k2,
195 west2, east2, south, north, bottom, top) ; }
196 # if DEBUG > 0
197 fprintf(stdout, "\n SplitX : 3. calling ND") ;
198 fprintf(stdout, "\n m1 = %d, south = %d, north = %d",
199 m1, south, north) ;
200 # endif
201 m = k3 ;
202 for ( k = bottom ; k <= top ; k++ ) {
203 for ( j = south ; j <= north ; j++ ) {
204 newToOld[m++] = m1 + j * n1 + k * n1 * n2 ; } }
205
206 return ; }
207
208 /*--------------------------------------------------------------------*/
209 /*
210 ------------------------------------------------
211 purpose -- to split a subregion with a X-Z plane
212
213 created -- 95nov16, cca
214 ------------------------------------------------
215 */
216 static void
SplitY(int n1,int n2,int n3,int newToOld[],int west,int east,int south,int north,int bottom,int top)217 SplitY (
218 int n1,
219 int n2,
220 int n3,
221 int newToOld[],
222 int west,
223 int east,
224 int south,
225 int north,
226 int bottom,
227 int top
228 ) {
229 int i, k, k1, k2, k3, m, m1, north1, north2, south1, south2 ;
230
231 m1 = MidPoint(south, north, n2/2) ;
232 south1 = south ;
233 north1 = m1 - 1 ;
234 south2 = m1 + 1 ;
235 north2 = north ;
236 k1 = 0 ;
237 k2 = k1 + (m1 - south) * (east - west + 1) * (top - bottom + 1) ;
238 k3 = k2 + (north - m1) * (east - west + 1) * (top - bottom + 1) ;
239 if ( m1 > south ) {
240 # if DEBUG > 0
241 fprintf(stdout, "\n SplitY : calling ND(%d:%d,%d:%d,%d:%d)",
242 west, east, south1, north1, bottom, top) ;
243 # endif
244 mkNDperm(n1, n2, n3, newToOld + k1,
245 west, east, south1, north1, bottom, top) ; }
246 if ( north > m1 ) {
247 # if DEBUG > 0
248 fprintf(stdout, "\n SplitY : calling ND(%d:%d,%d:%d,%d:%d)",
249 west, east, south2, north2, bottom, top) ;
250 # endif
251 mkNDperm(n1, n2, n3, newToOld + k2,
252 west, east, south2, north2, bottom, top) ; }
253 # if DEBUG > 0
254 fprintf(stdout, "\n SplitY : ordering (%d:%d,%d:%d,%d:%d)",
255 west, east, m1, m1, bottom, top) ;
256 # endif
257 m = k3 ;
258 for ( k = bottom ; k <= top ; k++ ) {
259 for ( i = west ; i <= east ; i++ ) {
260 # if DEBUG > 0
261 fprintf(stdout, "\n SplitY : newToOld[%d] = %d",
262 m, i + m1 * n1 + k * n1 * n2) ;
263 # endif
264 newToOld[m++] = i + m1 * n1 + k * n1 * n2 ; } }
265
266 return ; }
267
268 /*--------------------------------------------------------------------*/
269 /*
270 ------------------------------------------------
271 purpose -- to split a subregion with a X-Y plane
272
273 created -- 95nov16, cca
274 ------------------------------------------------
275 */
276 static void
SplitZ(int n1,int n2,int n3,int newToOld[],int west,int east,int south,int north,int bottom,int top)277 SplitZ (
278 int n1,
279 int n2,
280 int n3,
281 int newToOld[],
282 int west,
283 int east,
284 int south,
285 int north,
286 int bottom,
287 int top
288 ) {
289 int bottom1, bottom2, i, j, k1, k2, k3, m, m1, top1, top2 ;
290
291 m1 = MidPoint(bottom, top, n2/2) ;
292 bottom1 = bottom ;
293 top1 = m1 - 1 ;
294 bottom2 = m1 + 1 ;
295 top2 = top ;
296 k1 = 0 ;
297 k2 = k1 + (m1 - bottom) * (east - west + 1) * (north - south + 1) ;
298 k3 = k2 + (top - m1) * (east - west + 1) * (north - south + 1) ;
299 if ( m1 > bottom ) {
300 # if DEBUG > 0
301 fprintf(stdout, "\n SplitZ : 1. calling ND") ;
302 # endif
303 mkNDperm(n1, n2, n3, newToOld + k1,
304 west, east, south, north, bottom1, top1) ; }
305 if ( top > m1 ) {
306 # if DEBUG > 0
307 fprintf(stdout, "\n SplitZ : 2. calling ND") ;
308 # endif
309 mkNDperm(n1, n2, n3, newToOld + k2,
310 west, east, south, north, bottom2, top2) ; }
311 # if DEBUG > 0
312 fprintf(stdout, "\n SplitZ : 3. calling ND") ;
313 # endif
314 m = k3 ;
315 for ( j = south ; j <= north ; j++ ) {
316 for ( i = west ; i <= east ; i++ ) {
317 newToOld[m++] = i + j * n1 + m1 * n1 * n2 ; } }
318
319 return ; }
320
321 /*--------------------------------------------------------------------*/
322 /*
323 ----------------------------------------------
324 purpose -- to compute the midpoint of a region
325
326 input --
327
328 left -- left coordinate
329 right -- right coordinate
330 glob_center -- glob_center coordinate
331 ----------------------------------------------
332 */
333 static int
MidPoint(int left,int right,int glob_center)334 MidPoint (
335 int left,
336 int right,
337 int glob_center
338 ) {
339 int mid ;
340
341 mid = (left + right)/2 ;
342 if ( (left + right) % 2 == 0 ) {
343 return(mid) ; }
344 else if ( mid >= glob_center ) {
345 return(mid) ; }
346 else {
347 return(mid+1) ; } }
348
349 /*--------------------------------------------------------------------*/
350