Line data Source code
1 : module spacecurve_mod
2 : use cam_logfile, only: iulog
3 :
4 : implicit none
5 : private
6 :
7 : type, public :: factor_t
8 : integer :: numfact
9 : integer, dimension(:),pointer :: factors => NULL()
10 : end type factor_t
11 :
12 :
13 : integer,public, dimension(:,:), allocatable :: ordered
14 : integer,public, dimension(:,:), allocatable :: dir ! direction to move along each level
15 : integer,public, dimension(:) , allocatable :: pos ! position along each of the axes
16 :
17 : integer,public :: maxdim ! dimensionality of entire space
18 : integer,public :: vcnt ! visitation count
19 : logical,private :: verbose=.FALSE.
20 :
21 : type (factor_t), public :: fact
22 :
23 : SAVE:: fact
24 : public :: map
25 : public :: hilbert_old
26 : public :: PeanoM,hilbert, Cinco
27 : public :: GenCurve
28 : public :: GenSpaceCurve
29 : public :: log2,Factor
30 : public :: PrintCurve
31 : public :: IsFactorable,IsLoadBalanced
32 : public :: genspacepart
33 : contains
34 : !---------------------------------------------------------
35 1536 : recursive function Cinco(l,type,ma,md,ja,jd) result(ierr)
36 :
37 : implicit none
38 : integer,intent(in) :: l,type,ma,md,ja,jd
39 :
40 : integer :: lma,lmd,lja,ljd,ltype
41 : integer :: ll
42 : integer :: ierr
43 : logical :: debug = .FALSE.
44 :
45 1536 : ll = l
46 1536 : if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
47 :
48 : !--------------------------------------------------------------
49 : ! Position [0,0]
50 : !--------------------------------------------------------------
51 1536 : lma = ma
52 1536 : lmd = md
53 1536 : lja = lma
54 1536 : ljd = lmd
55 :
56 1536 : if(ll .gt. 1) then
57 1536 : if(debug) write(iulog,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
58 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
59 1536 : if(debug) call PrintCurve(ordered)
60 : else
61 0 : ierr = IncrementCurve(lja,ljd)
62 0 : if(debug) write(iulog,*)'Cinco: After Position [0,0] ',pos
63 : endif
64 :
65 : !--------------------------------------------------------------
66 : ! Position [1,0]
67 : !--------------------------------------------------------------
68 1536 : lma = ma
69 1536 : lmd = md
70 1536 : lja = lma
71 1536 : ljd = lmd
72 :
73 1536 : if(ll .gt. 1) then
74 1536 : if(debug) write(iulog,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
75 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
76 1536 : if(debug) call PrintCurve(ordered)
77 : else
78 0 : ierr = IncrementCurve(lja,ljd)
79 0 : if(debug) write(iulog,*)'After Position [1,0] ',pos
80 : endif
81 :
82 : !--------------------------------------------------------------
83 : ! Position [2,0]
84 : !--------------------------------------------------------------
85 1536 : lma = MOD(ma+1,maxdim)
86 1536 : lmd = md
87 1536 : lja = lma
88 1536 : ljd = lmd
89 :
90 1536 : if(ll .gt. 1) then
91 1536 : if(debug) write(iulog,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
92 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
93 1536 : if(debug) call PrintCurve(ordered)
94 : else
95 0 : ierr = IncrementCurve(lja,ljd)
96 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
97 : endif
98 :
99 : !--------------------------------------------------------------
100 : ! Position [2,1]
101 : !--------------------------------------------------------------
102 1536 : lma = MOD(ma+1,maxdim)
103 1536 : lmd = md
104 1536 : lja = lma
105 1536 : ljd = lmd
106 :
107 1536 : if(ll .gt. 1) then
108 1536 : if(debug) write(iulog,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
109 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
110 1536 : if(debug) call PrintCurve(ordered)
111 : else
112 0 : ierr = IncrementCurve(lja,ljd)
113 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
114 : endif
115 :
116 : !--------------------------------------------------------------
117 : ! Position [2,2]
118 : !--------------------------------------------------------------
119 1536 : lma = MOD(ma+1,maxdim)
120 1536 : lmd = md
121 1536 : lja = ma
122 1536 : ljd = -md
123 :
124 1536 : if(ll .gt. 1) then
125 1536 : if(debug) write(iulog,25) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
126 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
127 1536 : if(debug) call PrintCurve(ordered)
128 : else
129 0 : ierr = IncrementCurve(lja,ljd)
130 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
131 : endif
132 :
133 :
134 : !--------------------------------------------------------------
135 : ! Position [1,2]
136 : !--------------------------------------------------------------
137 1536 : lma = MOD(ma+1,maxdim)
138 1536 : lmd = -md
139 1536 : lja = lma
140 1536 : ljd = lmd
141 :
142 1536 : if(ll .gt. 1) then
143 1536 : if(debug) write(iulog,26) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
144 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
145 1536 : if(debug) call PrintCurve(ordered)
146 : else
147 0 : ierr = IncrementCurve(lja,ljd)
148 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
149 : endif
150 :
151 : !--------------------------------------------------------------
152 : ! Position [1,1]
153 : !--------------------------------------------------------------
154 1536 : lma = ma
155 1536 : lmd = -md
156 1536 : lja = lma
157 1536 : ljd = lmd
158 :
159 1536 : if(ll .gt. 1) then
160 1536 : if(debug) write(iulog,27) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
161 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
162 1536 : if(debug) call PrintCurve(ordered)
163 : else
164 0 : ierr = IncrementCurve(lja,ljd)
165 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
166 : endif
167 : !--------------------------------------------------------------
168 : ! Position [0,1]
169 : !--------------------------------------------------------------
170 1536 : lma = ma
171 1536 : lmd = -md
172 1536 : lja = MOD(ma+1,maxdim)
173 1536 : ljd = md
174 :
175 1536 : if(ll .gt. 1) then
176 1536 : if(debug) write(iulog,28) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
177 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
178 1536 : if(debug) call PrintCurve(ordered)
179 : else
180 0 : ierr = IncrementCurve(lja,ljd)
181 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
182 : endif
183 :
184 : !--------------------------------------------------------------
185 : ! Position [0,2]
186 : !--------------------------------------------------------------
187 1536 : lma = MOD(ma+1,maxdim)
188 1536 : lmd = md
189 1536 : lja = lma
190 1536 : ljd = lmd
191 :
192 1536 : if(ll .gt. 1) then
193 1536 : if(debug) write(iulog,29) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
194 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
195 1536 : if(debug) call PrintCurve(ordered)
196 : else
197 0 : ierr = IncrementCurve(lja,ljd)
198 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
199 : endif
200 :
201 : !--------------------------------------------------------------
202 : ! Position [0,3]
203 : !--------------------------------------------------------------
204 1536 : lma = MOD(ma+1,maxdim)
205 1536 : lmd = md
206 1536 : lja = lma
207 1536 : ljd = lmd
208 :
209 1536 : if(ll .gt. 1) then
210 1536 : if(debug) write(iulog,30) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
211 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
212 1536 : if(debug) call PrintCurve(ordered)
213 : else
214 0 : ierr = IncrementCurve(lja,ljd)
215 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
216 : endif
217 :
218 : !--------------------------------------------------------------
219 : ! Position [0,4]
220 : !--------------------------------------------------------------
221 1536 : lma = ma
222 1536 : lmd = md
223 1536 : lja = lma
224 1536 : ljd = lmd
225 :
226 1536 : if(ll .gt. 1) then
227 1536 : if(debug) write(iulog,31) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
228 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
229 1536 : if(debug) call PrintCurve(ordered)
230 : else
231 0 : ierr = IncrementCurve(lja,ljd)
232 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
233 : endif
234 :
235 : !--------------------------------------------------------------
236 : ! Position [1,4]
237 : !--------------------------------------------------------------
238 1536 : lma = ma
239 1536 : lmd = md
240 1536 : lja = MOD(ma+1,maxdim)
241 1536 : ljd = -md
242 :
243 1536 : if(ll .gt. 1) then
244 1536 : if(debug) write(iulog,32) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
245 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
246 1536 : if(debug) call PrintCurve(ordered)
247 : else
248 0 : ierr = IncrementCurve(lja,ljd)
249 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
250 : endif
251 :
252 : !--------------------------------------------------------------
253 : ! Position [1,3]
254 : !--------------------------------------------------------------
255 1536 : lma = MOD(ma+1,maxdim)
256 1536 : lmd = -md
257 1536 : lja = ma
258 1536 : ljd = md
259 :
260 1536 : if(ll .gt. 1) then
261 1536 : if(debug) write(iulog,33) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
262 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
263 1536 : if(debug) call PrintCurve(ordered)
264 : else
265 0 : ierr = IncrementCurve(lja,ljd)
266 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
267 : endif
268 :
269 : !--------------------------------------------------------------
270 : ! Position [2,3]
271 : !--------------------------------------------------------------
272 1536 : lma = MOD(ma+1,maxdim)
273 1536 : lmd = md
274 1536 : lja = lma
275 1536 : ljd = lmd
276 :
277 1536 : if(ll .gt. 1) then
278 1536 : if(debug) write(iulog,34) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
279 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
280 1536 : if(debug) call PrintCurve(ordered)
281 : else
282 0 : ierr = IncrementCurve(lja,ljd)
283 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
284 : endif
285 :
286 : !--------------------------------------------------------------
287 : ! Position [2,4]
288 : !--------------------------------------------------------------
289 1536 : lma = ma
290 1536 : lmd = md
291 1536 : lja = lma
292 1536 : ljd = lmd
293 :
294 1536 : if(ll .gt. 1) then
295 1536 : if(debug) write(iulog,35) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
296 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
297 1536 : if(debug) call PrintCurve(ordered)
298 : else
299 0 : ierr = IncrementCurve(lja,ljd)
300 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
301 : endif
302 :
303 : !--------------------------------------------------------------
304 : ! Position [3,4]
305 : !--------------------------------------------------------------
306 1536 : lma = ma
307 1536 : lmd = md
308 1536 : lja = lma
309 1536 : ljd = lmd
310 :
311 1536 : if(ll .gt. 1) then
312 1536 : if(debug) write(iulog,36) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
313 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
314 1536 : if(debug) call PrintCurve(ordered)
315 : else
316 0 : ierr = IncrementCurve(lja,ljd)
317 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
318 : endif
319 :
320 : !--------------------------------------------------------------
321 : ! Position [4,4]
322 : !--------------------------------------------------------------
323 1536 : lma = ma
324 1536 : lmd = md
325 1536 : lja = MOD(ma+1,maxdim)
326 1536 : ljd = -md
327 :
328 1536 : if(ll .gt. 1) then
329 1536 : if(debug) write(iulog,37) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
330 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
331 1536 : if(debug) call PrintCurve(ordered)
332 : else
333 0 : ierr = IncrementCurve(lja,ljd)
334 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
335 : endif
336 :
337 : !--------------------------------------------------------------
338 : ! Position [4,3]
339 : !--------------------------------------------------------------
340 1536 : lma = ma
341 1536 : lmd = -md
342 1536 : lja = lma
343 1536 : ljd = lmd
344 :
345 1536 : if(ll .gt. 1) then
346 1536 : if(debug) write(iulog,38) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
347 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
348 1536 : if(debug) call PrintCurve(ordered)
349 : else
350 0 : ierr = IncrementCurve(lja,ljd)
351 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
352 : endif
353 :
354 : !--------------------------------------------------------------
355 : ! Position [3,3]
356 : !--------------------------------------------------------------
357 1536 : lma = MOD(ma+1,maxdim)
358 1536 : lmd = -md
359 1536 : lja = lma
360 1536 : ljd = lmd
361 :
362 1536 : if(ll .gt. 1) then
363 1536 : if(debug) write(iulog,39) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
364 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
365 1536 : if(debug) call PrintCurve(ordered)
366 : else
367 0 : ierr = IncrementCurve(lja,ljd)
368 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
369 : endif
370 :
371 : !--------------------------------------------------------------
372 : ! Position [3,2]
373 : !--------------------------------------------------------------
374 1536 : lma = MOD(ma+1,maxdim)
375 1536 : lmd = -md
376 1536 : lja = ma
377 1536 : ljd = md
378 :
379 1536 : if(ll .gt. 1) then
380 1536 : if(debug) write(iulog,40) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
381 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
382 1536 : if(debug) call PrintCurve(ordered)
383 : else
384 0 : ierr = IncrementCurve(lja,ljd)
385 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
386 : endif
387 :
388 : !--------------------------------------------------------------
389 : ! Position [4,2]
390 : !--------------------------------------------------------------
391 1536 : lma = ma
392 1536 : lmd = md
393 1536 : lja = MOD(ma+1,maxdim)
394 1536 : ljd = -md
395 :
396 1536 : if(ll .gt. 1) then
397 1536 : if(debug) write(iulog,41) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
398 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
399 1536 : if(debug) call PrintCurve(ordered)
400 : else
401 0 : ierr = IncrementCurve(lja,ljd)
402 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
403 : endif
404 :
405 : !--------------------------------------------------------------
406 : ! Position [4,1]
407 : !--------------------------------------------------------------
408 1536 : lma = ma
409 1536 : lmd = -md
410 1536 : lja = lma
411 1536 : ljd = lmd
412 :
413 1536 : if(ll .gt. 1) then
414 1536 : if(debug) write(iulog,42) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
415 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
416 1536 : if(debug) call PrintCurve(ordered)
417 : else
418 0 : ierr = IncrementCurve(lja,ljd)
419 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
420 : endif
421 :
422 : !--------------------------------------------------------------
423 : ! Position [3,1]
424 : !--------------------------------------------------------------
425 1536 : lma = MOD(ma+1,maxdim)
426 1536 : lmd = -md
427 1536 : lja = lma
428 1536 : ljd = lmd
429 :
430 1536 : if(ll .gt. 1) then
431 1536 : if(debug) write(iulog,43) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
432 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
433 1536 : if(debug) call PrintCurve(ordered)
434 : else
435 0 : ierr = IncrementCurve(lja,ljd)
436 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
437 : endif
438 :
439 : !--------------------------------------------------------------
440 : ! Position [3,0]
441 : !--------------------------------------------------------------
442 1536 : lma = MOD(ma+1,maxdim)
443 1536 : lmd = -md
444 1536 : lja = ma
445 1536 : ljd = md
446 :
447 1536 : if(ll .gt. 1) then
448 1536 : if(debug) write(iulog,44) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
449 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
450 1536 : if(debug) call PrintCurve(ordered)
451 : else
452 0 : ierr = IncrementCurve(lja,ljd)
453 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
454 : endif
455 :
456 : !--------------------------------------------------------------
457 : ! Position [4,0]
458 : !--------------------------------------------------------------
459 1536 : lma = ma
460 1536 : lmd = md
461 1536 : lja = ja
462 1536 : ljd = jd
463 :
464 1536 : if(ll .gt. 1) then
465 1536 : if(debug) write(iulog,45) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
466 1536 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
467 1536 : if(debug) call PrintCurve(ordered)
468 : else
469 0 : ierr = IncrementCurve(lja,ljd)
470 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
471 : endif
472 :
473 : 21 format('Call Cinco Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
474 : 22 format('Call Cinco Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
475 : 23 format('Call Cinco Pos [2,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
476 : 24 format('Call Cinco Pos [2,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
477 : 25 format('Call Cinco Pos [2,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
478 : 26 format('Call Cinco Pos [1,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
479 : 27 format('Call Cinco Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
480 : 28 format('Call Cinco Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
481 : 29 format('Call Cinco Pos [0,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
482 : 30 format('Call Cinco Pos [0,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
483 : 31 format('Call Cinco Pos [0,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
484 : 32 format('Call Cinco Pos [1,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
485 : 33 format('Call Cinco Pos [1,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
486 : 34 format('Call Cinco Pos [2,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
487 : 35 format('Call Cinco Pos [2,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
488 : 36 format('Call Cinco Pos [3,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
489 : 37 format('Call Cinco Pos [4,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
490 : 38 format('Call Cinco Pos [4,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
491 : 39 format('Call Cinco Pos [3,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
492 : 40 format('Call Cinco Pos [3,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
493 : 41 format('Call Cinco Pos [4,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
494 : 42 format('Call Cinco Pos [4,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
495 : 43 format('Call Cinco Pos [3,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
496 : 44 format('Call Cinco Pos [3,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
497 : 45 format('Call Cinco Pos [4,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
498 :
499 1536 : end function Cinco
500 :
501 : !---------------------------------------------------------
502 38400 : recursive function PeanoM(l,type,ma,md,ja,jd) result(ierr)
503 :
504 : implicit none
505 : integer,intent(in) :: l,type,ma,md,ja,jd
506 :
507 : integer :: lma,lmd,lja,ljd,ltype
508 : integer :: ll
509 : integer :: ierr
510 : logical :: debug = .FALSE.
511 :
512 38400 : ll = l
513 38400 : if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
514 : !--------------------------------------------------------------
515 : ! Position [0,0]
516 : !--------------------------------------------------------------
517 38400 : lma = MOD(ma+1,maxdim)
518 38400 : lmd = md
519 38400 : lja = lma
520 38400 : ljd = lmd
521 :
522 38400 : if(ll .gt. 1) then
523 38400 : if(debug) write(iulog,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
524 38400 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
525 38400 : if(debug) call PrintCurve(ordered)
526 : else
527 0 : ierr = IncrementCurve(lja,ljd)
528 0 : if(debug) write(iulog,*)'After Position [0,0] ',pos
529 : endif
530 :
531 :
532 : !--------------------------------------------------------------
533 : ! Position [0,1]
534 : !--------------------------------------------------------------
535 38400 : lma = MOD(ma+1,maxdim)
536 38400 : lmd = md
537 38400 : lja = lma
538 38400 : ljd = lmd
539 38400 : if(ll .gt. 1) then
540 38400 : if(debug) write(iulog,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
541 38400 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
542 38400 : if(debug) call PrintCurve(ordered)
543 : else
544 0 : ierr = IncrementCurve(lja,ljd)
545 0 : if(debug) write(iulog,*)'After Position [0,1] ',pos
546 : endif
547 :
548 : !--------------------------------------------------------------
549 : ! Position [0,2]
550 : !--------------------------------------------------------------
551 38400 : lma = ma
552 38400 : lmd = md
553 38400 : lja = lma
554 38400 : ljd = lmd
555 38400 : if(ll .gt. 1) then
556 38400 : if(debug) write(iulog,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
557 38400 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
558 38400 : if(debug) call PrintCurve(ordered)
559 : else
560 0 : ierr = IncrementCurve(lja,ljd)
561 0 : if(debug) write(iulog,*)'After Position [0,2] ',pos
562 : endif
563 :
564 : !--------------------------------------------------------------
565 : ! Position [1,2]
566 : !--------------------------------------------------------------
567 38400 : lma = ma
568 38400 : lmd = md
569 38400 : lja = lma
570 38400 : ljd = lmd
571 38400 : if(ll .gt. 1) then
572 38400 : if(debug) write(iulog,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
573 38400 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
574 38400 : if(debug) call PrintCurve(ordered)
575 : else
576 0 : ierr = IncrementCurve(lja,ljd)
577 0 : if(debug) write(iulog,*)'After Position [1,2] ',pos
578 : endif
579 :
580 :
581 : !--------------------------------------------------------------
582 : ! Position [2,2]
583 : !--------------------------------------------------------------
584 38400 : lma = ma
585 38400 : lmd = md
586 38400 : lja = MOD(lma+1,maxdim)
587 38400 : ljd = -lmd
588 :
589 38400 : if(ll .gt. 1) then
590 38400 : if(debug) write(iulog,25) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
591 38400 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
592 38400 : if(debug) call PrintCurve(ordered)
593 : else
594 0 : ierr = IncrementCurve(lja,ljd)
595 0 : if(debug) write(iulog,*)'After Position [2,2] ',pos
596 : endif
597 :
598 : !--------------------------------------------------------------
599 : ! Position [2,1]
600 : !--------------------------------------------------------------
601 38400 : lma = ma
602 38400 : lmd = -md
603 38400 : lja = lma
604 38400 : ljd = lmd
605 :
606 38400 : if(ll .gt. 1) then
607 38400 : if(debug) write(iulog,26) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
608 38400 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
609 38400 : if(debug) call PrintCurve(ordered)
610 : else
611 0 : ierr = IncrementCurve(lja,ljd)
612 0 : if(debug) write(iulog,*)'After Position [2,1] ',pos
613 : endif
614 :
615 : !--------------------------------------------------------------
616 : ! Position [1,1]
617 : !--------------------------------------------------------------
618 38400 : lma = MOD(ma+1,maxdim)
619 38400 : lmd = -md
620 38400 : lja = lma
621 38400 : ljd = lmd
622 :
623 38400 : if(ll .gt. 1) then
624 38400 : if(debug) write(iulog,27) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
625 38400 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
626 38400 : if(debug) call PrintCurve(ordered)
627 : else
628 0 : ierr = IncrementCurve(lja,ljd)
629 0 : if(debug) write(iulog,*)'After Position [1,1] ',pos
630 : endif
631 :
632 :
633 : !--------------------------------------------------------------
634 : ! Position [1,0]
635 : !--------------------------------------------------------------
636 38400 : lma = MOD(ma+1,maxdim)
637 38400 : lmd = -md
638 38400 : lja = MOD(lma+1,maxdim)
639 38400 : ljd = -lmd
640 :
641 38400 : if(ll .gt. 1) then
642 38400 : if(debug) write(iulog,28) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
643 38400 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
644 38400 : if(debug) call PrintCurve(ordered)
645 : else
646 0 : ierr = IncrementCurve(lja,ljd)
647 0 : if(debug) write(iulog,*)'After Position [1,0] ',pos
648 : endif
649 :
650 : !--------------------------------------------------------------
651 : ! Position [2,0]
652 : !--------------------------------------------------------------
653 38400 : lma = ma
654 38400 : lmd = md
655 38400 : lja = ja
656 38400 : ljd = jd
657 :
658 38400 : if(ll .gt. 1) then
659 38400 : if(debug) write(iulog,29) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
660 38400 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
661 38400 : if(debug) call PrintCurve(ordered)
662 : else
663 0 : ierr = IncrementCurve(lja,ljd)
664 0 : if(debug) write(iulog,*)'After Position [2,0] ',pos
665 : endif
666 :
667 : 21 format('Call PeanoM Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
668 : 22 format('Call PeanoM Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
669 : 23 format('Call PeanoM Pos [0,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
670 : 24 format('Call PeanoM Pos [1,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
671 : 25 format('Call PeanoM Pos [2,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
672 : 26 format('Call PeanoM Pos [2,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
673 : 27 format('Call PeanoM Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
674 : 28 format('Call PeanoM Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
675 : 29 format('Call PeanoM Pos [2,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
676 :
677 38400 : end function PeanoM
678 : !---------------------------------------------------------
679 345600 : recursive function hilbert(l,type,ma,md,ja,jd) result(ierr)
680 :
681 : implicit none
682 : integer,intent(in) :: l,type,ma,md,ja,jd
683 :
684 : integer :: lma,lmd,lja,ljd,ltype
685 : integer :: ll
686 : integer :: ierr
687 : logical :: debug = .FALSE.
688 :
689 345600 : ll = l
690 345600 : if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
691 : !--------------------------------------------------------------
692 : ! Position [0,0]
693 : !--------------------------------------------------------------
694 345600 : lma = MOD(ma+1,maxdim)
695 345600 : lmd = md
696 345600 : lja = lma
697 345600 : ljd = lmd
698 :
699 345600 : if(ll .gt. 1) then
700 0 : if(debug) write(iulog,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
701 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
702 0 : if(debug) call PrintCurve(ordered)
703 : else
704 345600 : ierr = IncrementCurve(lja,ljd)
705 345600 : if(debug) write(iulog,*)'After Position [0,0] ',pos
706 : endif
707 :
708 :
709 : !--------------------------------------------------------------
710 : ! Position [0,1]
711 : !--------------------------------------------------------------
712 345600 : lma = ma
713 345600 : lmd = md
714 345600 : lja = lma
715 345600 : ljd = lmd
716 345600 : if(ll .gt. 1) then
717 0 : if(debug) write(iulog,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
718 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
719 0 : if(debug) call PrintCurve(ordered)
720 : else
721 345600 : ierr = IncrementCurve(lja,ljd)
722 345600 : if(debug) write(iulog,*)'After Position [0,1] ',pos
723 : endif
724 :
725 :
726 : !--------------------------------------------------------------
727 : ! Position [1,1]
728 : !--------------------------------------------------------------
729 345600 : lma = ma
730 345600 : lmd = md
731 345600 : lja = MOD(ma+1,maxdim)
732 345600 : ljd = -md
733 :
734 345600 : if(ll .gt. 1) then
735 0 : if(debug) write(iulog,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
736 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
737 0 : if(debug) call PrintCurve(ordered)
738 : else
739 345600 : ierr = IncrementCurve(lja,ljd)
740 345600 : if(debug) write(iulog,*)'After Position [1,1] ',pos
741 : endif
742 :
743 : !--------------------------------------------------------------
744 : ! Position [1,0]
745 : !--------------------------------------------------------------
746 345600 : lma = MOD(ma+1,maxdim)
747 345600 : lmd = -md
748 345600 : lja = ja
749 345600 : ljd = jd
750 :
751 345600 : if(ll .gt. 1) then
752 0 : if(debug) write(iulog,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
753 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
754 0 : if(debug) call PrintCurve(ordered)
755 : else
756 345600 : ierr = IncrementCurve(lja,ljd)
757 345600 : if(debug) write(iulog,*)'After Position [1,0] ',pos
758 : endif
759 :
760 : 21 format('Call Hilbert Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
761 : 22 format('Call Hilbert Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
762 : 23 format('Call Hilbert Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
763 : 24 format('Call Hilbert Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
764 :
765 345600 : end function hilbert
766 : !---------------------------------------------------------
767 1382400 : function IncrementCurve(ja,jd) result(ierr)
768 :
769 : implicit none
770 :
771 : integer :: ja,jd
772 : integer :: ierr
773 :
774 1382400 : ordered(pos(0)+1,pos(1)+1) = vcnt
775 1382400 : vcnt = vcnt + 1
776 1382400 : pos(ja) = pos(ja) + jd
777 :
778 1382400 : ierr = 0
779 1382400 : end function IncrementCurve
780 : !---------------------------------------------------------
781 0 : recursive function hilbert_old(l,d,ma,md,ja,jd) result(ierr)
782 :
783 : integer :: l,d ! log base 2 of levels and dimensions left
784 : integer :: ma,md ! main axis and direction
785 : integer :: ja,jd ! joiner axis and direction
786 :
787 : integer :: ierr
788 : integer :: axis
789 : integer :: ll
790 :
791 0 : if(verbose) write(iulog,10) l,d,ma,md,ja,jd,pos(0),pos(1)
792 0 : ll = l ! Copy this to a temporary variable
793 0 : if(d == 0) then
794 0 : ll=ll-1
795 0 : if(ll == 0) then
796 : return
797 : endif
798 0 : axis = ja
799 0 : if(dir(ll,axis) /= jd) then ! do not move away from joiner plane
800 0 : axis = MOD(axis+1,maxdim) ! next axis
801 : endif
802 0 : if(verbose) write(iulog,*)'hilbert_old: call hilbert_old(l,d) #1:'
803 0 : ierr = hilbert_old(ll,maxdim,axis,dir(ll,axis),ja,jd)
804 0 : dir(ll,ja) = -dir(ll,ja)
805 0 : return
806 : endif
807 0 : axis = MOD(ma+1,maxdim)
808 0 : if(verbose) write(iulog,*)'hilbert_old: before call hilbert_old(l,d) #2:'
809 0 : ierr = hilbert_old(ll,d-1,axis,dir(ll,axis),ma,md)
810 0 : if(verbose) write(iulog,*)'hilbert_old: after call hilbert_old(l,d) #2:'
811 0 : if(verbose) write(iulog,30) l,d,ma,md,ja,jd,pos(0),pos(1)
812 :
813 :
814 0 : pos(ma) = pos(ma) + md
815 0 : dir(ll,ma) = - dir(ll,ma)
816 :
817 : !----------------------------------
818 : ! Mark this node as visited
819 : !----------------------------------
820 0 : if(verbose) write(iulog,20) l,d,ma,md,ja,jd,pos(0),pos(1)
821 0 : vcnt=vcnt+1
822 0 : if(verbose) write(iulog,15) pos(0)+1,pos(1)+1,vcnt
823 0 : if(verbose) write(iulog,*)' '
824 0 : if(verbose) write(iulog,*)' '
825 0 : ordered(pos(0)+1,pos(1)+1)=vcnt
826 :
827 0 : if(verbose) write(iulog,*)'hilbert_old: before call hilbert_old(l,d) #3:'
828 0 : ierr = hilbert_old(ll,d-1,axis,dir(ll,axis),ja,jd)
829 0 : if(verbose) write(iulog,*)'hilbert_old: after call hilbert_old(l,d) #3:'
830 :
831 : 10 format('hilbert_old: Entering hilbert_old (l,d,ma,md,ja,jd) are: ', &
832 : 2(i4),' [',2(i3),'][',2(i3),']',2(i3))
833 : 15 format('hilbert_old: mark element {x,y,ordered}:',3(i4))
834 : 20 format('hilbert_old: Before visit code (l,d,ma,md,ja,jd) are:', &
835 : 2(i4),' [',2(i3),'][',2(i3),']',2(i3))
836 :
837 : 30 format('hilbert_old: after call hilbert_old(l,d) #2: (l,d,ma,md,ja,jd are:', &
838 : 2(i4),' [',2(i3),'][',2(i3),']',2(i3))
839 :
840 : end function hilbert_old
841 : !---------------------------------------------------------
842 3072 : function log2( n)
843 :
844 : implicit none
845 :
846 : integer :: n
847 :
848 : integer :: log2,tmp
849 : !
850 : ! Find the log2 of input value
851 : !
852 3072 : log2 = 1
853 3072 : tmp =n
854 12288 : do while (tmp/2 .ne. 1)
855 9216 : tmp=tmp/2
856 9216 : log2=log2+1
857 : enddo
858 :
859 3072 : end function log2
860 : !---------------------------------------------------------
861 0 : function IsLoadBalanced(nelem,npart)
862 :
863 : implicit none
864 :
865 : integer :: nelem,npart
866 :
867 : logical :: IsLoadBalanced
868 :
869 : integer :: tmp1
870 :
871 0 : tmp1 = nelem/npart
872 :
873 0 : if(npart*tmp1 == nelem ) then
874 : IsLoadBalanced=.TRUE.
875 : else
876 0 : IsLoadBalanced=.FALSE.
877 : endif
878 :
879 0 : end function IsLoadBalanced
880 : !---------------------------------------------------------
881 385536 : recursive function GenCurve(l,type,ma,md,ja,jd) result(ierr)
882 :
883 : implicit none
884 : integer,intent(in) :: l,type,ma,md,ja,jd
885 : integer :: ierr
886 :
887 385536 : if(type == 2) then
888 345600 : ierr = hilbert(l,type,ma,md,ja,jd)
889 39936 : elseif ( type == 3) then
890 38400 : ierr = PeanoM(l,type,ma,md,ja,jd)
891 1536 : elseif ( type == 5) then
892 1536 : ierr = Cinco(l,type,ma,md,ja,jd)
893 : endif
894 :
895 385536 : end function GenCurve
896 : !---------------------------------------------------------
897 3072 : function Factor(num) result(res)
898 :
899 : implicit none
900 : integer,intent(in) :: num
901 :
902 : type (factor_t) :: res
903 : integer :: tmp,tmp2,tmp3,tmp5
904 : integer :: i,n
905 : logical :: found
906 :
907 : ! --------------------------------------
908 : ! Allocate for max # of factors
909 : ! --------------------------------------
910 3072 : tmp = num
911 3072 : tmp2 = log2(num)
912 9216 : allocate(res%factors(tmp2))
913 :
914 3072 : n=0
915 : !-----------------------
916 : ! Look for factors of 2
917 : !-----------------------
918 3072 : found=.TRUE.
919 9216 : do while (found)
920 6144 : found = .FALSE.
921 6144 : tmp2 = tmp/2
922 9216 : if( tmp2*2 == tmp ) then
923 3072 : n = n + 1
924 3072 : res%factors(n) = 2
925 3072 : found = .TRUE.
926 3072 : tmp = tmp2
927 : endif
928 : enddo
929 :
930 : !-----------------------
931 : ! Look for factors of 3
932 : !-----------------------
933 : found=.TRUE.
934 9216 : do while (found)
935 6144 : found = .FALSE.
936 6144 : tmp3 = tmp/3
937 9216 : if( tmp3*3 == tmp ) then
938 3072 : n = n + 1
939 3072 : res%factors(n) = 3
940 3072 : found = .TRUE.
941 3072 : tmp = tmp3
942 : endif
943 : enddo
944 :
945 : !-----------------------
946 : ! Look for factors of 5
947 : !-----------------------
948 : found=.TRUE.
949 9216 : do while (found)
950 6144 : found = .FALSE.
951 6144 : tmp5 = tmp/5
952 9216 : if( tmp5*5 == tmp ) then
953 3072 : n = n + 1
954 3072 : res%factors(n) = 5
955 3072 : found = .TRUE.
956 3072 : tmp = tmp5
957 : endif
958 : enddo
959 :
960 : tmp=1
961 12288 : do i=1,n
962 12288 : tmp = tmp * res%factors(i)
963 : enddo
964 3072 : if(tmp == num) then
965 3072 : res%numfact = n
966 : else
967 0 : res%numfact = -1
968 : endif
969 :
970 3072 : end function Factor
971 : !---------------------------------------------------------
972 :
973 1536 : function IsFactorable(n)
974 : use cam_abortutils, only: endrun
975 :
976 : integer,intent(in) :: n
977 : type (factor_t) :: fact
978 :
979 : logical :: IsFactorable
980 :
981 : if (associated(fact%factors)) then
982 : call endrun("fact already allocated!!!")
983 : end if
984 1536 : fact = Factor(n)
985 1536 : if(fact%numfact .ne. -1) then
986 : IsFactorable = .TRUE.
987 : else
988 0 : IsFactorable = .FALSE.
989 : endif
990 :
991 1536 : end function IsFactorable
992 : !------------------------------------------------
993 :
994 1536 : subroutine map(l)
995 :
996 : implicit none
997 : integer :: l,d
998 : integer :: type, ierr
999 :
1000 1536 : d = SIZE(pos)
1001 :
1002 4608 : pos=0
1003 1536 : maxdim=d
1004 1536 : vcnt=0
1005 :
1006 1536 : type = fact%factors(l)
1007 1536 : ierr = GenCurve(l,type,0,1,0,1)
1008 :
1009 1536 : end subroutine map
1010 : !---------------------------------------------------------
1011 1536 : subroutine GenSpaceCurve(Mesh)
1012 :
1013 : implicit none
1014 :
1015 : integer,target,intent(inout) :: Mesh(:,:)
1016 : integer :: level,dim
1017 :
1018 : integer :: gridsize
1019 :
1020 : ! Setup the size of the grid to traverse
1021 :
1022 1536 : dim = 2
1023 1536 : gridsize = SIZE(Mesh,dim=1)
1024 1536 : fact = factor(gridsize)
1025 1536 : level = fact%numfact
1026 :
1027 1536 : if(verbose) write(iulog,*)'GenSpacecurve: level is ',level
1028 6144 : allocate(ordered(gridsize,gridsize))
1029 :
1030 : ! Setup the working arrays for the traversal
1031 1536 : allocate(pos(0:dim-1))
1032 :
1033 : ! The array ordered will contain the visitation order
1034 1430016 : ordered(:,:) = 0
1035 :
1036 1536 : call map(level)
1037 :
1038 1430016 : Mesh(:,:) = ordered(:,:)
1039 :
1040 1536 : end subroutine GenSpaceCurve
1041 : !-------------------------------------------------------------------------------------------------------
1042 0 : subroutine PrintCurve(Mesh)
1043 : implicit none
1044 : integer,target :: Mesh(:,:)
1045 : integer :: gridsize,i
1046 :
1047 0 : gridsize = SIZE(Mesh,dim=1)
1048 :
1049 0 : if(gridsize == 2) then
1050 0 : write (iulog,*) "A Level 1 Hilbert Curve:"
1051 0 : write (iulog,*) "------------------------"
1052 0 : do i=1,gridsize
1053 0 : write(iulog,2) Mesh(1,i),Mesh(2,i)
1054 : enddo
1055 : else if(gridsize == 3) then
1056 0 : write (iulog,*) "A Level 1 Peano Meandering Curve:"
1057 0 : write (iulog,*) "---------------------------------"
1058 0 : do i=1,gridsize
1059 0 : write(iulog,3) Mesh(1,i),Mesh(2,i),Mesh(3,i)
1060 : enddo
1061 : else if(gridsize == 4) then
1062 0 : write (iulog,*) "A Level 2 Hilbert Curve:"
1063 0 : write (iulog,*) "------------------------"
1064 0 : do i=1,gridsize
1065 0 : write(iulog,4) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i)
1066 : enddo
1067 : else if(gridsize == 5) then
1068 0 : write (iulog,*) "A Level 1 Cinco Curve:"
1069 0 : write (iulog,*) "------------------------"
1070 0 : do i=1,gridsize
1071 0 : write(iulog,5) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i),Mesh(5,i)
1072 : enddo
1073 : else if(gridsize == 6) then
1074 0 : write (iulog,*) "A Level 1 Hilbert and Level 1 Peano Curve:"
1075 0 : write (iulog,*) "------------------------------------------"
1076 0 : do i=1,gridsize
1077 0 : write(iulog,6) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i),Mesh(5,i),Mesh(6,i)
1078 : enddo
1079 : else if(gridsize == 8) then
1080 0 : write (iulog,*) "A Level 3 Hilbert Curve:"
1081 0 : write (iulog,*) "------------------------"
1082 0 : do i=1,gridsize
1083 0 : write(iulog,8) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1084 0 : Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i)
1085 : enddo
1086 : else if(gridsize == 9) then
1087 0 : write (iulog,*) "A Level 2 Peano Meandering Curve:"
1088 0 : write (iulog,*) "---------------------------------"
1089 0 : do i=1,gridsize
1090 0 : write(iulog,9) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1091 0 : Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
1092 0 : Mesh(9,i)
1093 : enddo
1094 : else if(gridsize == 10) then
1095 0 : write (iulog,*) "A Level 1 Hilbert and Level 1 Cinco Curve:"
1096 0 : write (iulog,*) "---------------------------------"
1097 0 : do i=1,gridsize
1098 0 : write(iulog,10) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1099 0 : Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
1100 0 : Mesh(9,i),Mesh(10,i)
1101 : enddo
1102 : else if(gridsize == 12) then
1103 0 : write (iulog,*) "A Level 2 Hilbert and Level 1 Peano Curve:"
1104 0 : write (iulog,*) "------------------------------------------"
1105 0 : do i=1,gridsize
1106 0 : write(iulog,12) Mesh(1,i),Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1107 0 : Mesh(5,i),Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1108 0 : Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i)
1109 : enddo
1110 : else if(gridsize == 15) then
1111 0 : write (iulog,*) "A Level 1 Peano and Level 1 Cinco Curve:"
1112 0 : write (iulog,*) "------------------------"
1113 0 : do i=1,gridsize
1114 0 : write(iulog,15) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1115 0 : Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
1116 0 : Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1117 0 : Mesh(13,i),Mesh(14,i),Mesh(15,i)
1118 : enddo
1119 : else if(gridsize == 16) then
1120 0 : write (iulog,*) "A Level 4 Hilbert Curve:"
1121 0 : write (iulog,*) "------------------------"
1122 0 : do i=1,gridsize
1123 0 : write(iulog,16) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1124 0 : Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
1125 0 : Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1126 0 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i)
1127 : enddo
1128 : else if(gridsize == 18) then
1129 0 : write (iulog,*) "A Level 1 Hilbert and Level 2 Peano Curve:"
1130 0 : write (iulog,*) "------------------------------------------"
1131 0 : do i=1,gridsize
1132 0 : write(iulog,18) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1133 0 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1134 0 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1135 0 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1136 0 : Mesh(17,i),Mesh(18,i)
1137 : enddo
1138 : else if(gridsize == 20) then
1139 0 : write (iulog,*) "A Level 2 Hilbert and Level 1 Cinco Curve:"
1140 0 : write (iulog,*) "------------------------------------------"
1141 0 : do i=1,gridsize
1142 0 : write(iulog,20) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1143 0 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1144 0 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1145 0 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1146 0 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i)
1147 : enddo
1148 : else if(gridsize == 24) then
1149 0 : write (iulog,*) "A Level 3 Hilbert and Level 1 Peano Curve:"
1150 0 : write (iulog,*) "------------------------------------------"
1151 0 : do i=1,gridsize
1152 0 : write(iulog,24) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1153 0 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1154 0 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1155 0 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1156 0 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
1157 0 : Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i)
1158 : enddo
1159 : else if(gridsize == 25) then
1160 0 : write (iulog,*) "A Level 2 Cinco Curve:"
1161 0 : write (iulog,*) "------------------------------------------"
1162 0 : do i=1,gridsize
1163 0 : write(iulog,25) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1164 0 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1165 0 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1166 0 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1167 0 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
1168 0 : Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
1169 0 : Mesh(25,i)
1170 : enddo
1171 : else if(gridsize == 27) then
1172 0 : write (iulog,*) "A Level 3 Peano Meandering Curve:"
1173 0 : write (iulog,*) "---------------------------------"
1174 0 : do i=1,gridsize
1175 0 : write(iulog,27) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1176 0 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1177 0 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1178 0 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1179 0 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
1180 0 : Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
1181 0 : Mesh(25,i),Mesh(26,i),Mesh(27,i)
1182 : enddo
1183 : else if(gridsize == 32) then
1184 0 : write (iulog,*) "A Level 5 Hilbert Curve:"
1185 0 : write (iulog,*) "------------------------"
1186 0 : do i=1,gridsize
1187 0 : write(iulog,32) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1188 0 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
1189 0 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
1190 0 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
1191 0 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
1192 0 : Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
1193 0 : Mesh(25,i),Mesh(26,i),Mesh(27,i),Mesh(28,i), &
1194 0 : Mesh(29,i),Mesh(30,i),Mesh(31,i),Mesh(32,i)
1195 : enddo
1196 : endif
1197 : 2 format('|',2(i2,'|'))
1198 : 3 format('|',3(i2,'|'))
1199 : 4 format('|',4(i2,'|'))
1200 : 5 format('|',5(i2,'|'))
1201 : 6 format('|',6(i2,'|'))
1202 : 8 format('|',8(i2,'|'))
1203 : 9 format('|',9(i2,'|'))
1204 : 10 format('|',10(i2,'|'))
1205 : 12 format('|',12(i3,'|'))
1206 : 15 format('|',15(i3,'|'))
1207 : 16 format('|',16(i3,'|'))
1208 : 18 format('|',18(i3,'|'))
1209 : 20 format('|',20(i3,'|'))
1210 : 24 format('|',24(i3,'|'))
1211 : 25 format('|',25(i3,'|'))
1212 : 27 format('|',27(i3,'|'))
1213 : 32 format('|',32(i4,'|'))
1214 :
1215 0 : end subroutine PrintCurve
1216 :
1217 : !-------------------------------------------------------------------------------------------------------
1218 1536 : subroutine genspacepart(GridVertex)
1219 : use dimensions_mod, only: npart
1220 : use gridgraph_mod, only: gridedge_t, gridvertex_t
1221 :
1222 : type (GridVertex_t), intent(inout) :: GridVertex(:)
1223 :
1224 : integer :: nelem, nelemd
1225 : integer :: k, tmp1, id, s1, extra
1226 :
1227 1536 : nelem = SIZE(GridVertex(:))
1228 :
1229 1536 : nelemd = nelem / npart
1230 : ! every cpu gets nelemd elements, but the first 'extra' get nelemd+1
1231 1536 : extra = mod(nelem,npart)
1232 1536 : s1 = extra*(nelemd+1)
1233 :
1234 : ! split curve into two curves:
1235 : ! 1 ... s1 s2 ... nelem
1236 : !
1237 : ! s1 = extra*(nelemd+1) (count be 0)
1238 : ! s2 = s1+1
1239 : !
1240 : ! First region gets nelemd+1 elements per Processor
1241 : ! Second region gets nelemd elements per Processor
1242 :
1243 : ! ===========================================
1244 : ! Add the partitioning information into the
1245 : ! Grid Vertex and Grid Edge structures
1246 : ! ===========================================
1247 :
1248 8295936 : do k = 1, nelem
1249 8294400 : id = GridVertex(k)%SpaceCurve
1250 8295936 : if (id <= s1) then
1251 296448 : tmp1 = id/(nelemd+1)
1252 296448 : GridVertex(k)%processor_number = tmp1 + 1
1253 : else
1254 7997952 : id = id - s1
1255 7997952 : tmp1 = id / nelemd
1256 7997952 : GridVertex(k)%processor_number = extra + tmp1+1
1257 : end if
1258 : end do
1259 : #if 0
1260 : if (masterproc) then
1261 : write(iulog, *)'Space-Filling Curve Parititioning: '
1262 : write(iulog, '(2(a,i0))') 'npart = ',npart,', nelem = ',nelem
1263 : write(iulog, '(2(a,i0))') 'nelemd = ',npart,', extra = ',extra
1264 : write(iulog, '(a)') ' elem task#'
1265 : do k = 1, nelem
1266 : write(iulog,'(i6," ",i6)') k, GridVertex(k)%processor_number
1267 : end do
1268 : end if
1269 : call mpi_barrier(mpicom, tmp1)
1270 : #endif
1271 :
1272 1536 : end subroutine genspacepart
1273 :
1274 0 : end module spacecurve_mod
|
