1 /* ND2.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 note, all separators are double wide
12
13 created -- 96feb01, cca
14 ---------------------------------------------------------------
15 */
16
17 #define DEBUG 0
18
19 /*--------------------------------------------------------------------*/
20 /*
21 -----------------
22 static procedures
23 -----------------
24 */
25 static void SplitX ( int n1, int n2, int n3, int newToOld[],
26 int west, int east, int south, int north,
27 int bottom, int top ) ;
28 static void SplitY ( int n1, int n2, int n3, int newToOld[],
29 int west, int east, int south, int north,
30 int bottom, int top ) ;
31 static void SplitZ ( int n1, int n2, int n3, int newToOld[],
32 int west, int east, int south, int north,
33 int bottom, int top ) ;
34 static int WhichCut ( int west, int east, int south, int north,
35 int bottom, int top ) ;
36 static int MidPoint ( int left, int right, int glob_center ) ;
37
38 /*--------------------------------------------------------------------*/
39 /*
40 ------------------------------------------------------------
41 purpose -- main procedure. if the input region is m1 x m2 x m3,
42 where 0 < m1, m2, m3 < 2, then
43 the node is put into the permutation vector.
44 otherwise the region is split into three pieces,
45 two subregions and a separator, and recursive
46 calls are made to order the subregions
47
48 input --
49
50 n1 -- number of points in the first direction
51 n2 -- number of points in the second direction
52 n3 -- number of points in the third direction
53 newToOld -- pointer to the permutation vector
54 west -- west coordinate
55 east -- east coordinate
56 south -- south coordinate
57 north -- north coordinate
58 bottom -- bottom coordinate
59 top -- top coordinate
60
61 created -- 95nov15, cca
62 ------------------------------------------------------------
63 */
64 void
mkNDperm2(int n1,int n2,int n3,int newToOld[],int west,int east,int south,int north,int bottom,int top)65 mkNDperm2 (
66 int n1,
67 int n2,
68 int n3,
69 int newToOld[],
70 int west,
71 int east,
72 int south,
73 int north,
74 int bottom,
75 int top
76 ) {
77 int count, i, j, k ;
78 if ( n1 <= 0 || n2 <= 0 || n3 <= 0 || newToOld == NULL
79 || west < 0 || east >= n1
80 || south < 0 || north >= n2
81 || bottom < 0 || top >= n3 ) {
82 fprintf(stderr,
83 "\n fatal error in mkND2perm(%d,%d,%d,%p,%d,%d,%d,%d,%d,%d)"
84 "\n bad input data\n",
85 n1, n2, n3, newToOld,
86 west, east, south, north, bottom, top) ;
87 exit(-1) ;
88 }
89
90 if ( east - west <= 1 && north - south <= 1 && top - bottom <= 1 ) {
91 count = 0 ;
92 for ( i = west ; i <= east ; i++ ) {
93 for ( j = south ; j <= north ; j++ ) {
94 for ( k = bottom ; k <= top ; k++ ) {
95 # if DEBUG > 0
96 fprintf(stdout, "\n ND : ordering %d",
97 i + j * n1 + k * n1 * n2) ;
98 # endif
99 newToOld[count++] = i + j * n1 + k * n1 * n2 ;
100 }
101 }
102 }
103 } else {
104 switch ( WhichCut(west, east, south, north, bottom, top) ) {
105 case 1 :
106 # if DEBUG > 0
107 fprintf(stdout, "\n ND : calling SplitX(%d,%d,%d,%d,%d,%d)",
108 west, east, south, north, bottom, top) ;
109 # endif
110 SplitX(n1, n2, n3, newToOld,
111 west, east, south, north, bottom, top) ;
112 break ;
113 case 2 :
114 # if DEBUG > 0
115 fprintf(stdout, "\n ND : calling SplitY(%d,%d,%d,%d,%d,%d)",
116 west, east, south, north, bottom, top) ;
117 # endif
118 SplitY(n1, n2, n3, newToOld,
119 west, east, south, north, bottom, top) ;
120 break ;
121 case 3 :
122 # if DEBUG > 0
123 fprintf(stdout, "\n ND : calling SplitZ(%d,%d,%d,%d,%d,%d)",
124 west, east, south, north, bottom, top) ;
125 # endif
126 SplitZ(n1, n2, n3, newToOld,
127 west, east, south, north, bottom, top) ;
128 break ; } }
129
130 return ; }
131
132 /*--------------------------------------------------------------------*/
133 /*
134 -------------------------------------------------------
135 purpose -- to decide which way to partition a subregion
136
137 created -- 95nov16, cca
138 -------------------------------------------------------
139 */
140 static int
WhichCut(int west,int east,int south,int north,int bottom,int top)141 WhichCut (
142 int west,
143 int east,
144 int south,
145 int north,
146 int bottom,
147 int top
148 ) {
149 int d1, d2, d3 ;
150
151 d1 = east - west + 1 ;
152 d2 = north - south + 1 ;
153 d3 = top - bottom + 1 ;
154 #if DEBUG > 0
155 fprintf(stdout, "\n WhichCut : (d1,d2,d3) = (%d,%d,%d)", d1,d2,d3) ;
156 #endif
157 if ( d1 > d2 && d1 > d3 ) {
158 return(1) ; }
159 else if ( d2 > d1 && d2 > d3 ) {
160 return(2) ; }
161 else if ( d3 > d1 && d3 > d2 ) {
162 return(3) ; }
163 else if ( d1 == d2 && d1 > d3 ) {
164 return(1) ; }
165 else if ( d1 == d3 && d1 > d2 ) {
166 return(1) ; }
167 else if ( d2 == d3 && d2 > d1 ) {
168 return(2) ; }
169 else {
170 return(3) ; } }
171
172 /*--------------------------------------------------------------------*/
173 /*
174 ------------------------------------------------
175 purpose -- to split a subregion with a Y-Z plane
176 ------------------------------------------------
177 */
178 static void
SplitX(int n1,int n2,int n3,int newToOld[],int west,int east,int south,int north,int bottom,int top)179 SplitX (
180 int n1,
181 int n2,
182 int n3,
183 int newToOld[],
184 int west,
185 int east,
186 int south,
187 int north,
188 int bottom,
189 int top
190 ) {
191 int depth, east1, east2, height, j, k, k1, k2, k3, m, mid,
192 west1, west2, width1, width2 ;
193
194 # if DEBUG > 0
195 fprintf(stdout, "\n inside SplitX(%d:%d,%d:%d,%d:%d)",
196 west, east, south, north, bottom, top) ;
197 # endif
198 mid = MidPoint(west, east, n1/2) ;
199 west1 = west ;
200 east1 = mid - 1 ;
201 west2 = mid + 2 ;
202 east2 = east ;
203 width1 = east1 - west1 + 1 ;
204 width2 = east2 - west2 + 1 ;
205 height = north - south + 1 ;
206 depth = top - bottom + 1 ;
207 k1 = 0 ;
208 k2 = k1 + width1 * height * depth ;
209 k3 = k2 + width2 * height * depth ;
210 if ( width1 > 0 ) {
211 # if DEBUG > 0
212 fprintf(stdout, "\n SplitX : 1. calling ND(%d,%d,%d,%d,%d,%d",
213 west1, east1, south, north, bottom, top) ;
214 # endif
215 mkNDperm2(n1, n2, n3, newToOld + k1,
216 west1, east1, south, north, bottom, top) ;
217 }
218 if ( width2 > 0 ) {
219 # if DEBUG > 0
220 fprintf(stdout, "\n SplitX : 2. calling ND(%d,%d,%d,%d,%d,%d",
221 west2, east2, south, north, bottom, top) ;
222 # endif
223 mkNDperm2(n1, n2, n3, newToOld + k2,
224 west2, east2, south, north, bottom, top) ;
225 }
226 m = k3 ;
227 for ( k = bottom ; k <= top ; k++ ) {
228 for ( j = south ; j <= north ; j++ ) {
229 # if DEBUG > 0
230 fprintf(stdout, "\n newToOld[%d] = %d",
231 m, mid + j * n1 + k * n1 * n2) ;
232 fprintf(stdout, "\n newToOld[%d] = %d",
233 m+1, mid + 1 + j * n1 + k * n1 * n2) ;
234 # endif
235 newToOld[m++] = mid + j * n1 + k * n1 * n2 ;
236 newToOld[m++] = mid + 1 + j * n1 + k * n1 * n2 ;
237 }
238 }
239 # if DEBUG > 0
240 fprintf(stdout, "\n leaving SplitX(%d:%d,%d:%d,%d:%d)",
241 west, east, south, north, bottom, top) ;
242 # endif
243
244 return ; }
245
246 /*--------------------------------------------------------------------*/
247 /*
248 ------------------------------------------------
249 purpose -- to split a subregion with a X-Z plane
250
251 created -- 95nov16, cca
252 ------------------------------------------------
253 */
254 static void
SplitY(int n1,int n2,int n3,int newToOld[],int west,int east,int south,int north,int bottom,int top)255 SplitY (
256 int n1,
257 int n2,
258 int n3,
259 int newToOld[],
260 int west,
261 int east,
262 int south,
263 int north,
264 int bottom,
265 int top
266 ) {
267 int depth, height1, height2, i, k, k1, k2, k3, m, mid,
268 north1, north2, south1, south2, width ;
269
270 # if DEBUG > 0
271 fprintf(stdout, "\n inside SplitY(%d:%d,%d:%d,%d:%d)",
272 west, east, south, north, bottom, top) ;
273 # endif
274 mid = MidPoint(south, north, n2/2) ;
275 south1 = south ;
276 north1 = mid - 1 ;
277 south2 = mid + 2 ;
278 north2 = north ;
279 # if DEBUG > 0
280 fprintf(stdout,
281 "\n mid = %d, south1 = %d, north1 = %d, south2 = %d, north2 = %d",
282 mid, south1, north1, south2, north2) ;
283 # endif
284 width = east - west + 1 ;
285 height1 = north1 - south1 + 1 ;
286 height2 = north2 - south2 + 1 ;
287 depth = top - bottom + 1 ;
288 k1 = 0 ;
289 k2 = k1 + width * height1 * depth ;
290 k3 = k2 + width * height2 * depth ;
291 if ( height1 > 0 ) {
292 # if DEBUG > 0
293 fprintf(stdout, "\n SplitY : calling ND(%d:%d,%d:%d,%d:%d)",
294 west, east, south1, north1, bottom, top) ;
295 # endif
296 mkNDperm2(n1, n2, n3, newToOld + k1,
297 west, east, south1, north1, bottom, top) ;
298 }
299 if ( height2 > 0 ) {
300 # if DEBUG > 0
301 fprintf(stdout, "\n SplitY : calling ND(%d:%d,%d:%d,%d:%d)",
302 west, east, south2, north2, bottom, top) ;
303 # endif
304 mkNDperm2(n1, n2, n3, newToOld + k2,
305 west, east, south2, north2, bottom, top) ;
306 }
307 m = k3 ;
308 for ( k = bottom ; k <= top ; k++ ) {
309 for ( i = west ; i <= east ; i++ ) {
310 # if DEBUG > 0
311 fprintf(stdout, "\n SplitY : newToOld[%d] = %d",
312 m, i + mid * n1 + k * n1 * n2) ;
313 fprintf(stdout, "\n SplitY : newToOld[%d] = %d",
314 m+1, i + (mid + 1) * n1 + k * n1 * n2) ;
315 # endif
316 newToOld[m++] = i + mid*n1 + k*n1*n2 ;
317 newToOld[m++] = i + (mid + 1)*n1 + k*n1*n2 ;
318 }
319 }
320 # if DEBUG > 0
321 fprintf(stdout, "\n leaving SplitY(%d:%d,%d:%d,%d:%d)",
322 west, east, south, north, bottom, top) ;
323 # endif
324
325 return ; }
326
327 /*--------------------------------------------------------------------*/
328 /*
329 ------------------------------------------------
330 purpose -- to split a subregion with a X-Y plane
331
332 created -- 95nov16, cca
333 ------------------------------------------------
334 */
335 static void
SplitZ(int n1,int n2,int n3,int newToOld[],int west,int east,int south,int north,int bottom,int top)336 SplitZ (
337 int n1,
338 int n2,
339 int n3,
340 int newToOld[],
341 int west,
342 int east,
343 int south,
344 int north,
345 int bottom,
346 int top
347 ) {
348 int bottom1, bottom2, depth1, depth2, height, i, j, k1, k2, k3,
349 m, mid, top1, top2, width ;
350
351 # if DEBUG > 0
352 fprintf(stdout, "\n inside SplitZ(%d:%d,%d:%d,%d:%d)",
353 west, east, south, north, bottom, top) ;
354 # endif
355 mid = MidPoint(bottom, top, n2/2) ;
356 bottom1 = bottom ;
357 top1 = mid - 1 ;
358 bottom2 = mid + 2 ;
359 top2 = top ;
360 width = east - west + 1 ;
361 height = north - south + 1 ;
362 depth1 = top1 - bottom1 + 1 ;
363 depth2 = top2 - bottom2 + 1 ;
364 k1 = 0 ;
365 k2 = k1 + width * height * depth1 ;
366 k3 = k2 + width * height * depth2 ;
367 if ( depth1 > 0 ) {
368 # if DEBUG > 0
369 fprintf(stdout, "\n SplitZ : 1. calling ND(%d:%d,%d:%d,%d:%d)",
370 west, east, south, north, bottom1, top1) ;
371 # endif
372 mkNDperm2(n1, n2, n3, newToOld + k1,
373 west, east, south, north, bottom1, top1) ;
374 }
375 if ( depth2 > 0 ) {
376 # if DEBUG > 0
377 fprintf(stdout, "\n SplitZ : 2. calling ND(%d:%d,%d:%d,%d:%d)",
378 west, east, south, north, bottom2, top2) ;
379 # endif
380 mkNDperm2(n1, n2, n3, newToOld + k2,
381 west, east, south, north, bottom2, top2) ;
382 }
383 m = k3 ;
384 for ( j = south ; j <= north ; j++ ) {
385 for ( i = west ; i <= east ; i++ ) {
386 # if DEBUG > 0
387 fprintf(stdout, "\n newToOld[%d] = %d",
388 m, i + j*n1 + mid*n1*n2) ;
389 fprintf(stdout, "\n newToOld[%d] = %d",
390 m+1, i + j*n1 + (mid+1)*n1*n2) ;
391 # endif
392 newToOld[m++] = i + j*n1 + mid*n1*n2 ;
393 newToOld[m++] = i + j*n1 + (mid + 1)*n1*n2 ;
394 }
395 }
396 # if DEBUG > 0
397 fprintf(stdout, "\n leaving SplitZ(%d:%d,%d:%d,%d:%d)",
398 west, east, south, north, bottom, top) ;
399 # endif
400
401 return ; }
402
403 /*--------------------------------------------------------------------*/
404 /*
405 ----------------------------------------------
406 purpose -- to compute the midpoint of a region
407
408 input --
409
410 left -- left coordinate
411 right -- right coordinate
412 glob_center -- glob_center coordinate
413 ----------------------------------------------
414 */
415 static int
MidPoint(int left,int right,int glob_center)416 MidPoint (
417 int left,
418 int right,
419 int glob_center
420 ) {
421 int mid ;
422
423 /*
424 mid = (left + right)/2 ;
425 if ( (left + right) % 2 == 0 ) {
426 return(mid) ; }
427 else if ( mid >= glob_center ) {
428 return(mid) ; }
429 else {
430 return(mid+1) ; }
431 */
432 mid = left + (right - left - 1)/2 ;
433
434 return(mid) ; }
435
436 /*--------------------------------------------------------------------*/
437