rcb.F [SRC] [CPP] [JOB] [SCAN] SOURCES / FUNCTIONS

```   1 | include(dom.inc)
2 |
3 | !     ===============================================================
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:
 Makefile [SOURCES] - 89 meshpartition.F [SOURCES/MAIN/MASTER] - 78