1 | include(dom.inc)
2 |
3 | ! ===============================================================
4 | ! Copyright (c) CERFACS (all rights reserved)
5 | ! ===============================================================
6 |
7 |
8 | subroutine rcb ( ndim, nnode, npartition, x, weight, sort, &
9 | & rwork, newlist, ioldlist, ibegpart, ilenpart )
10 |
11 | ! !
12 | ! ******************************************************************** !
13 | ! !
14 | ! Purpose: Partitions the graph of a grid using the Recursive !
15 | ! Coordinate Bisection algorithm. !
16 | ! !
17 | ! Input: The number of space dimensions, `ndim'; !
18 | ! the number of nodes `nnode'; the number of !
19 | ! partitions `npartition'; the coordinates of !
20 | ! the nodes of the graph `x'; an optional !
21 | ! real-valued and positive weighting function !
22 | ! associated with each node `weight'; an external !
23 | ! subroutine `sort' for sorting an array ( see !
24 | ! accompanying subroutine library ). !
25 | ! !
26 | ! Output: A list `newlist' whose nth element contains the !
27 | ! new number associated with the nth node ; a list !
28 | ! `ioldlist' that gives the correspondance between !
29 | ! old and new nodes; begin addresses `ibegpart' and !
30 | ! the lengths `ilenpart' of the partitions !
31 | ! of `newlist'. !
32 | ! !
33 | ! Local: A real workspace `rwork' spanning the nodes. !
34 | ! !
35 | ! Authors: M. Rudgyard. !
36 | ! !
37 | ! Checked: 21/4/94 !
38 | ! !
39 | ! Last modified: M. Garcia (Bug ID 152) - 26/05/08 !
40 | ! !
41 | ! Changed to f90(fixed): M. Garcia (10/07/2003) !
42 | ! !
43 | ! ******************************************************************** !
44 | ! !
45 | implicit none
46 |
47 | include 'dom_constants.h'
48 |
49 | ! IN/OUT
50 | ! ------
51 | DOM_INT, intent(in) :: ndim
52 | DOM_INT, intent(in) :: nnode
53 | DOM_INT, intent(inout) :: npartition
54 | DOM_INT, dimension(1:nnode), intent(inout) :: newlist
55 | DOM_INT, dimension(1:nnode), intent(inout) :: ioldlist
56 | DOM_INT, dimension(1:npartition), intent(inout) :: ibegpart
57 | DOM_INT, dimension(1:npartition), intent(inout) :: ilenpart
58 |
59 | DOM_REAL, dimension(1:*), intent(inout) :: weight
60 | DOM_REAL, dimension(1:nnode), intent(inout) :: rwork
61 | DOM_REAL, dimension(1:ndim,1:nnode), intent(in) :: x
62 |
63 | ! LOCAL
64 | ! -----
65 | DOM_INT :: ibp, icount, icut, idiv, idummy, ilp
66 | DOM_INT :: in_list, ipart, ireal_sort, ishift, isize, isort
67 | DOM_INT :: level, n, nd, ndir, new, npart, nlevel, maxcount
68 | DOM_INT, dimension(1:3) :: nmin
69 | DOM_INT, dimension(1:3) :: nmax
70 |
71 | DOM_REAL :: cut, dxmax, totalw, tweight
72 | DOM_REAL, dimension(1:3) :: xmin
73 | DOM_REAL, dimension(1:3) :: xmax
74 |
75 | DOM_INT, dimension(:), allocatable :: idivide
76 |
77 | external sort
78 |
79 | ! !
80 | ! ******************************************************************** !
81 |
82 | !----------------------------------------------------------------------!
83 | ! Define number of partitions required: !
84 | !----------------------------------------------------------------------!
85 |
86 | if (npartition > nnode) then
87 |
88 | write(6,*) '****************** Warning ******************'
89 | write(6,*) 'Unable to cut into required no. of partitions'
90 |
91 | npartition = npartition / 2
92 |
93 | do while (npartition > nnode)
94 |
95 | npartition = npartition / 2
96 |
97 | end do
98 |
99 | write(6,*) 'Number of partitions set to :', npartition
100 |
101 | end if
102 |
103 | !----------------------------------------------------------------------!
104 | ! Initialise arrays: !
105 | !----------------------------------------------------------------------!
106 |
107 | do n = 1, nnode
108 |
109 | newlist(n) = n
110 |
111 | end do
112 |
113 | !----------------------------------------------------------------------!
114 | ! Normalise weighting functions if required: !
115 | !----------------------------------------------------------------------!
116 |
117 | if ( weight(1) > 0.0 ) then
118 |
119 | totalw = 0.0
120 |
121 | do n = 1, nnode
122 |
123 | if (weight(n) < 0.0) then
124 |
125 | write(6,*) '*************** ERROR *****************'
126 | write(6,*) 'Negative weight --- program terminated'
127 | stop
128 |
129 | end if
130 |
131 | totalw = totalw + weight(n)
132 |
133 | end do
134 |
135 | do n = 1, nnode
136 |
137 | weight(n) = weight(n) / totalw
138 |
139 | end do
140 |
141 | end if
142 |
143 | !----------------------------------------------------------------------!
144 | ! Define number of levels required: !
145 | !----------------------------------------------------------------------!
146 |
147 | nlevel = 1
148 |
149 | do while (2**nlevel < npartition)
150 |
151 | nlevel = nlevel + 1
152 |
153 | end do
154 |
155 | !----------------------------------------------------------------------!
156 | ! Allocate local dynamic array: !
157 | !----------------------------------------------------------------------!
158 |
159 | isize = 2**nlevel
160 |
161 | allocate(idivide(isize))
162 |
163 | do n = 1, isize
164 |
165 | idivide(n) = 0
166 |
167 | end do
168 |
169 | !----------------------------------------------------------------------!
170 | ! Set initial level: !
171 | !----------------------------------------------------------------------!
172 |
173 | npart = 1
174 | ibegpart(1) = 1
175 | ilenpart(1) = nnode
176 | idivide(1) = npartition
177 |
178 | !----------------------------------------------------------------------!
179 | ! For each level, loop over existing partitions !
180 | ! and define new partitions: !
181 | !----------------------------------------------------------------------!
182 |
183 | do level = 1, nlevel
184 |
185 | ipart = 1
186 | maxcount = npart
187 |
188 | do icount = 1, maxcount
189 |
190 | idiv = idivide(ipart)
191 | !
192 | ! Do not subdivide unit partitions:
193 | ! ---------------------------------
194 |
195 | if (idivide(ipart) == 1) then
196 |
197 | ipart = ipart + idiv
198 | cycle
199 |
200 | end if
201 |
202 | !
203 | ! Find geometric bounds of the partition:
204 | ! ---------------------------------------
205 |
206 | do nd = 1, ndim
207 |
208 | xmin(nd) = big
209 | xmax(nd) = -big
210 |
211 | do new = ibegpart(ipart), &
212 | & ibegpart(ipart) + ilenpart(ipart) - 1
213 |
214 | n = newlist(new)
215 |
216 | if (x(nd,n) < xmin(nd)) then
217 |
218 | xmin(nd) = x(nd,n)
219 | nmin(nd) = nnode
220 |
221 | end if
222 |
223 | if (x(nd,n) > xmax(nd)) then
224 |
225 | xmax(nd) = x(nd,n)
226 | nmax(nd) = nnode
227 |
228 | end if
229 |
230 | end do
231 |
232 | end do
233 |
234 | !
235 | ! Find the longest coordinate direction:
236 | ! --------------------------------------
237 |
238 | dxmax = -big
239 |
240 | do nd = 1, ndim
241 |
242 | if (xmax(nd)-xmin(nd) > dxmax) then
243 |
244 | dxmax = xmax(nd) - xmin(nd)
245 | ndir = nd
246 |
247 | end if
248 |
249 | end do
250 |
251 | !
252 | ! Fill the nodal workspace to be used for sorting:
253 | ! ------------------------------------------------
254 |
255 | do new = ibegpart(ipart), &
256 | & ibegpart(ipart) + ilenpart(ipart) - 1
257 |
258 | n = newlist(new)
259 | rwork(new) = x(ndir,n)
260 |
261 | end do
262 |
263 | !
264 | ! Sort the graph in this direction:
265 | ! ---------------------------------
266 |
267 | ireal_sort = 1
268 | in_list = 1
269 | isort = 1
270 | ilp = ilenpart(ipart)
271 | ibp = ibegpart(ipart)
272 |
273 | call sort ( ireal_sort, isort, in_list, ilp, &
274 | & idummy, rwork(ibp), newlist(ibp) )
275 |
276 | !
277 | ! Set up new partition:
278 | ! ---------------------
279 |
280 | npart = npart + 1
281 | ishift = idivide(ipart) / 2
282 | idivide(ipart+ishift) = idivide(ipart) - idivide(ipart) / 2
283 |
284 | !
285 | ! Define disection point of old partition:
286 | ! ----------------------------------------
287 |
288 | icut = ( ilenpart(ipart)/idivide(ipart) ) * ishift
289 |
290 | if ( weight(1) > 0.0 ) then
291 |
292 | !
293 | ! Disect graph when the sum of the weighted coordinates exceeds
294 | ! the average value:
295 | ! -------------------------------------------------------------
296 |
297 | totalw = 0.0
298 |
299 | do new = ibegpart(ipart), &
300 | & ibegpart(ipart) + ilenpart(ipart) - 1
301 |
302 | n = newlist(new)
303 | totalw = totalw + weight(n)
304 |
305 | end do
306 |
307 | cut = ( totalw * icut ) / ilenpart(ipart)
308 |
309 | new = ibegpart(ipart) - 1
310 |
311 | tweight = zero
312 | new = new + 1
313 | n = newlist(new)
314 | tweight = tweight + weight(n)
315 |
316 | do while (tweight < cut)
317 |
318 | new = new + 1
319 | n = newlist(new)
320 | tweight = tweight + weight(n)
321 |
322 | end do
323 |
324 | icut = n
325 |
326 | end if
327 |
328 | !
329 | ! Define start addresses and length of new partitions:
330 | ! ----------------------------------------------------
331 |
332 | ibegpart(ipart+ishift) = ibegpart(ipart) + icut
333 | ilenpart(ipart+ishift) = ilenpart(ipart) - icut
334 | ilenpart(ipart) = icut
335 |
336 | idivide(ipart) = ishift
337 |
338 | ipart = ipart + idiv
339 |
340 | end do
341 |
342 | end do
343 |
344 | !----------------------------------------------------------------------!
345 | ! Set up the reverse pointer: !
346 | !----------------------------------------------------------------------!
347 |
348 | do new = 1, nnode
349 |
350 | n = newlist(new)
351 | ioldlist(n) = new
352 |
353 | end do
354 |
355 | !----------------------------------------------------------------------!
356 | ! Deallocate local dynamic array: !
357 | !----------------------------------------------------------------------!
358 |
359 | deallocate(idivide)
360 |
361 |
362 | return
363 | end subroutine rcb
364 |
rcb.F could be called by: